STAT/Macros/kDTreeTest.C

   1 #include <malloc.h>
   2 #include "TMatrixD.h"
   3 #include "TRandom.h"
   4 #include "TGraph.h"
   5 #include "TStopwatch.h"
   6 #include "../src/TKDTree.h"
   7
   8
   9 Float_t Mem()
  10 {
  11   // size in KB
  12   static struct mallinfo memdebug;
  13   memdebug = mallinfo();
  14   return memdebug.uordblks/1024.;
  15 }
  16
  17
  18 void TestBuild(const Int_t npoints = 1000000, const Int_t bsize = 100);
  19 void TestSpeed(const Int_t npower2 = 20, const Int_t bsize = 10);
  20
  21 void testindexes(Int_t nloop, Int_t npoints);
  22 void  testkdtreeIF(Int_t npoints=1000, Int_t bsize=9, Int_t nloop=1000, Int_t mode = 2);
  23 void TestSizeIF(Int_t npoints=1000, Int_t nrows = 150, Int_t nsec=18, Int_t mode = 1, Int_t bsize=9);
  24
  25 //______________________________________________________________________
  26 void kDTreeTest()
  27 {
  28         printf("\n\tTesting kDTree memory usage ...\n");
  29         TestBuild();
  30
  31         printf("\n\tTesting kDTree speed ...\n");
  32         TestSpeed();
  33 }
  34
  35 //______________________________________________________________________
  36 void TestBuild(const Int_t npoints, const Int_t bsize){
  37   Float_t *data0 =  new Float_t[npoints*2];
  38   Float_t *data[2];
  39   data[0] = &data0[0];
  40   data[1] = &data0[npoints];
  41   for (Int_t i=0;i<npoints;i++) {
  42     data[1][i]= gRandom->Rndm();
  43     data[0][i]= gRandom->Rndm();
  44   }
  45   Float_t before =Mem();
  46   TKDTreeIF *kdtree = new TKDTreeIF(npoints, 2, bsize, data);
  47         Float_t after = Mem();
  48   printf("Memory usage %f KB\n",after-before);
  49         delete kdtree;
  50   Float_t end = Mem();
  51   printf("Memory leak %f KB\n", end-before);
  52         return;
  53 }
  54
  55 //______________________________________________________________________
  56 void TestSpeed(const Int_t npower2, const Int_t bsize)
  57 {
  58         if(npower2 < 10){
  59                 printf("Please specify a power of 2 greater than 10\n");
  60                 return;
  61         }
  62
  63         Int_t npoints = Int_t(pow(2., npower2))*bsize;
  64         Float_t *data0 =  new Float_t[npoints*2];
  65   Float_t *data[2];
  66   data[0] = &data0[0];
  67   data[1] = &data0[npoints];
  68   for (Int_t i=0;i<npoints;i++) {
  69     data[1][i]= gRandom->Rndm();
  70     data[0][i]= gRandom->Rndm();
  71   }
  72
  73         TGraph *g = new TGraph(npower2-10);
  74         g->SetMarkerStyle(7);
  75         TStopwatch timer;
  76         Int_t tpoints;
  77   TKDTreeIF *kdtree = 0x0;
  78         for(int i=10; i<npower2; i++){
  79                 tpoints = Int_t(pow(2., i))*bsize;
  80                 timer.Start(kTRUE);
  81                 kdtree = new TKDTreeIF(tpoints, 2, bsize, data);
  82                 timer.Stop();
  83                 g->SetPoint(i-10, i, timer.CpuTime());
  84                 printf("npoints [%d] nodes [%d] cpu time %f [s]\n", tpoints, kdtree->GetNNodes(), timer.CpuTime());
  85                 //timer.Print("u");
  86                 delete kdtree;
  87         }
  88         g->Draw("apl");
  89         return;
  90 }
  91
  92 //______________________________________________________________________
  93 void TestSizeIF(Int_t npoints, Int_t nrows, Int_t nsec, Int_t mode, Int_t bsize)
  94 {
  95   //
  96   // test size to build kdtree
  97   //
  98   Float_t before =Mem();
  99   for (Int_t isec=0; isec<nsec;isec++)
 100     for (Int_t irow=0;irow<nrows;irow++){
 101       testkdtreeIF(npoints,0,mode,bsize);
 102     }
 103   Float_t after = Mem();
 104   printf("Memory usage %f\n",after-before);
 105 }
 106
 107
 108 void testindexes(Int_t nloop, Int_t npoints){
 109   //
 110   // test indexing
 111   //
 112   TKDTree<Int_t, Float_t> *kdtree = new TKDTree<Int_t, Float_t>(0,0,0,0);
 113   Int_t row =0;
 114   Int_t collumn =0;
 115   TStopwatch timer;
 116   timer.Start();
 117   row = 10;
 118   for (Int_t iloop=0;iloop<nloop;iloop++)
 119   for (Int_t index=1;index<npoints;index++){
 120     //row = TMath::Log2(index);
 121     //row=0;
 122    //  if (index< (16<<row)) row=0;
 123 //     for (; index>=(32<<row);row+=5);
 124 //     for (; index>=(2<<row);row++);
 125 //     collumn= index-(1<<row);
 126       TKDTree<Int_t, Float_t>::GetCoord(index,row,collumn);
 127     //
 128     //Int_t index2=kdtree->GetIndex(row,collumn);
 129     //printf("%d\t%d\t%d\t%d\n",index,index2,row,collumn);
 130     if (kdtree->GetIndex(row,collumn)!=index || collumn<0) {
 131       printf("Problem\n");
 132     }
 133   }
 134   timer.Stop();
 135   timer.Print();
 136 }
 137
 138
 139 void  testkdtreeIF(Int_t npoints, Int_t bsize, Int_t nloop, Int_t mode)
 140 {
 141 //
 142 // Test speed and functionality of 2D kdtree.
 143 // Input parametrs:
 144 // npoints - number of data points
 145 // bsize   - bucket size
 146 // nloop   - number of loops
 147 // mode    - tasks to be performed by the kdTree
 148 //         - 0  : time building the tree
 149 //
 150
 151
 152   Float_t rangey  = 100;
 153   Float_t rangez  = 100;
 154   Float_t drangey = 0.1;
 155   Float_t drangez = 0.1;
 156
 157   //
 158   Float_t *data0 =  new Float_t[npoints*2];
 159   Float_t *data[2];
 160   data[0] = &data0[0];
 161   data[1] = &data0[npoints];
 162   Int_t i;
 163   for (i=0; i<npoints; i++){
 164     data[0][i]          = gRandom->Uniform(-rangey, rangey);
 165     data[1][i]          = gRandom->Uniform(-rangez, rangez);
 166   }
 167   TStopwatch timer;
 168
 169         // check time build
 170         printf("building kdTree ...\n");
 171         timer.Start(kTRUE);
 172   TKDTreeIF *kdtree = new TKDTreeIF(npoints, 2, bsize, data);
 173   timer.Stop();
 174   timer.Print();
 175         if(mode == 0) return;
 176
 177
 178
 179         Float_t countern=0;
 180         Float_t counteriter  = 0;
 181         Float_t counterfound = 0;
 182
 183
 184         if (mode ==2){
 185                 if (nloop) timer.Start(kTRUE);
 186                 Int_t res[npoints];
 187                 Int_t nfound = 0;
 188                 for (Int_t kloop = 0;kloop<nloop;kloop++){
 189                         if (kloop==0){
 190                                 counteriter = 0;
 191                                 counterfound= 0;
 192                                 countern    = 0;
 193                         }
 194                         for (Int_t i=0;i<npoints;i++){
 195                                 Float_t point[2]={data[0][i],data[1][i]};
 196                                 Float_t delta[2]={drangey,drangez};
 197                                 Int_t iter  =0;
 198                                 nfound =0;
 199                                 Int_t bnode =0;
 200                                 //kdtree->FindBNode(point,delta, bnode);
 201                                 //continue;
 202                                 kdtree->FindInRangeA(point,delta,res,nfound,iter,bnode);
 203                                 if (kloop==0){
 204                                         //Bool_t isOK = kTRUE;
 205                                         Bool_t isOK = kFALSE;
 206                                         for (Int_t ipoint=0;ipoint<nfound;ipoint++)
 207                                                 if (res[ipoint]==i) isOK =kTRUE;
 208                                         counteriter+=iter;
 209                                         counterfound+=nfound;
 210                                         if (isOK) {
 211                                                 countern++;
 212                                         }else{
 213                                                 printf("Bug\n");
 214                                         }
 215                                 }
 216                         }
 217                 }
 218
 219                 if (nloop){
 220                         timer.Stop();
 221                         timer.Print();
 222                 }
 223         }
 224         delete [] data0;
 225
 226         counteriter/=npoints;
 227         counterfound/=npoints;
 228         if (nloop) printf("Find nearest point:\t%f\t%f\t%f\n",countern, counteriter, counterfound);
 229 }