StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StMultyKeyMap.cxx
1 #include <stdlib.h>
2 #include <stdio.h>
3 #include <math.h>
4 #include <string.h>
5 #include <assert.h>
6 static int gMyId=0;
7 
8 #include "StMultyKeyMap.h"
9 #ifndef MAX
10 #define MAX(a,b) ((a) > (b) ? (a) : (b))
11 #define MIN(a,b) ((a) < (b) ? (a) : (b))
12 #endif
13 //______________________________________________________________________________
14 StMultyKeyMap::StMultyKeyMap(int nkeys,int nbins)
15 {
16  mNKey=nkeys;
17  mNBin=nbins;
18  mJKey = 1946;
19  mLink.resize(mNBin, 0);
20  mDow.resize (mNKey, 1e11);
21  mStp.resize (mNKey,-1e11);
22  SetMap(this);
23 }
24 //______________________________________________________________________________
25 StMultyKeyMap::~StMultyKeyMap()
26 {
27 
28 }
29 //______________________________________________________________________________
30 void StMultyKeyMap::Clear(const char *)
31 {
32  mSize = 0;
33  mJKey = 1946;
34 
35  StMultyKeyDivd::Clear();
36  mArr.clear();
37 }
38 //______________________________________________________________________________
39 void StMultyKeyMap::Add(const void *obj,const double *keys)
40 {
41  float buf[200];
42  for (int i=0;i<mNKey;i++) {buf[i]=(float)keys[i];}
43  Add(obj,buf);
44 }
45 //______________________________________________________________________________
46 void StMultyKeyMap::Add(const void *obj,const float *keys)
47 {
48 assert(obj);
49  StMultyKeyNode *node = MakePair();
50  for (int ik=0;ik<mNKey;ik++) {
51  if (mDow[ik]>keys[ik]) mDow[ik]=keys[ik];
52  if (mStp[ik]<keys[ik]) mStp[ik]=keys[ik];
53  }
54  node->Set(obj,keys);
55  mArr.push_back(node);
56  mSize = mArr.size();
57 }
58 //______________________________________________________________________________
59 int StMultyKeyMap::GetJKey()
60 {
61  mJKey+=1000000007; return mJKey%mNKey;
62 }
63 //______________________________________________________________________________
64 StMultyKeyPair *StMultyKeyMap::MakePair()
65 {
66  StMultyKeyPair *pair = 0;
67  if (!mPairs.empty()) {// reuse
68  pair = mPairs.front(); mPairs.pop_front();}
69  else {
70  pair = new StMultyKeyPair();
71  pair->mKeys.resize(mNKey,0);
72  pair->SetMap(this);
73  }
74  return pair;
75 }
76 //______________________________________________________________________________
77 StMultyKeyDivd *StMultyKeyMap::MakeNode()
78 {
79  StMultyKeyDivd *node = 0;
80  if (!mNodes.empty()) {// reuse
81  node = mNodes.front();mNodes.pop_front();}
82  else {
83  node = new StMultyKeyDivd();
84  node->SetMap(this);
85  }
86  node->mLink.resize(mNBin,0);
87 
88  return node;
89 }
90 //______________________________________________________________________________
91 int StMultyKeyMap::MakeTree(int keepArray)
92 {
93  int nNodes = mArr.size();
94  if (!nNodes) return 0;
95  for (int ik=0;ik<mNKey;ik++) {mStp[ik] = 1.001*(mStp[ik]-mDow[ik])/mNBin;}
96 
97  SetIKey(0);
98 
99  for (int i=0;i<nNodes;i++) {Add(mArr[i]);}
100 
101  if (keepArray) return nNodes;
102 std::vector<StMultyKeyNode*> tmp(0);
103  mArr.swap(tmp); //destroy internal array completely;
104  return nNodes;
105 }
106 //______________________________________________________________________________
107 StMultyKeyDivd::StMultyKeyDivd()
108 {
109  mIKey = 0;
110 }
111 //______________________________________________________________________________
112 StMultyKeyDivd::~StMultyKeyDivd()
113 {
114 //delete [] mKeys;
115  delete mLink[0];
116  delete mLink[1];
117 }
118 //______________________________________________________________________________
119 //______________________________________________________________________________
120 int StMultyKeyNode::GetNumb() const
121 {
122  int nb = 0;
123  if (GetKeys()) return 1;
124  int nk = GetNKey();
125  for (int ik=0;ik<nk;ik++) {if (Link(ik)) nb += Link(ik)->GetNumb();}
126  return nb;
127 }
128 //______________________________________________________________________________
129 //______________________________________________________________________________
130 void StMultyKeyDivd::Clear(const char *)
131 {
132  int nk = GetNKey();
133  for (int ik=0;ik<nk;ik++) {if (mLink[ik]) mLink[ik]->Clear();mLink[ik]=0;}
134  mIKey = 0;
135  if (this == mMap) return;
136  mMap->mNodes.push_front(this);
137 }
138 //______________________________________________________________________________
139 StMultyKeyPair::StMultyKeyPair()
140 {
141  mId = ++gMyId; mObj = 0;
142 }
143 //______________________________________________________________________________
144 void StMultyKeyPair::Clear(const char *)
145 {
146  mKeys.clear(); mObj = 0;
147  mMap->mPairs.push_front(this);
148 }
149 //______________________________________________________________________________
150 void StMultyKeyPair::Set(const void *obj,const float *keys)
151 {
152  mObj = obj;;
153  mKeys.clear();
154  int nk = mMap->GetNKey();
155  for (int i=0;i<nk;i++) {mKeys.push_back(keys[i]);}
156 }
157 //______________________________________________________________________________
158 void StMultyKeyDivd::Add(StMultyKeyNode *pair)
159 {
160 static int nCall=0; nCall++;
161  const float *keyA = pair->GetKeys();
162  int nBin = mMap->GetNBin();
163  int ibin = (keyA[mIKey]-mDow[mIKey])/mStp[mIKey];
164 assert(ibin>=0 && ibin<nBin);
165  StMultyKeyNode *node = mLink[ibin];
166  const float *keyB = 0;
167  if (!node) { mLink[ibin] = pair; return;}
168  else if (!(keyB=node->GetKeys())) { node->Add(pair); return;}
169 
170 // Node is Pair
171  StMultyKeyDivd *bode = mMap->MakeNode();
172  mLink[ibin]=bode;
173  bode->mDow = mDow; bode->mStp = mStp;
174  bode->mDow[mIKey] += mStp[mIKey]*ibin;
175  bode->mStp[mIKey] /= nBin;
176  bode->mIKey=mMap->GetJKey();;
177  bode->Add(pair);
178  bode->Add(node);
179  return;
180 }
181 //______________________________________________________________________________
182 double StMultyKeyNode::Quality(int *numb) const
183 {
184 static int nCall=0; nCall++;
185  if (GetKeys()) { *numb = 1; return 0;}
186 
187  int ni=0,nt=0,ni2=0;
188  double qi=0,qa=0;
189  int nb = GetNBin();
190  for (int ib=0;ib<nb;ib++) {
191  if (!Link(ib)) continue;
192  qi = Link(ib)->Quality(&ni);
193  nt +=ni; ni2 += ni*ni; qa +=qi*ni;
194  }
195  qa = (qa + nt*nt-ni2)/nt;
196  if (numb) {*numb = nt;} else { qa/=nt;}
197  return qa;
198 }
199 //______________________________________________________________________________
200 int StMultyKeyNode::MaxDeep(int *deep) const
201 {
202  if (GetKeys()) { *deep = 1; return 0;}
203  int nk = GetNKey();int maxi=0;
204  for (int ik=0;ik<nk;ik++) {
205  int ni = (deep)? *deep+1:1;
206 
207  int mxi = ni;
208  if (Link(ik)) mxi = Link(ik)->MaxDeep(&ni);
209  if (maxi<mxi) maxi=mxi;
210  }
211 
212  return maxi;
213 }
214 //______________________________________________________________________________
215 int StMultyKeyNode::ls(const char *file) const
216 {
217 
218  FILE *out=stdout;
219  if (!file) file="";
220  if (file && file[0]) out=fopen(file,"w");
221 
222  StMultyKeyMapIter iter(this);
223 
224  int n=0;
225  for (const StMultyKeyNode *node=0;(node = *iter);++iter) {
226  n++;
227  if (!out) continue;
228  int lev = iter.Level();
229  const float *keys = node->GetKeys();
230  fprintf(out,"%4d - ",n);
231  fprintf(out,"Lev(%d) \t(%10p keys=)",lev,(void*)node);
232  int nk = node->GetNKey();
233  for (int k=0;k<nk;k++) {fprintf(out,"%g ",keys[k]);}
234  fprintf(out,"\n");
235  }
236  return n;
237 }
238 
239 
240 //______________________________________________________________________________
241 //______________________________________________________________________________
242 StMultyKeyMapIter::StMultyKeyMapIter(const StMultyKeyNode *node,const float *kMin,const float *kMax)
243  :mMinMax(0),mStk(32)
244 {
245  if (!node) return;
246  Set(node,kMin,kMax);
247 }
248 //______________________________________________________________________________
249 void StMultyKeyMapIter::Reset()
250 {
251  memset(mTouched,0,sizeof(mTouched));
252  mStk.resize(32);
253  mLev = 1;
254  int ini = InitLev();
255  if (ini<0 && FullCheck()==0) return;
256  mLev = 2; ++(*this);
257 }
258 //______________________________________________________________________________
259 void StMultyKeyMapIter::Set(const StMultyKeyNode *node,const float *kMin,const float *kMax)
260 {
261  memset(mTouched,0,sizeof(mTouched));
262  mTop = node;
263  mStk.resize(100);
264  mNK = node->GetNKey();
265  mNB = node->GetNBin();
266  mKMin=0;mKMax=0;
267  if (kMin) {
268  mMinMax.resize(2*mNK);
269  mKMin = &mMinMax[0];
270  mKMax = mKMin+mNK;
271  int sk = mNK*sizeof(mKMin[0]);
272  memcpy(mKMin,kMin,sk);
273  memcpy(mKMax,kMax,sk);
274  }
275  mLev = 0; mStk[0].node=0;
276  mLev = 1; mStk[1].node=node;
277  int ini = InitLev();
278  mLev = 2;
279  if (ini<0 && !FullCheck()) return;
280  mLev = 2; ++(*this);
281 }
282 //______________________________________________________________________________
283 int StMultyKeyMapIter::InitLev()
284 {
285  const StMultyKeyNode *node = mStk[mLev].node;
286  mStk[mLev].rite = 0;
287  mStk[mLev].ibin = 0;
288  if (node->GetKeys()) return -1;
289  int iKey = node->GetIKey();
290  mStk[mLev].rite = mNB-1;
291  mStk[mLev].ibin =-1;
292  if (mKMin) {
293  float dow = node->GetDow()[iKey];
294  float stp = node->GetStp()[iKey];
295  int left = (mKMin[iKey]-dow)/stp; if (left< 0) {left= 0;}; if(left>=mNB) return 1;
296  int rite = (mKMax[iKey]-dow)/stp; if (rite>=mNB) {rite=mNB-1;}; if(rite< 0) return 2;
297  mStk[mLev].rite = rite;
298  mStk[mLev].ibin = left-1;
299  }
300  return 0;
301 }
302 //______________________________________________________________________________
303 void StMultyKeyMapIter::Update(const float *kMin,const float *kMax)
304 {
305  int sk = mNK*sizeof(mKMin[0]);
306  if (kMin) memcpy(mKMin,kMin,sk);
307  if (kMax) memcpy(mKMax,kMax,sk);
308 }
309 //______________________________________________________________________________
310 StMultyKeyMapIter::~StMultyKeyMapIter()
311 {
312 }
313 //______________________________________________________________________________
314 StMultyKeyMapIter &StMultyKeyMapIter::operator++()
315 {
316  --mLev;
317  while (mLev) {
318  myStk_t &stk =mStk[mLev];
319  const StMultyKeyNode* node = stk.node;
320  const StMultyKeyNode* binNode = 0;
321 // Find non zero bin
322  while (++stk.ibin<=stk.rite) {
323  binNode = node->Link(stk.ibin); if (binNode) break;
324  }
325  if (!binNode) { mLev--; continue;} //did not find
326  mStk[++mLev].node = binNode;
327  int ini = InitLev();
328  if (ini==0) { continue;} //normal case
329  if (ini>0) { mLev--; continue;}
330  if (FullCheck()) { mLev--; continue;}
331  break; // found the good one
332  }
333  return *this;
334 }
335 //______________________________________________________________________________
336 int StMultyKeyMapIter::FullCheck()
337 {
338 // check node 0=good,1=bad
339 
340  const StMultyKeyNode *node = mStk[mLev].node;
341  if (!mKMin) return 0;
342  const float *fk = node->GetKeys();
343  mTouched[2]++;
344  for (int k=0;k<mNK;k++) {
345  if (mKMin[k]>fk[k] || fk[k] >= mKMax[k]) return 1;
346  }
347  return 0;
348 }
349 //______________________________________________________________________________
350 //______________________________________________________________________________
351 //______________________________________________________________________________
352 #include "TRandom.h"
353 #include "TStopwatch.h"
354 #include <map>
355 
356 void StMultyKeyMap::Test()
357 {
358 printf("StMultyKeyMap::Test() started\n");
359  float rnd;
360  int nEVTS=50000;
361  StMultyKeyMap map(1);
362  for (int i=0;i<nEVTS;i++) {
363  rnd = 1 - (i+1.)/nEVTS;
364  map.Add((void*)(-1),&rnd);
365  }
366  map.MakeTree();
367  map.ls();
368  double qa = map.Quality();
369  printf(" Quality of tree = %g\n\n",qa);
370 
371  StMultyKeyMapIter iter(map.GetTop());
372  int n = 0;
373  float pre=0;
374 
375  printf("\n%d evts No bounds\n",nEVTS);
376  for (StMultyKeyNode *node=0;(node = *iter);++iter)
377  {
378  const float * keys = node->GetKeys();
379  if(!keys) continue;
380  n++;
381  rnd = keys[0];
382 // printf("%4d - %g %p\n",n,rnd,node);
383  assert(pre<=rnd);
384  pre = rnd;
385  }
386 assert(n==nEVTS);
387 printf("\nNo bounds OK, touched %d %d %d\n"
388  ,iter.Touched()[0],iter.Touched()[1],iter.Touched()[2]);
389 
390 
391  float kMin=0.5,kMax=0.6;
392  iter.Set(map.GetTop(),&kMin,&kMax);
393  pre=0;n=0;
394  int nEst = int((nEVTS)*(kMax-kMin)+0.5);
395 printf("\n%d ~evts bounds=%g %g\n",nEst,kMin,kMax);
396  for (StMultyKeyNode *node=0;(node = *iter);++iter)
397  {
398  rnd = node->GetKeys()[0];
399  n++; rnd = node->GetKeys()[0];
401  assert(pre<=rnd);
402  assert((kMin<=rnd) && (rnd < kMax));
403  pre = rnd;
404  }
405 printf("\nGot %d. Bounds OK, Touched %d %d %d\n",n
406  ,iter.Touched()[0],iter.Touched()[1],iter.Touched()[2]);
407 
408 }
409 //______________________________________________________________________________
410 void StMultyKeyMap::Test2()
411 {
412 printf("StMultyKeyMap::Test2() started\n");
413  int nBin = 50;
414  StMultyKeyMap map(4,nBin);
415  float key[4];
416  int nEvts = 50000;
417 
418  TStopwatch TW;
419  TW.Start();
420 // Fill the map
421  for (int iEv=0;iEv<nEvts ;iEv++) {
422  for (int ik=0;ik<4;ik++) { key[ik]= gRandom->Rndm();}
423  map.Add((void*)1,key);
424  }
425  map.MakeTree(1946);
426  TW.Stop();
427  printf ( "MakeTree Cpu = %g\n",TW.CpuTime());
428  TW.Print();
429 // map.ls();
430 // Check map quality
431  double qa = map.Quality();
432  int maxDeep = map.MaxDeep();
433  int etaDeep = log(nEvts)/log(2.)+0.5;
434  printf(" Quality of tree = %g maxDeep=%d etaDeep %d",qa,maxDeep,etaDeep);
435 
436 
437  float dow[4]={0, 0.1,0.2,0.3};
438  float upp[4]={0.2,0.3,0.4,0.5};
439  double ev = nEvts;for (int i=0;i<4;i++){ev*=(upp[i]-dow[i]);};
440 printf("\n%d ~evts \n",int(ev+0.5));
441  int nk = map.GetNKey();
442  int nSel = 0,nBad=0;
443  StMultyKeyMapIter iter(map.GetTop(),dow,upp);
444 
445  TW.Start(1);
446  int nkl = 10000;
447 for (int jkl=0;jkl<nkl;jkl++) {
448  iter.Reset();
449  nSel = 0;nBad=0;
450  for (StMultyKeyNode *node=0;(node = *iter);++iter)
451  {
452  nSel++;
453 // int good = 0;
454 // const float *key = node->GetKeys();
455 // for (int j=0;j<nk;j++) {if (key[j]>=dow[j] && key[j]<upp[j]) good++;}
456 // nBad += (good!=nk);
457 // printf("%4d - %g %g %g %g \n",nSel,key[0],key[1],key[2],key[3]);
458  }
459 }//end jkl
460  TW.Stop();
461  printf ( "Search Cpu = %g\n",TW.CpuTime()/nkl);
462  int nb = map.Size();
463  StMultyKeyNode **nodes = map.GetArr();
464  int nMust=0;
465  for (int i=0;i<nb;i++)
466  {
467  StMultyKeyNode *node = nodes[i];
468  const float *key = node->GetKeys();
469  int good = 1;
470  for (int k=0;k<nk;k++) {
471  if (dow[k]<=key[k] && key[k] < upp[k]) continue;
472  good=0;break;
473  }
474  if (!good) continue;
475 
476  nMust++;
477  break;
478  }
479  map.Clear();
480 printf("\nSelected %d bad %d and must be %d\n",nSel,nBad,nMust);
481 printf("Touched %d %d %d\n"
482  ,iter.Touched()[0],iter.Touched()[1],iter.Touched()[2]);
483 
484 }
485 
486 
487 
Generated by   doxygen 1.8.5