ATLAS Offline Software
ErrorMatrixCompressor.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 /*
6  * ErrorMatrixCompressor.cxx
7  *
8  * Created by Dmitry Emeliyanov on 12/11/08.
9  * <Dmitry.Emeliyanov@cern.ch>
10  *
11  */
12 
13 #include <stdio.h>
14 #include <stdlib.h>
15 #include <math.h>
16 
17 #include "EventCommonTPCnv/ErrorMatrixCompressor.h"
18 
19 #include <fstream>
20 #include <cstring>
21 
22 using std::memset;
23 
24 void DecoderFloat_IEEE754::print()
25 {
26  // printf("0x%X\n",m_data.l);
27  unsigned int mask = 0x00000001;
28  for(int i=31;i>=0;i--)
29  {
30  unsigned int m = mask << i;
31  //printf("0x%X\n",m);
32  if(m_data.l & m) std::cout<<"1";
33  else std::cout<<"0";
34  }
35  std::cout<<std::endl;
36  /*
37  mask = 0x80000000;
38  int sign = (m_data.l & mask) >> 31;
39  mask = 0x7f800000;
40  short int exponent = (m_data.l & mask) >> 23;
41  mask = 0x007fffff;
42  unsigned int mant = (m_data.l & mask);
43  printf("Sign=0x%X Exponent = 0x%X Mantissa = 0x%X\n",sign,exponent,mant);
44  printf("Sign=%d Exponent = %d Mantissa = %ld\n",sign,exponent,mant);
45  */
46 }
47 
48 short int DecoderFloat_IEEE754::getExponent()
49 {
50  unsigned int mask = 0x7f800000;
51  return ((m_data.l & mask) >> 23);
52 }
53 
54 unsigned int DecoderFloat_IEEE754::getMantissa()
55 {
56  unsigned int mask = 0x007fffff;
57  return (m_data.l & mask);
58 }
59 
60 
61 void DecoderFloat_IEEE754::setExponent(short int ex)
62 {
63  unsigned int buf=ex;
64  buf=buf<<23;
65  unsigned int mask=0x7F800000;
66  m_data.l=m_data.l | (buf & mask);
67 }
68 
69 void DecoderFloat_IEEE754::setSign(int s)
70 {
71  if(s>0)
72  {
73  m_data.l=m_data.l & 0x7FFFFFFF;
74  }
75  else
76  {
77  m_data.l=m_data.l | 0x80000000;
78  }
79 }
80 
81 
82 ErrorMatrixCompressor::ErrorMatrixCompressor(int b)
83 {
84  short int biasArray[6] = {116,120,118,116,116,102};
85 
86 
87  double T[5]={0.2,0.5,100.0,100.0,20000.0};
88 
89  m_bitStrip=b;
90  for(int i=0;i<6;i++)
91  {
92  m_biases[i]=biasArray[i];
93  }
94  for(int i=0;i<5;i++)
95  {
96  m_scales[i]=T[i];
97  }
98  m_limits[0]=15;
99  m_limits[1]=31;
100  m_srcMasks[0]= 0x80000000;
101  m_srcMasks[1]= 0xC0000000;
102  m_srcMasks[2]= 0xE0000000;
103  m_srcMasks[3]= 0xF0000000;
104  m_srcMasks[4]= 0xF8000000;
105  m_srcMasks[5]= 0xFC000000;
106  m_srcMasks[6]= 0xFE000000;
107  m_srcMasks[7]= 0xFF000000;
108  m_srcMasks[8]= 0xFF800000;
109  m_srcMasks[9]= 0xFFC00000;
110  m_srcMasks[10]=0xFFE00000;
111  m_srcMasks[11]=0xFFF00000;
112  m_srcMasks[12]=0xFFF80000;
113  m_srcMasks[13]=0xFFFC0000;
114  m_srcMasks[14]=0xFFFE0000;
115  m_srcMasks[15]=0xFFFF0000;
116  m_tripVec.clear();
117  m_tripVec.push_back(Triplet(0,1,3));
118  m_tripVec.push_back(Triplet(2,4,6));
119  m_tripVec.push_back(Triplet(5,7,8));
120  m_tripVec.push_back(Triplet(9,10,11));
121  m_tripVec.push_back(Triplet(14,12,13));
122 }
123 
124 void ErrorMatrixCompressor::setBiases(const int b[6])
125 {
126  for(int i=0;i<6;i++)
127  {
128  m_biases[i]=b[i];
129  }
130 }
131 
132 void ErrorMatrixCompressor::setUpperLimits(const int l[2])
133 {
134  for(int i=0;i<2;i++)
135  {
136  m_limits[i]=l[i];
137  }
138 }
139 
140 
141 bool ErrorMatrixCompressor::CholeskyDecomposition(double a[5][5], double L[5][5])
142 {
143 
144  int i,j,k;
145  double sum;
146  double p[5];
147 
148  for(i=0;i<5;i++)
149  {
150  for(j=i;j<5;j++)
151  {
152  sum=a[i][j];
153  for(k=i-1;k>=0;k--)
154  sum-=a[i][k]*a[j][k];
155  if(i==j)
156  {
157  if(sum<=0.0)
158  {
159  return false;
160  }
161  p[i]=sqrt(sum);L[i][i]=p[i];
162  }
163  else
164  {
165  a[j][i]=sum/p[i];
166  L[j][i]=a[j][i];
167  }
168  }
169  }
170  return true;
171 }
172 
173 bool ErrorMatrixCompressor::compress(const std::vector<double>& src, std::vector<unsigned int>& dest)
174 {
175  int i,j;
176  double L[5][5],C0[5][5],C[5][5];
177  float S[5][5];
178 
179  dest.clear();
180 
181  int idx=0;
182  for (i=0;i<5;i++)
183  for(j=0;j<=i;j++)
184  {
185  C[i][j]=src[idx];
186  C[j][i]=C0[i][j]=C0[j][i]=C[i][j];
187  idx++;
188  }
189 
190  memset(&L[0][0],0,sizeof(L));
191 
192  for(i=0;i<5;i++)
193  for(j=0;j<=i;j++)
194  {
195  C[i][j]*=(m_scales[i]*m_scales[j]);
196  C[j][i]=C[i][j];
197  }
198  if(!CholeskyDecomposition(C,L)) return false;
199  memset(&S[0][0],0,sizeof(S));
200  for(i=0;i<5;i++)
201  for(j=0;j<=i;j++) S[i][j]=L[i][j];
202 
203  std::vector<FloatRep> vecFR;
204  vecFR.clear();
205 
206  for(i=0;i<5;i++)
207  {
208  for(j=0;j<=i;j++)
209  {
210  m_decoder.setF(S[i][j]);
211  char sign=(S[i][j]<0)?1:0;
212  unsigned int mant=m_decoder.getMantissa();
213  unsigned short int ex=m_decoder.getExponent();
214  vecFR.push_back(FloatRep(sign,ex,mant));
215  }
216  }
217  std::vector<unsigned short> vShorts;
218  vShorts.clear();
219  if(!compressFR(vecFR,vShorts)) return false;
220 
221  //re-packing as ints
222 
223  // std::cout<<"Total "<<vShorts.size()<<" unsigned short"<<std::endl;
224 
225  int nBits=0;
226  unsigned int buffer = 0x00000000;
227 
228  for(std::vector<unsigned short>::iterator it = vShorts.begin(); it != vShorts.end();++it) {
229  // printf("short = 0x%X\n",(*it));
230  if(nBits==2) {
231  dest.push_back(buffer);
232  nBits=0;
233  buffer = 0x00000000;
234  }
235  if(nBits==0) {
236  buffer = (*it);
237  buffer = buffer << 16;
238  nBits++;
239  }
240  else {
241  buffer = buffer | (*it);
242  nBits++;
243  }
244  }
245  dest.push_back(buffer);
246  return true;
247 }
248 
249 bool ErrorMatrixCompressor::restore(const std::vector<unsigned int>& src, std::vector<double>& dest)
250 {
251  int i,j;
252  float S[5][5];
253  double L[5][5],C[5][5];
254 
255  dest.clear();
256 
257  std::vector<FloatRep> vecFR;
258  vecFR.clear();
259 
260  std::vector<unsigned short> vShorts;
261 
262  vShorts.clear();
263 
264  for (unsigned int ii : src) {
265  unsigned short s1,s2;
266 
267  s1 = (unsigned short)((0xFFFF0000 & ii) >> 16);
268  s2 = (unsigned short)(0x0000FFFF & ii);
269  vShorts.push_back(s1);
270  //if(!((s2==0) && ((it+1)==src.end()))) //do not store last zero
271  vShorts.push_back(s2);
272  }
273 
274  if(!restoreFR(vShorts,vecFR)) return false;
275 
276  std::vector<FloatRep>::iterator fIt(vecFR.begin());
277  memset(&S[0][0],0,sizeof(S));
278  for(i=0;i<5;i++)
279  for(j=0;j<=i;j++)
280  {
281  if(fIt==vecFR.end()) break;
282  S[i][j]=(*fIt).restore();
283  ++fIt;
284  }
285 
286  memset(&L[0][0],0,sizeof(L));
287  for(i=0;i<5;i++)
288  for(j=0;j<=i;j++) L[i][j]=S[i][j];
289  for(i=0;i<5;i++)
290  for(j=i;j<5;j++)
291  {
292  C[i][j]=0.0;
293  for(int k=0;k<5;k++)
294  C[i][j]+=L[i][k]*L[j][k];
295  C[j][i]=C[i][j];
296  }
297  for(i=0;i<5;i++)
298  for(j=0;j<=i;j++)
299  {
300  C[i][j]/=(m_scales[i]*m_scales[j]);
301  C[j][i]=C[i][j];
302  }
303  for(i=0;i<5;i++) for(j=0;j<=i;j++) dest.push_back(C[i][j]);
304  return true;
305 }
306 
307 
308 bool ErrorMatrixCompressor::compressFR(const std::vector<FloatRep>& src, std::vector<unsigned short>& dest)
309 {
310  unsigned short buf=0x0000;
311  dest.clear();
312  int nMantLength=23-m_bitStrip;
313  if(nMantLength<8)
314  {
315  std::cout<<"Requested mantissa reduction is too large: 23->"<<nMantLength<<std::endl;
316  return false;
317  }
318  std::vector<FloatRep>::const_iterator fIt;
319  /*
320  for(fIt=src.begin();fIt!=src.end();++fIt)
321  {
322  printf("%d %d %d -> %2.8f\n",(*fIt).sign(),
323  (*fIt).exponent(),(*fIt).mantissa(),(*fIt).restore());
324  }
325  */
326  // 1. Check limits
327  fIt=src.begin();
328  int i,j;
329  for (i=0;i<5;i++)
330  for (j=0;j<=i;j++)
331  {
332  unsigned short int ex=(*fIt).exponent();
333  unsigned int mant=(*fIt).mantissa();
334  int bias,limit;
335  if(i==j)
336  {
337  bias=m_biases[i];
338  limit=m_limits[0];
339  }
340  else
341  {
342  bias=m_biases[5];
343  limit=m_limits[1];
344  }
345  if((ex!=0)&&(mant!=0))
346  {
347  if((ex<bias)||(ex-bias>limit))
348  {
349  // std::cout<<"i="<<i<<" j="<<j<<" ex="<<ex<<" bias="<<bias<<" lim="<<limit<<std::endl;
350  return false;
351  }
352  }
353  ++fIt;
354  }
355  // 2. Pack exponents: diag e + non-diag s1 + non-diag e1 + non-diag s2 + non-diag e2
356 
357  i=0;
358  for(std::vector<Triplet>::iterator trIt=m_tripVec.begin();trIt!=m_tripVec.end();++trIt)
359  {
360  int i1,i2[2];
361  i1=(*trIt).m_d;
362  i2[0]=(*trIt).m_nd1;i2[1]=(*trIt).m_nd2;
363  //printf("%d %d %d\n",i1,i2[0],i2[1]);
364  buf=0x0000;
365  unsigned short e;
366 
367  e=src[i1].exponent()-m_biases[i];
368  buf = buf | ((e<<12) & 0xF000);
369  if((src[i2[0]].exponent()==0)&&(src[i2[0]].mantissa()==0))
370  {
371  e=src[i2[0]].exponent();
372  buf = buf | 0x0800;//-0
373  }
374  else
375  {
376  if(src[i2[0]].sign())
377  buf = buf | 0x0800;
378  e=src[i2[0]].exponent()-m_biases[5];
379  }
380  buf = buf | ((e<<6) & 0x07C0);
381  if((src[i2[1]].exponent()==0)&&(src[i2[1]].mantissa()==0))
382  {
383  buf = buf | 0x0020;//-0
384  e=src[i2[1]].exponent();
385  }
386  else
387  {
388  e=src[i2[1]].exponent()-m_biases[5];
389  if(src[i2[1]].sign())
390  buf = buf | 0x0020;
391  }
392  buf = buf | (e & 0x001F);
393  // printf("Exponent 0x%X\n",buf);
394  dest.push_back(buf);
395  i++;
396  }
397  // 3. Pack reduced mantissas
398  fIt=src.begin();
399  unsigned int nPacked=0;
400  int nFreeBits=0,nBitsToStore=0,nBufferLength=0;
401  unsigned int srcBuffer=0x00000000;
402  // printf("L=%d\n",nMantLength);
403  while (nPacked<=src.size()+1)
404  {
405  if(nFreeBits==0)
406  {
407  if(nBufferLength!=0)
408  {
409  dest.push_back(buf);
410  // printf("Storing 0x%X\n",buf);
411  }
412  buf=0x0000;
413  nFreeBits=16;
414  nBufferLength++;
415  }
416  if(nBitsToStore==0)
417  {
418  if(nPacked!=0)
419  ++fIt;
420  if(fIt==src.end()) break;
421  nPacked++;
422  srcBuffer=((*fIt).mantissa()<<9);
423  // printf("Packing 0x%X = %ud\n",(*fIt).mantissa(),(*fIt).mantissa());
424  nBitsToStore=nMantLength;
425  }
426  int Np=(nBitsToStore>nFreeBits)?nFreeBits:nBitsToStore;
427  unsigned int mask=m_srcMasks[Np-1];
428  unsigned int slice = srcBuffer & mask;
429  slice = (slice >> (32-nFreeBits)) & 0x0000FFFF;
430  unsigned int tmp=(unsigned int)(slice);
431  // printf("Np=%d Tmp=0x%X\n",Np,tmp);
432  buf = buf | tmp;
433  srcBuffer = srcBuffer << Np;
434  // printf("dest=0x%X src=0x%X\n",buf,srcBuffer);
435  nFreeBits-=Np;
436  nBitsToStore-=Np;
437  //cout<<"F="<<nFreeBits<<" TS="<<nBitsToStore<<endl;
438  }
439  // printf("Storing 0x%X\n",buf);
440  dest.push_back(buf);
441  return true;
442 }
443 
444 bool ErrorMatrixCompressor::restoreFR(const std::vector<unsigned short>& src, std::vector<FloatRep>& dest)
445 {
446  int i,nRestored,nFreeBits,nBitsToStore;
447  unsigned short buf=0x0000;
448  unsigned int destBuffer=0x00000000;
449  std::vector<unsigned short>::const_iterator uIt(src.begin());
450  dest.clear();
451 
452  for(i=0;i<15;i++)
453  {
454  dest.push_back(FloatRep(0,0,0));
455  }
456 
457  i=0;
458  for(;uIt!=src.end();++uIt)
459  {
460  int i0=m_tripVec[i].m_d;
461  int i1=m_tripVec[i].m_nd1;
462  int i2=m_tripVec[i].m_nd2;
463  buf=(*uIt);
464  char s=((buf & 0x0800)==0)?0:1;
465  dest[i1].sign(s);
466  s=((buf & 0x0020)==0)?0:1;
467  dest[i2].sign(s);
468  unsigned short e = ((buf & 0xF000) >> 12);
469  dest[i0].exponent(e+m_biases[i]);
470  e = ((buf & 0x07C0) >> 6);
471  dest[i1].exponent(e+m_biases[5]);
472  e = (buf & 0x001F);
473  dest[i2].exponent(e+m_biases[5]);
474  i++;
475  if(i==5) break;
476  }
477  if(i<5) return false;
478  nRestored=0;
479  nBitsToStore=0;
480  nFreeBits=0;
481  int nMantLength=23-m_bitStrip;
482  while(nRestored<=15)
483  {
484  if(nFreeBits==0)
485  {
486  if(nRestored!=0)
487  {
488  //printf("Dest 0x%X ",destBuffer);
489  destBuffer = destBuffer << m_bitStrip;
490  //printf("<< 0x%X = %lu\n",destBuffer,destBuffer);
491  dest[nRestored-1].mantissa(destBuffer);
492  }
493  ++nRestored;
494  destBuffer=0x00000000;
495  nFreeBits=nMantLength;
496  }
497  if(nBitsToStore==0)
498  {
499  if(uIt==src.end()) break;
500  ++uIt;
501  if(uIt==src.end())
502  {
503  //printf("Breaking ... nR=%d\n",nRestored);
504  break;
505  }
506  nBitsToStore=16;buf=(*uIt);//printf("Source 0x%X\n",buf);
507  }
508  int Np=(nFreeBits>nBitsToStore) ? nBitsToStore : nFreeBits;
509  unsigned int tmp = buf;
510  tmp = tmp >> (16-Np);
511  buf = (unsigned int)((buf << Np) & 0x0000FFFF);
512  //printf("F=%d TS=%d Np=%d 0x%X 0x%X\n",nFreeBits,nBitsToStore,Np,tmp,buf);
513  nBitsToStore-=Np;
514  //printf("Copy 0x%X ",destBuffer);
515  destBuffer = destBuffer << Np;
516  //printf("<< 0x%X ",destBuffer);
517  destBuffer = destBuffer | tmp;
518  //printf(" | 0x%X\n",destBuffer);
519  nFreeBits-=Np;
520  //printf("nR=%d\n",nRestored);
521  }
522  return true;
523 }
524 
xAOD::iterator
JetConstituentVector::iterator iterator
Definition: JetConstituentVector.cxx:68
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
ErrorMatrixCompressor::restoreFR
bool restoreFR(const std::vector< unsigned short > &, std::vector< FloatRep > &)
Definition: ErrorMatrixCompressor.cxx:444
ReadCellNoiseFromCoolCompare.s1
s1
Definition: ReadCellNoiseFromCoolCompare.py:378
ErrorMatrixCompressor::compress
bool compress(const std::vector< double > &, std::vector< unsigned int > &)
Definition: ErrorMatrixCompressor.cxx:173
python.SystemOfUnits.s
int s
Definition: SystemOfUnits.py:131
ErrorMatrixCompressor::m_bitStrip
int m_bitStrip
Definition: ErrorMatrixCompressor.h:157
xAOD::short
short
Definition: Vertex_v1.cxx:165
DecoderFloat_IEEE754::setExponent
void setExponent(short int)
Definition: ErrorMatrixCompressor.cxx:61
DecoderFloat_IEEE754::getExponent
short int getExponent()
Definition: ErrorMatrixCompressor.cxx:48
python.SystemOfUnits.m
int m
Definition: SystemOfUnits.py:91
CaloCellPos2Ntuple.int
int
Definition: CaloCellPos2Ntuple.py:24
WriteCellNoiseToCool.src
src
Definition: WriteCellNoiseToCool.py:513
ErrorMatrixCompressor::m_biases
short int m_biases[6]
Definition: ErrorMatrixCompressor.h:153
DecoderFloat_IEEE754::setSign
void setSign(int)
Definition: ErrorMatrixCompressor.cxx:69
DMTest::C
C_v1 C
Definition: C.h:26
skel.it
it
Definition: skel.GENtoEVGEN.py:396
ErrorMatrixCompressor.h
ErrorMatrixCompressor::CholeskyDecomposition
bool CholeskyDecomposition(double a[5][5], double L[5][5])
Definition: ErrorMatrixCompressor.cxx:141
ErrorMatrixCompressor::m_srcMasks
unsigned int m_srcMasks[16]
Definition: ErrorMatrixCompressor.h:161
UploadAMITag.l
list l
Definition: UploadAMITag.larcaf.py:158
perfmonmt-printer.dest
dest
Definition: perfmonmt-printer.py:189
ErrorMatrixCompressor::setUpperLimits
void setUpperLimits(const int l[2])
Definition: ErrorMatrixCompressor.cxx:132
JetTiledMap::S
@ S
Definition: TiledEtaPhiMap.h:44
ErrorMatrixCompressor::m_decoder
DecoderFloat_IEEE754 m_decoder
Definition: ErrorMatrixCompressor.h:160
python.utils.AtlRunQueryLookup.mask
string mask
Definition: AtlRunQueryLookup.py:460
DecoderFloat_IEEE754::FloatLongIntUnion::l
unsigned int l
Definition: ErrorMatrixCompressor.h:26
DecoderFloat_IEEE754::getMantissa
unsigned int getMantissa()
Definition: ErrorMatrixCompressor.cxx:54
ErrorMatrixCompressor::setBiases
void setBiases(const int b[6])
Definition: ErrorMatrixCompressor.cxx:124
DecoderFloat_IEEE754::m_data
FloatLongIntUnion m_data
Definition: ErrorMatrixCompressor.h:60
ErrorMatrixCompressor::Triplet
Definition: ErrorMatrixCompressor.h:126
ErrorMatrixCompressor::m_tripVec
std::vector< Triplet > m_tripVec
Definition: ErrorMatrixCompressor.h:162
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:210
createCoolChannelIdFile.buffer
buffer
Definition: createCoolChannelIdFile.py:12
convertTimingResiduals.sum
sum
Definition: convertTimingResiduals.py:55
perfmonmt-refit.slice
slice
Definition: perfmonmt-refit.py:52
lumiFormat.i
int i
Definition: lumiFormat.py:85
ErrorMatrixCompressor::restore
bool restore(const std::vector< unsigned int > &, std::vector< double > &)
Definition: ErrorMatrixCompressor.cxx:249
ErrorMatrixCompressor::m_scales
double m_scales[5]
Definition: ErrorMatrixCompressor.h:155
sign
int sign(int a)
Definition: TRT_StrawNeighbourSvc.h:107
DecoderFloat_IEEE754::print
void print()
Definition: ErrorMatrixCompressor.cxx:24
ErrorMatrixCompressor::ErrorMatrixCompressor
ErrorMatrixCompressor(int)
Definition: ErrorMatrixCompressor.cxx:82
DeMoUpdate.tmp
string tmp
Definition: DeMoUpdate.py:1167
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:77
FloatRep
Definition: ErrorMatrixCompressor.h:64
a
TList * a
Definition: liststreamerinfos.cxx:10
ReadCellNoiseFromCoolCompare.s2
s2
Definition: ReadCellNoiseFromCoolCompare.py:379
LArNewCalib_DelayDump_OFC_Cali.idx
idx
Definition: LArNewCalib_DelayDump_OFC_Cali.py:69
ErrorMatrixCompressor::m_limits
short int m_limits[2]
Definition: ErrorMatrixCompressor.h:154
ErrorMatrixCompressor::compressFR
bool compressFR(const std::vector< FloatRep > &, std::vector< unsigned short > &)
Definition: ErrorMatrixCompressor.cxx:308
DecoderFloat_IEEE754::setF
void setF(float f)
Definition: ErrorMatrixCompressor.h:35
updateCoolNtuple.limit
int limit
Definition: updateCoolNtuple.py:45
TSU::T
unsigned long long T
Definition: L1TopoDataTypes.h:35
fitman.k
k
Definition: fitman.py:528
Generated on Thu Nov 7 2024 21:14:50 for ATLAS Offline Software by   doxygen 1.8.18