new macros for testing the STAT package
[u/mrichter/AliRoot.git] / 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 }