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