STAT/Macros/testkdtreeIF.C

   1 // #include <malloc.h>
   2 // #include "TMatrixD.h"
   3 // #include "TRandom.h"
   4 // #include "TGraph.h"
   5 // #include "TStopwatch.h"
   6 //
   7
   8
   9 Int_t Mem()
  10 {
  11   // size in 1000 bytes
  12   static struct mallinfo memdebug;
  13   memdebug = mallinfo();
  14   return memdebug.uordblks/1000;
  15 }
  16
  17
  18
  19 void testindexes(Int_t nloop, Int_t npoints);
  20 void  testkdtreeIF(Int_t npoints=1000, Int_t nloop=1000, Int_t mode = 1, Int_t bsize=9);
  21 void TestSizeIF(Int_t npoints=1000, Int_t nrows = 150, Int_t nsec=18, Int_t mode = 1, Int_t bsize=9);
  22
  23
  24
  25 void testindexes(Int_t nloop, Int_t npoints){
  26   //
  27   // test indexing
  28   //
  29   TKDTree<Int_t, Float_t> *kdtree = new TKDTree<Int_t, Float_t>(0,0,0,0);
  30   Int_t row =0;
  31   Int_t collumn =0;
  32   TStopwatch timer;
  33   timer.Start();
  34   row = 10;
  35   for (Int_t iloop=0;iloop<nloop;iloop++)
  36   for (Int_t index=1;index<npoints;index++){
  37     //row = TMath::Log2(index);
  38     //row=0;
  39    //  if (index< (16<<row)) row=0;
  40 //     for (; index>=(32<<row);row+=5);
  41 //     for (; index>=(2<<row);row++);
  42 //     collumn= index-(1<<row);
  43       TKDTree<Int_t, Float_t>::GetCoord(index,row,collumn);
  44     //
  45     //Int_t index2=kdtree->GetIndex(row,collumn);
  46     //printf("%d\t%d\t%d\t%d\n",index,index2,row,collumn);
  47     if (kdtree->GetIndex(row,collumn)!=index || collumn<0) {
  48       printf("Problem\n");
  49     }
  50   }
  51   timer.Stop();
  52   timer.Print();
  53 }
  54
  55 void TestBuild(Int_t npoints, Int_t bsize){
  56   Float_t *data0 =  new Float_t[npoints*2];
  57   Float_t *data[2];
  58   data[0] = &data0[0];
  59   data[1] = &data0[npoints];
  60   for (Int_t i=0;i<npoints;i++) {
  61     data[1][i]= gRandom->Rndm();
  62     data[0][i]= gRandom->Rndm();
  63     //data[1][i]= 0;
  64     //data[0][i]= -i;
  65   }
  66   Float_t dataf[2];
  67   TKDTreeIF *kdtree = new TKDTreeIF(npoints, 2, bsize, data);
  68   //kdtree->Build();
  69   for (Int_t i=0; i<npoints;i++){
  70     Int_t index = kdtree->fIndPoints[i];
  71     //printf("%d\t%d\t%f\n",i,index,data[0][index]);
  72   }
  73   for (Int_t i=0; i<kdtree->fNnodes;i++){
  74     //printf("%d\t%f\n",i,kdtree->fNodes[i].fValue);
  75   }
  76   Float_t sumiter = 0;
  77   for (Int_t i=0;i<npoints;i++){
  78     dataf[0] = data[0][i];
  79     dataf[1] = data[1][i];
  80     Int_t index=-1;
  81     Int_t iter =0;
  82     kdtree->FindPoint(dataf,index, iter);
  83     if (i!=index){
  84       printf("%d\t%d\t%f\t%f\n",i,index,dataf[0],data[0][index]);
  85     }
  86     sumiter+=iter;
  87   }
  88   printf("Mean iter = %f\n",float(sumiter)/float(npoints));
  89 }
  90
  91
  92 void TestSizeIF(Int_t npoints, Int_t nrows, Int_t nsec, Int_t mode, Int_t bsize)
  93 {
  94   //
  95   // test size to build kdtree
  96   //
  97   Int_t before =Mem();
  98   for (Int_t isec=0; isec<nsec;isec++)
  99     for (Int_t irow=0;irow<nrows;irow++){
 100       testkdtreeIF(npoints,0,mode,bsize);
 101     }
 102   Int_t after = Mem();
 103   printf("Memory usage %d\n",after-before);
 104 }
 105
 106
 107 void  testkdtreeIF(Int_t npoints, Int_t nloop, Int_t mode, Int_t bsize)
 108 {
 109   //
 110   // test speed and functionality of kdtree
 111   //
 112   Float_t rangey  = 100;
 113   Float_t rangez  = 100;
 114   Float_t drangey = 0.1;
 115   Float_t drangez = 0.1;
 116 //   Float_t rangey  = 20;
 117 //   Float_t rangez  = 250;
 118 //   Float_t drangey = 1;
 119 //   Float_t drangez = 1;
 120
 121   //
 122   Float_t *data0 =  new Float_t[npoints*2];
 123   Float_t *data[2];
 124   data[0] = &data0[0];
 125   data[1] = &data0[npoints];
 126   Int_t i;
 127   for (i=0; i<npoints; i++){
 128     data[0][i]          = gRandom->Uniform(-rangey, rangey);
 129     data[1][i]          = gRandom->Uniform(-rangez, rangez);
 130     //     data[i+npoints]  = TMath::Nint(gRandom->Uniform(-rangez, rangez)/10.)*10.;
 131     //printf("%d %f  %f\n", i, data[i], data[i+npoints]);
 132   }
 133   TStopwatch timerbuild;
 134   TKDTree<Int_t, Float_t> *kdtree = new TKDTree<Int_t, Float_t>(npoints,2,bsize,data);
 135    kdtree->Build();
 136    timerbuild.Stop();
 137    timerbuild.Print();
 138   TStopwatch timer;
 139
 140    Float_t countern=0;
 141    Float_t counteriter  = 0;
 142    Float_t counterfound = 0;
 143
 144
 145    if (mode ==2){
 146      if (nloop) timer.Start();
 147      Int_t res[npoints];
 148      Int_t nfound = 0;
 149      for (Int_t kloop = 0;kloop<nloop;kloop++){
 150        if (kloop==0){
 151          counteriter = 0;
 152          counterfound= 0;
 153          countern    = 0;
 154        }
 155        for (Int_t i=0;i<npoints;i++){
 156          Float_t point[2]={data[0][i],data[1][i]};
 157          Float_t delta[2]={drangey,drangez};
 158          Int_t iter  =0;
 159          nfound =0;
 160          Int_t bnode =0;
 161          //kdtree->FindBNode(point,delta, bnode);
 162          //continue;
 163          kdtree->FindInRangeA(point,delta,res,nfound,iter,bnode);
 164          if (kloop==0){
 165            //Bool_t isOK = kTRUE;
 166            Bool_t isOK = kFALSE;
 167            for (Int_t ipoint=0;ipoint<nfound;ipoint++)
 168              if (res[ipoint]==i) isOK =kTRUE;
 169            counteriter+=iter;
 170            counterfound+=nfound;
 171            if (isOK) {
 172              countern++;
 173            }else{
 174              printf("Bug\n");
 175            }
 176          }
 177        }
 178      }
 179      if (nloop){
 180        timer.Stop();
 181        timer.Print();
 182      }
 183    }
 184    delete [] data0;
 185
 186    counteriter/=npoints;
 187    counterfound/=npoints;
 188    if (nloop) printf("Find nearest point:\t%f\t%f\t%f\n",countern, counteriter, counterfound);
 189 }