StRoot: StarRoot/StMultiKeyMap.cxx Source File

00001 #include <stdlib.h>
00002 #include <stdio.h>
00003 #include <math.h>
00004 #include <string.h>
00005 #include <assert.h>
00006 
00007 #include <algorithm>
00008 #include <numeric>
00009 static int gMyId=0,gMyInst=0;
00010 
00011 #include "StMultiKeyMap.h"
00012 #ifndef MAX
00013 #define MAX(a,b) ((a) > (b) ? (a) : (b))
00014 #define MIN(a,b) ((a) < (b) ? (a) : (b))
00015 #endif
00016 static void random_shuffle(std::vector<StMultiKeyNode*> &arr); // shuffle elements 
00017 //______________________________________________________________________________
00018 StMultiKeyMap::StMultiKeyMap(int nkeys) 
00019 {
00020   mNKey=nkeys;
00021   mTop = 0;
00022 }
00023 //______________________________________________________________________________
00024 StMultiKeyMap::~StMultiKeyMap() 
00025 {
00026   delete mTop;mTop=0;
00027 }
00028 //______________________________________________________________________________
00029 void StMultiKeyMap::Clear(const char *) 
00030 {
00031   delete mTop;mTop=0;
00032   mArr.clear();
00033 }
00034 //______________________________________________________________________________
00035 void StMultiKeyMap::Add(const void *obj,const double *keys)
00036 {
00037    if (!obj) obj = (void*)(-1);
00038    float buf[200];
00039    for (int i=0;i<mNKey;i++) {buf[i]=(float)keys[i];}
00040    Add(obj,buf);
00041 }
00042 //______________________________________________________________________________
00043 void StMultiKeyMap::Add(const void *obj,const float *keys)
00044 {
00045   if (!obj) obj = (void*)(-1);
00046   assert(!mTop);
00047   StMultiKeyNode *node = new StMultiKeyNode(mNKey);
00048   node->Set(obj,keys);
00049   mArr.push_back(node);
00050 }
00051 //______________________________________________________________________________
00052 double StMultiKeyMap::StMultiKeyMap::Quality() 
00053 {
00054   assert(mTop);
00055   return mTop->Quality();
00056 }
00057 //______________________________________________________________________________
00058 int StMultiKeyMap::MakeTree()
00059 {
00060    assert(!mTop);
00061    int nNodes = mArr.size();
00062    if (!nNodes) return 0;
00063 //   std::random_shuffle( mArr.begin(),mArr.end() ); // shuffle elements 
00064    random_shuffle(mArr); 
00065    mTop =  mArr[0];
00066    for (int i=1;i<nNodes;i++) {mTop->Add(mArr[i]);}
00067 
00068 std::vector<StMultiKeyNode*> tmp(0);
00069    assert(nNodes == mTop->Size());
00070    mArr.swap(tmp);      //destroy internal array completely;
00071    return nNodes;
00072 }
00073 //______________________________________________________________________________
00074 int StMultiKeyMap::ls(const char *file) const
00075 {
00076   return mTop->ls(file);
00077 }
00078 //______________________________________________________________________________
00079 int StMultiKeyMap::Size() const
00080 {
00081   if (mTop) return mTop->Size();
00082   else      return mArr.size();
00083 }
00084 //______________________________________________________________________________
00085 void StMultiKeyMap::Test()
00086 {
00087 printf("StMultiKeyMap::Test() started\n");
00088    float rnd; 
00089    int nEVTS=1000;
00090    StMultiKeyMap map(1);
00091    for (int i=0;i<nEVTS;i++) {
00092      rnd = (i+1.)/nEVTS; 
00093      map.Add(0,&rnd);
00094    }
00095    map.MakeTree();
00096    double qa = map.Quality();
00097    printf(" Quality of tree = %g\n\n",qa);
00098    map.ls();
00099 
00100    StMultiKeyMapIter iter(map.GetTop());
00101    int n = 0;
00102    float pre=0;
00103 
00104 printf("\n%d evts No bounds\n",nEVTS);
00105    for (StMultiKeyNode *node=0;(node = *iter);++iter)
00106    {
00107      n++; rnd = node->GetKeys()[0];
00109      assert(pre<=rnd);
00110      pre = rnd;
00111    }
00112 assert(n==nEVTS);
00113 printf("\nNo bounds OK, touched %d %d\n",iter.Touched()[0],iter.Touched()[1]);
00114 
00115 
00116    float kMin=0.5,kMax=0.6;
00117    iter.Set(map.GetTop(),&kMin,&kMax);
00118    pre=0;n=0;
00119    int nEst = int((nEVTS)*(kMax-kMin)+0.5);
00120 printf("\n%d ~evts bounds=%g %g\n",nEst,kMin,kMax);
00121    for (StMultiKeyNode *node=0;(node = *iter);++iter)
00122    {
00123      n++; rnd = node->GetKeys()[0];
00125      assert(pre<=rnd);
00126      assert((kMin<=rnd) && (rnd < kMax));
00127      pre = rnd;
00128    }
00129 printf("\nGot %d. Bounds OK, Touched %d %d\n",n,iter.Touched()[0],iter.Touched()[1]);
00130 
00131 }
00132 //______________________________________________________________________________
00133 void StMultiKeyMap::Test2()
00134 {
00135 printf("StMultiKeyMap::Test2() started\n");
00136    StMultiKeyMap map(4);
00137    float key[4];
00138    for (int ix=0;ix<100;ix+=10) {
00139      key[0]=ix;key[1]=ix+10;
00140      for (int iy=0;iy<200;iy+=10) {
00141      key[2]=iy;key[3]=iy+10;
00142      map.Add(0,key);
00143    } }
00144    map.MakeTree();
00145    double qa = map.Quality();
00146    printf(" Quality of tree = %g\n\n",qa);
00147    map.ls();
00148 
00149 
00150    float sel[4]={15,25,105,115};
00151    StMultiKeyMapIter iter(map.GetTop(),sel,0);
00152 printf("\n4 ~evts \n");
00153    int n = 0;
00154    for (StMultiKeyNode *node=0;(node = *iter);++iter)
00155    {
00156      n++; 
00157      const float *key = node->GetKeys();
00158      printf("%4d - %g %g %g %g \n",n,key[0],key[1],key[2],key[3]);
00159    }
00160 printf("\nGot %d. Bounds OK , Touched %d %d\n",n,iter.Touched()[0],iter.Touched()[1]);
00161 
00162 }
00163 
00164 //______________________________________________________________________________
00165 //______________________________________________________________________________
00166 StMultiKeyNode::StMultiKeyNode(int nKeys)
00167 {
00168   Init();
00169   mNKey=nKeys; 
00170 } 
00171 //______________________________________________________________________________
00172 StMultiKeyNode::StMultiKeyNode(const StMultiKeyNode &fr)
00173 {
00174   Init();
00175   mNKey = fr.mNKey;
00176   if (fr.mObj) Set(fr.mObj,fr.mKeys);
00177 } 
00178 //______________________________________________________________________________
00179 void StMultiKeyNode::Init()
00180 {
00181   memset(&mNKey,0,(char*)&mKeys-&mNKey+sizeof(mKeys));
00182   mId = ++gMyId; gMyInst++;
00183 }
00184 //______________________________________________________________________________
00185 void StMultiKeyNode::Clear()
00186 {
00187   memset(&mIKey,0,(char*)&mObj - &mIKey);
00188 }
00189 //______________________________________________________________________________
00190 StMultiKeyNode::~StMultiKeyNode()
00191 {
00192   if (mObj) delete [] mKeys;
00193   delete mLink[0];
00194   delete mLink[1];
00195   gMyInst--;
00196 } 
00197 //______________________________________________________________________________
00198 int StMultiKeyNode::GetNInst() 
00199 { return gMyInst;}
00200 //______________________________________________________________________________
00201 void StMultiKeyNode::Set(const void *obj,const float *keys)
00202 {
00203    Clear();
00204    if (!obj) obj = (void*)(-1);
00205    mObj = obj; mIKey=0;
00206    if (!mKeys) mKeys = new float[mNKey]; 
00207    memcpy(mKeys,keys,sizeof(mKeys[0])*mNKey);
00208 }
00209 //______________________________________________________________________________
00210 void StMultiKeyNode::Set(const void *obj,const double *keys)
00211 {
00212    float buf[200];
00213    for (int i=0;i<mNKey;i++) {buf[i]=(float)keys[i];}
00214    Set(obj,buf);
00215 }
00216 //______________________________________________________________________________
00217 void StMultiKeyNode::Add(const void *obj,const float *keys)
00218 {
00219   StMultiKeyNode *node = new StMultiKeyNode(mNKey);
00220   node->Set(obj,keys);
00221   Add(node);
00222 }
00223 //______________________________________________________________________________
00224 void StMultiKeyNode::Add(StMultiKeyNode *node)
00225 {
00226 static int nCall=0; nCall++;
00227   assert(this != node);
00228   node->mIKey = (mIKey+1000003)%mNKey;
00229   int way = (node->mKeys[int(mIKey)] <= GetKey())? 0:1;
00230   if (mLink[way])            { mLink[way]->Add(node);}
00231   else                       { mLink[way] = node    ;}
00232   mNumb[way]++;
00233   return;
00234 }  
00235 //______________________________________________________________________________
00236 double StMultiKeyNode::Quality() 
00237 {
00238    double qa = 0;
00239    int nTot = GetNumb(0)+GetNumb(1)+1;
00240    StMultiKeyMapIter iter(this);
00241    
00242    int n=0; 
00243    for (StMultiKeyNode *node=0;(node = *iter);++iter) {
00244      n++;
00245      int nL = node->GetNumb(0);
00246      int nR = node->GetNumb(1);
00247      if (!nL || !nR) continue;
00248      qa += (2.*nL*nR+nL+nR)/nTot;
00249    }
00250 printf("Quality() nodes %d\n",n);
00251    return qa/nTot;
00252 }
00253 //______________________________________________________________________________
00254 int StMultiKeyNode::ls(const char *file) const
00255 {
00256    FILE *out=stdout;
00257    if (!file) file="";
00258    if (file[0]=='-') out=0;
00259    else if (file[0]) out=fopen(file,"r");
00260 
00261    StMultiKeyMapIter iter((StMultiKeyNode*)this);
00262    
00263    int n=0; 
00264    for (const StMultiKeyNode *node=0;(node = *iter);++iter) {
00265      n++;
00266      if (!out) continue;
00267      int nL = node->GetNumb(0);
00268      int nR = node->GetNumb(1);
00269      int lev = iter.Level();
00270      if (node==this) {fprintf(out,"%4d * ",n);}
00271      else            {fprintf(out,"%4d - ",n);}
00272      fprintf(out,"Lev(%d) \t(%10p) \tL(%10p(%d)) \tR(%10p(%d)) \tDiv(%g)"
00273            ,lev,(void*)node
00274            ,(void*)node->mLink[0],nL
00275            ,(void*)node->mLink[1],nR
00276            ,node->GetKey());
00277      if (node->mObj) {
00278        for (int j=0;j<mNKey;j++) {fprintf(out," %g",node->mKeys[j]);}}
00279      fprintf(out,"\n");
00280    } 
00281    return n;
00282 }
00283 
00284 
00285 //______________________________________________________________________________
00286 //______________________________________________________________________________
00287 StMultiKeyMapIter::StMultiKeyMapIter(const StMultiKeyNode *node,const float *kMin,const float *kMax)
00288   :mMinMax(0),mStk(32)
00289 {
00290   if (!node) return;
00291   Set(node,kMin,kMax);
00292 }
00293 //______________________________________________________________________________
00294 void StMultiKeyMapIter::Set(const StMultiKeyNode *node,const float *kMin,const float *kMax)
00295 {
00296   mTouched[0]=0;mTouched[1]=0;
00297   mStk.resize(32);
00298   mBoundsOn=(kMin && !kMax);    //MapTree is a tree of boundaries. 
00299   mNK = node->GetNKey();
00300   mKMin=0;mKMax=0;
00301   if (kMin) {
00302     mMinMax.resize(2*mNK);
00303     mKMin = &mMinMax[0];
00304     mKMax  = (mBoundsOn) ? 0: mKMin+mNK;
00305     int sk = mNK*sizeof(mKMin[0]);
00306     memcpy(mKMin,kMin,sk);
00307     if (mKMax) memcpy(mKMax,kMax,sk);
00308   }
00309   mLev = 0; mStk[0]=0;
00310   Down(node);
00311   SelfCheck();
00312 }
00313 //______________________________________________________________________________
00314 void StMultiKeyMapIter::Update(const float *kMin,const float *kMax)
00315 {
00316   int sk = mNK*sizeof(mKMin[0]);
00317   if (kMin) memcpy(mKMin,kMin,sk);
00318   if (kMax) memcpy(mKMax,kMax,sk);
00319 }
00320 //______________________________________________________________________________
00321 StMultiKeyMapIter::~StMultiKeyMapIter()
00322 {
00323 }
00324 //______________________________________________________________________________
00325 StMultiKeyMapIter &StMultiKeyMapIter::operator++()
00326 {
00327   if (!mLev)    return *this;
00328   const StMultiKeyNode* node = mStk[mLev];
00329   const StMultiKeyNode* rLink  = node->RLink();
00330   mLev--;
00331   if (!rLink)                   goto RETN;
00332   if (mKMin && FilterRite(node))goto RETN;
00333   Down(rLink);
00334   return *this;
00335 RETN:
00336   SelfCheck();
00337   return *this;
00338 }  
00339 //______________________________________________________________________________
00340 void StMultiKeyMapIter::Down(const StMultiKeyNode *node) 
00341 {
00342   while(node) {
00343     if ((int)mStk.size() <=mLev) mStk.resize(mLev*2);
00344     mStk[++mLev] = node;
00345     if (!node->LLink()) break;
00346     if (mKMin && FilterLeft(node)) break;
00347     node = node->LLink();
00348   }
00349   SelfCheck();
00350 }
00351 //______________________________________________________________________________
00352 void StMultiKeyMapIter::SelfCheck() 
00353 {
00354   const StMultiKeyNode *node = mStk[mLev];
00355   if (!node ) return;
00356   if (!node->GetObj()) {++(*this); return;}
00357   if (!mKMin) return;
00358   const float *fk = node->GetKeys();
00359   switch (mBoundsOn) {
00360     case 0:
00361       mTouched[1]++;
00362       for (int k=0;k<mNK;k++) {  
00363         if (mKMin[k]<=fk[k] && fk[k] < mKMax[k]) continue;
00364          ++(*this); return;
00365       }
00366     return;
00367 
00368     case 1:
00369       mTouched[1]++;
00370       for (int k=0;k<mNK;k++) {  
00371         float *lim = mKMin+(k&(-2));
00372         switch (k&1) {
00373           case 0: if (fk[k]>lim[1]) {++(*this); return;} break;
00374           case 1: if (fk[k]<lim[0]) {++(*this); return;} break;
00375         }//end switch
00376       }//end for k
00377     }//end boundOn switch 
00378   return;
00379 }  
00380 //______________________________________________________________________________
00381 int StMultiKeyMapIter::FilterLeft(const StMultiKeyNode *node) const
00382 {
00383   int   ik = node->GetIKey();
00384   float fk = node->GetKey();
00385   assert(ik>=0);
00386 
00387   switch (mBoundsOn) {
00388   
00389   case 0: //Normal case
00390     if ( mKMin[ik]>fk) return 1; return 0;
00391 
00392   case 1:
00393     float *lim = mKMin+(ik&(-2));
00394     switch(ik&1) {
00395     
00396       case 0: //key is a lower boundary
00397         return 0; 
00398         
00399       case 1://key is a upper boundary
00400         mTouched[0]++;
00401         if (lim[0]>fk) return 1; return 0;
00402     }//end ik switch
00403   }//end bounds switch
00404   return 0;
00405 }
00406 //______________________________________________________________________________
00407 int StMultiKeyMapIter::FilterRite(const StMultiKeyNode *node) const
00408 {
00409   int ik = node->GetIKey();
00410   assert(ik>=0);
00411   float fk = node->GetKey();
00412   switch (mBoundsOn) {
00413 
00414    case 0:   mTouched[0]++;
00415              if (mKMax[ik]<fk) return 1; return 0;
00416 
00417    case 1: float *lim = mKMin+(ik&(-2));
00418     switch(ik&1) {
00419     
00420       case 0: //key is a lower boundary
00421          mTouched[0]++; if (lim[1]<fk) return 1; return 0;
00422         
00423       case 1://key is a upper boundary
00424         return 0;
00425     }//end ik switch
00426   }//end bounds switch
00427   return 0;
00428 }
00429 //______________________________________________________________________________
00430 void random_shuffle(std::vector<StMultiKeyNode*> &arr)
00431 {
00432   int n = arr.size(); if (n<=3) return;
00433   unsigned int u=n/2,us=1000000007;
00434   int jr=n-1;
00435   while (0<jr) {
00436     int jj = (u+=us)%jr;
00437     StMultiKeyNode *v = arr[jj]; arr[jj]=arr[jr]; arr[jr]=v; jr--;
00438   }
00439 }
00440 //______________________________________________________________________________
00441 //______________________________________________________________________________