StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StMultiKeyMapM.cxx.C
1 #include <stdlib.h>
2 #include <stdio.h>
3 #include <math.h>
4 #include <string.h>
5 #include <assert.h>
6 
7 #include <algorithm>
8 #include <numeric>
9 static int gMyId=0,gMyInst=0;
10 
11 #include "StMultiKeyMapM.h"
12 #ifndef MAX
13 #define MAX(a,b) ((a) > (b) ? (a) : (b))
14 #define MIN(a,b) ((a) < (b) ? (a) : (b))
15 #endif
16 static void random_shuffle(std::vector<StMultiKeyNode*> &arr); // shuffle elements
17 //______________________________________________________________________________
18 StMultiKeyMapM::StMultiKeyMapM(int nkeys)
19 {
20  mNKey=nkeys;
21  mTop = 0;
22 }
23 //______________________________________________________________________________
24 StMultiKeyMapM::~StMultiKeyMapM()
25 {
26  delete mTop;mTop=0;
27 }
28 //______________________________________________________________________________
29 void StMultiKeyMapM::Clear(const char *)
30 {
31  delete mTop;mTop=0;
32  mArr.clear();
33 }
34 //______________________________________________________________________________
35 void StMultiKeyMapM::Add(const void *obj,const double *keys)
36 {
37  float buf[200];
38  for (int i=0;i<mNKey;i++) {buf[i]=(float)keys[i];}
39  Add(obj,buf);
40 }
41 //______________________________________________________________________________
42 void StMultiKeyMapM::Add(const void *obj,const float *keys)
43 {
44 assert(obj);
45  assert(!mTop);
46  StMultiKeyNode *node = new StMultiKeyNode(mNKey);
47  node->Set(obj,keys);
48  mArr.push_back(node);
49 }
50 //______________________________________________________________________________
51 double StMultiKeyMapM::StMultiKeyMapM::Quality()
52 {
53  assert(mTop);
54  return mTop->Quality();
55 }
56 //______________________________________________________________________________
57 int StMultiKeyMapM::MakeTree(int keepArray)
58 {
59  int nNodes = mArr.size();
60  if (!nNodes) return 0;
61 // std::random_shuffle( mArr.begin(),mArr.end() ); // shuffle elements
62  random_shuffle(mArr);
63  int jl = 0;
64  if (!mTop ) {mTop = mArr[0];jl=1;mTop->SetIKey(0);}
65  unsigned int myRnd = 1946;
66  for (int i=jl;i<nNodes;i++) {
67  myRnd+=1000000007; mArr[i]->SetIKey(myRnd%mNKey);
68  mTop->Add(mArr[i]);
69  }
70 
71  if (keepArray) return nNodes;
72 std::vector<StMultiKeyNode*> tmp(0);
73 // assert(nNodes == mTop->Size());
74  mArr.swap(tmp); //destroy internal array completely;
75  return nNodes;
76 }
77 //______________________________________________________________________________
78 int StMultiKeyMapM::ls(const char *file) const
79 {
80  return mTop->ls(file);
81 }
82 //______________________________________________________________________________
83 int StMultiKeyMapM::Size() const
84 {
85  if (mTop) return mTop->Size();
86  else return mArr.size();
87 }
88 //______________________________________________________________________________
89 //______________________________________________________________________________
90 StMultiKeyNode::StMultiKeyNode(int nKeys)
91 {
92  Init();
93  mNKey=nKeys;
94 }
95 //______________________________________________________________________________
96 StMultiKeyNode::StMultiKeyNode(const StMultiKeyNode &fr)
97 {
98  Init();
99  mNKey = fr.mNKey;
100  if (fr.mObj) Set(fr.mObj,fr.mKeys);
101 }
102 //______________________________________________________________________________
103 void StMultiKeyNode::Init()
104 {
105  memset(&mNKey,0,(char*)&mKeys-&mNKey+sizeof(mKeys));
106  mId = ++gMyId; gMyInst++; mIKey = -1;
107 }
108 //______________________________________________________________________________
109 void StMultiKeyNode::Clear()
110 {
111  memset(&mIKey,0,(char*)&mObj - &mIKey); mIKey = -1;
112 }
113 //______________________________________________________________________________
114 StMultiKeyNode::~StMultiKeyNode()
115 {
116  delete [] mKeys;
117  delete mLink[0];
118  delete mLink[1];
119  gMyInst--;
120 }
121 //______________________________________________________________________________
122 int StMultiKeyNode::GetNInst()
123 { return gMyInst;}
124 //______________________________________________________________________________
125 void StMultiKeyNode::Set(const void *obj,const float *keys)
126 {
127  Clear();
128  mObj = obj; mIKey=-1;
129  if (!mKeys) mKeys = new float[mNKey];
130  memcpy(mKeys,keys,sizeof(mKeys[0])*mNKey);
131 }
132 //______________________________________________________________________________
133 void StMultiKeyNode::Set(const void *obj,const double *keys)
134 {
135  float buf[200];
136  for (int i=0;i<mNKey;i++) {buf[i]=(float)keys[i];}
137  Set(obj,buf);
138 }
139 //______________________________________________________________________________
140 void StMultiKeyNode::Add(const void *obj,const float *keys)
141 {
142  StMultiKeyNode *node = new StMultiKeyNode(mNKey);
143  node->Set(obj,keys);
144  Add(node);
145 }
146 //______________________________________________________________________________
147 void StMultiKeyNode::Add(StMultiKeyNode *node)
148 {
149 static int nCall=0; nCall++;
150  assert(this != node);
151  int way = (node->mKeys[int(mIKey)] > GetKey());
152  if (mLink[way]) { mLink[way]->Add(node);}
153  else { mLink[way] = node ;}
154  mNumb[way]++;
155  return;
156 }
157 //______________________________________________________________________________
158 double StMultiKeyNode::Quality()
159 {
160  double qa = 0;
161  int nTot = GetNumb(0)+GetNumb(1)+1;
162  StMultiKeyMapMIter iter(this);
163 
164  int n=0;
165  for (StMultiKeyNode *node=0;(node = *iter);++iter) {
166  n++;
167  int nL = node->GetNumb(0);
168  int nR = node->GetNumb(1);
169  if (!nL || !nR) continue;
170  qa += (2.*nL*nR+nL+nR)/nTot;
171  }
172 printf("Quality() nodes %d\n",n);
173  return qa/nTot;
174 }
175 //______________________________________________________________________________
176 int StMultiKeyNode::ls(const char *file) const
177 {
178  FILE *out=stdout;
179  if (!file) file="";
180  if (file && file[0]) out=fopen(file,"w");
181 
182  StMultiKeyMapMIter iter((StMultiKeyNode*)this);
183 
184  int n=0;
185  for (const StMultiKeyNode *node=0;(node = *iter);++iter) {
186  n++;
187  if (!out) continue;
188  int nL = node->GetNumb(0);
189  int nR = node->GetNumb(1);
190  int lev = iter.Level();
191  if (node==this) {fprintf(out,"%4d * ",n);}
192  else {fprintf(out,"%4d - ",n);}
193  fprintf(out,"Lev(%d) \t(%10p) \tL(%10p(%d)) \tR(%10p(%d)) \tDiv(%g)"
194  ,lev,(void*)node
195  ,(void*)node->mLink[0],nL
196  ,(void*)node->mLink[1],nR
197  ,node->GetKey());
198  if (node->mObj) {
199  for (int j=0;j<mNKey;j++) {fprintf(out," %g",node->mKeys[j]);}}
200  fprintf(out,"\n");
201  }
202  return n;
203 }
204 
205 
206 //______________________________________________________________________________
207 //______________________________________________________________________________
208 StMultiKeyMapMIter::StMultiKeyMapMIter(const StMultiKeyNode *node,const float *kMin,const float *kMax)
209  :mMinMax(0),mStk(32)
210 {
211  if (!node) return;
212  Set(node,kMin,kMax);
213 }
214 //______________________________________________________________________________
215 void StMultiKeyMapMIter::Set(const StMultiKeyNode *node,const float *kMin,const float *kMax)
216 {
217  memset(mTouched,0,sizeof(mTouched));
218  mStk.resize(32);
219  mNK = node->GetNKey();
220  mKMin=0;mKMax=0;
221  if (kMin) {
222  mMinMax.resize(2*mNK);
223  mKMin = &mMinMax[0];
224  mKMax = mKMin+mNK;
225  int sk = mNK*sizeof(mKMin[0]);
226  memcpy(mKMin,kMin,sk);
227  memcpy(mKMax,kMax,sk);
228  }
229  mLev = 0; mStk[0]=0;
230 
231  Left(node);
232  if (FullCheck()) ++(*this);
233 }
234 //______________________________________________________________________________
235 void StMultiKeyMapMIter::Update(const float *kMin,const float *kMax)
236 {
237  int sk = mNK*sizeof(mKMin[0]);
238  if (kMin) memcpy(mKMin,kMin,sk);
239  if (kMax) memcpy(mKMax,kMax,sk);
240 }
241 //______________________________________________________________________________
242 StMultiKeyMapMIter::~StMultiKeyMapMIter()
243 {
244 }
245 //______________________________________________________________________________
246 StMultiKeyMapMIter &StMultiKeyMapMIter::operator++()
247 {
248  while (mLev) {
249  const StMultiKeyNode* node = mStk[mLev];
250  const StMultiKeyNode* rLink = node->RLink();
251  mLev--;
252  if (rLink && (!mKMin || !FilterRite(node))) Left(rLink);
253  if (!FullCheck()) break;
254  }
255  return *this;
256 }
257 //______________________________________________________________________________
258 void StMultiKeyMapMIter::Left(const StMultiKeyNode *node)
259 {
260  while(node) {
261  if ((int)mStk.size() <=mLev) mStk.resize(mLev*2);
262  mStk[++mLev] = node;
263  if (!node->LLink()) break;
264  if (mKMin && FilterLeft(node)) break;
265  node = node->LLink();
266  }
267  return;
268 }
269 //______________________________________________________________________________
270 int StMultiKeyMapMIter::FullCheck()
271 {
272  const StMultiKeyNode *node = mStk[mLev];
273 
274  if (!node ) return 0;
275  if (!mKMin) return 0;
276 // if (!node->GetObj()) { ++(*this); return;}
277  const float *fk = node->GetKeys();
278  mTouched[2]++;
279  for (int k=0;k<mNK;k++) {
280  if (mKMin[k]>fk[k] || fk[k] >= mKMax[k]) return 1;
281  }
282  return 0;
283 }
284 //______________________________________________________________________________
285 int StMultiKeyMapMIter::FilterLeft(const StMultiKeyNode *node) const
286 {
287  int ik = node->GetIKey();
288  float fk = node->GetKey();
289  mTouched[0]++;
290  return ( mKMin[ik]>fk);
291 }
292 //______________________________________________________________________________
293 int StMultiKeyMapMIter::FilterRite(const StMultiKeyNode *node) const
294 {
295  int ik = node->GetIKey();
296  float fk = node->GetKey();
297  mTouched[1]++;
298  return (mKMax[ik]<=fk);
299 }
300 //______________________________________________________________________________
301 void random_shuffle(std::vector<StMultiKeyNode*> &arr)
302 {
303 static const unsigned int us=1000000007;
304  int n = arr.size(); if (n<=3) return;
305  unsigned int u=n/2;
306  int jr=n-1;
307  while (0<jr) {
308  int jj = (u+=us)%jr;
309  StMultiKeyNode *v = arr[jj]; arr[jj]=arr[jr]; arr[jr]=v; jr--;
310  }
311 }
312 //______________________________________________________________________________
313 //______________________________________________________________________________
314 //______________________________________________________________________________
315 #include "TRandom.h"
316 void StMultiKeyMapM::Test()
317 {
318 printf("StMultiKeyMapM::Test() started\n");
319  float rnd;
320  int nEVTS=1000;
321  StMultiKeyMapM map(1);
322  for (int i=0;i<nEVTS;i++) {
323  rnd = 1 - (i+1.)/nEVTS;
324  map.Add((void*)(-1),&rnd);
325  }
326  map.MakeTree();
327  map.ls();
328  double qa = map.Quality();
329  printf(" Quality of tree = %g\n\n",qa);
330 
331  StMultiKeyMapMIter iter(map.GetTop());
332  int n = 0;
333  float pre=0;
334 
335 printf("\n%d evts No bounds\n",nEVTS);
336  for (StMultiKeyNode *node=0;(node = *iter);++iter)
337  {
338  n++; rnd = node->GetKeys()[0];
339 // printf("%4d - %g %p\n",n,rnd,node);
340  assert(pre<=rnd);
341  pre = rnd;
342  }
343 assert(n==nEVTS);
344 printf("\nNo bounds OK, touched %d %d %d\n"
345  ,iter.Touched()[0],iter.Touched()[1],iter.Touched()[2]);
346 
347 
348  float kMin=0.5,kMax=0.6;
349  iter.Set(map.GetTop(),&kMin,&kMax);
350  pre=0;n=0;
351  int nEst = int((nEVTS)*(kMax-kMin)+0.5);
352 printf("\n%d ~evts bounds=%g %g\n",nEst,kMin,kMax);
353  for (StMultiKeyNode *node=0;(node = *iter);++iter)
354  {
355  n++; rnd = node->GetKeys()[0];
357  assert(pre<=rnd);
358  assert((kMin<=rnd) && (rnd < kMax));
359  pre = rnd;
360  }
361 printf("\nGot %d. Bounds OK, Touched %d %d %d\n",n
362  ,iter.Touched()[0],iter.Touched()[1],iter.Touched()[2]);
363 
364 }
365 //______________________________________________________________________________
366 void StMultiKeyMapM::Test2()
367 {
368 printf("StMultiKeyMapM::Test2() started\n");
369  StMultiKeyMapM map(4);
370  float key[4];
371  int nEvts = 10000;
372  for (int iEv=0;iEv<nEvts ;iEv++) {
373  for (int ik=0;ik<4;ik++) { key[ik]= gRandom->Rndm();}
374  map.Add((void*)1,key);
375  }
376 
377  map.MakeTree(1946);
378 // map.ls();
379  double qa = map.Quality();
380  printf(" Quality of tree = %g\n\n",qa);
381 
382 
383  float dow[4]={0, 0.1,0.2,0.3};
384  float upp[4]={0.2,0.3,0.4,0.5};
385 
386  StMultiKeyMapMIter iter(map.GetTop(),dow,upp);
387  double ev = nEvts;for (int i=0;i<4;i++){ev*=(upp[i]-dow[i]);};
388 printf("\n%d ~evts \n",int(ev+0.5));
389  int nk = map.GetNKey();
390  int nSel = 0,nBad=0;
391  for (StMultiKeyNode *node=0;(node = *iter);++iter)
392  {
393  nSel++; int good = 0;
394  const float *key = node->GetKeys();
395  for (int j=0;j<nk;j++) {if (key[j]>=dow[j] && key[j]<upp[j]) good++;}
396  nBad += (good!=nk);
397 // printf("%4d - %g %g %g %g \n",nSel,key[0],key[1],key[2],key[3]);
398  }
399  int nb = map.Size();
400  StMultiKeyNode **nodes = map.GetArr();
401  int nMust=0;
402  for (int i=0;i<nb;i++)
403  {
404  StMultiKeyNode *node = nodes[i];
405  const float *key = node->GetKeys();
406  int good = 1;
407  for (int k=0;k<nk;k++) {
408  if (dow[k]<=key[k] && key[k] < upp[k]) continue;
409  good=0;break;
410  }
411  if (!good) continue;
412  nMust++;
413  }
414  map.Clear();
415 printf("\nSelected %d bad %d and must be %d\n",nSel,nBad,nMust);
416 printf("Touched %d %d %d\n"
417  ,iter.Touched()[0],iter.Touched()[1],iter.Touched()[2]);
418 
419 }
420