changes by Astrid commited
[u/mrichter/AliRoot.git] / STAT / Macros / kNNSpeed.C
1 #include "TSystem.h"
2 #include "TFile.h"
3 #include "TTree.h"
4 #include "TH2I.h"
5 #include "TLegend.h"
6 #include "TRandom.h"
7 #include "TString.h"
8 #include "TGraph.h"
9 #include "TMarker.h"
10 #include "TStopwatch.h"
11 #include "../src/TKDInterpolator.h"
12 #include "../src/TKDTree.h"
13
14 const int ntimes = 500000; // times to repeate kNN search
15 TStopwatch timer;
16
17 Float_t kdTreeNN(Float_t **p, const int ndim = 5, const int np = 10000);
18 Float_t standAloneNN(Float_t **p, const int ndim = 5, const int np = 10000);
19
20 //______________________________________________________________
21 TTree *db = 0x0;
22 void kNNSpeed(const Int_t method = 0)
23 {
24 // Estimate kNN algorithm speed:
25 // method = 0 : run kDTree kNN search
26 // method = 1 : run kNN loop search
27
28         const Int_t ndim  = 5;
29         const Int_t nstat = 9;
30
31
32         // generate test sample
33         Float_t **p = new Float_t*[ntimes];
34         for(int ip =0 ; ip<ntimes; ip++){
35                 p[ip] = new Float_t[ndim];
36                 for(int idim=0; idim<ndim; idim++) p[ip][idim] = gRandom->Uniform(-.5, .5);
37         }
38
39         Float_t time;
40         Int_t stat;
41         TLegend *leg = new TLegend(.7, .7, .9, .9);
42         TH1 *h1 = new TH1I("h1", "", 100, 9.5, 20.5);
43         h1->SetMaximum(1.E2);
44         h1->SetMinimum(1.E-1);
45         h1->GetXaxis()->SetTitle("Statistics");
46         h1->GetYaxis()->SetTitle("CPU [#mus]");
47         h1->SetStats(0);
48         h1->Draw();
49         TGraph **g=new TGraph*[ndim];
50         for(int idim=2; idim<3/*ndim*/; idim++){
51                 g[idim] = new TGraph(nstat);
52                 g[idim]->SetMarkerStyle(20+idim);
53                 g[idim]->SetMarkerColor(idim+1);
54                 leg->AddEntry(g[idim], Form("%dD", idim+1), "pl");
55                 stat = 16384;//1024;
56                 for(int istat = 0; istat<nstat; istat++){
57                         time = kdTreeNN(p, idim+1, stat);
58                         //time = standAloneNN(p, idim+1, stat);
59                         g[idim]->SetPoint(istat, TMath::Log2(float(stat)), time);
60                         stat *= 2;
61                 }
62                 g[idim]->Draw("pl");
63         }
64         leg->Draw();
65         
66         for(int ip=0; ip<ntimes; ip++) delete [] p[ip];
67         delete [] p;
68 }
69
70 //______________________________________________________________
71 Float_t kdTreeNN(Float_t **p, const int ndim, const int np)
72 {
73         Float_t **x = new Float_t*[ndim];
74         for(int idim =0 ; idim<ndim; idim++){
75                 x[idim] = new Float_t[np];
76                 for(int ip=0; ip<np; ip++) x[idim][ip] = gRandom->Gaus();
77         }
78         TKDTreeIF nnFinder(np, ndim, 1, x);
79         nnFinder.MakeBoundaries();
80         
81         Int_t *index, fNN = 1;
82         Float_t *d;
83         timer.Start(kTRUE);
84         for(int itime=0; itime<ntimes; itime++){
85                 nnFinder.FindNearestNeighbors(p[itime], fNN, index, d);
86         }
87         timer.Stop();
88         Float_t time = 1.E6*timer.CpuTime()/float(ntimes);
89
90         printf("np[%7d] nd[%d] time = %5.2f[mus]\n", np, ndim, time);
91         
92         // clean
93         for(int idim=0; idim<ndim; idim++) delete [] x[idim];
94         delete [] x;
95         
96
97         return time;
98 }
99
100 //______________________________________________________________
101 Float_t standAloneNN(Float_t **p, const int ndim, const int np)
102 {       
103         // STAND ALONE
104         const Int_t kNN = 1;
105         Float_t ftmp, gtmp, dist;
106         Int_t itmp, jtmp;
107         Int_t   fkNN[kNN];
108         Float_t fkNNdist[kNN];
109         for(int i=0; i<kNN; i++) fkNNdist[i] = 9999.;
110
111         Float_t **x = new Float_t*[ndim];
112         for(int idim =0 ; idim<ndim; idim++){
113                 x[idim] = new Float_t[np];
114                 for(int ip=0; ip<np; ip++) x[idim][ip] = gRandom->Gaus();
115         }
116         Int_t npoints = np;
117         
118         // calculate
119         Int_t ntimes = 100;
120         timer.Start(kTRUE);
121         for(int it=0; it<ntimes; it++){
122                 for(int idx=0; idx<npoints; idx++){
123                         // calculate distance in the L1 metric
124                         dist = 0.;
125                         for(int idim=0; idim<ndim; idim++) dist += TMath::Abs(p[it*10][idim] - x[idim][idx]);
126                         if(dist >= fkNNdist[kNN-1]) continue;
127         
128                         //insert neighbor
129                         int iNN=0;
130                         while(dist >= fkNNdist[iNN]) iNN++;
131                         itmp = fkNN[iNN]; ftmp = fkNNdist[iNN];
132                         fkNN[iNN]     = idx;
133                         fkNNdist[iNN] = dist;
134                         for(int ireplace=iNN+1; ireplace<kNN; ireplace++){
135                                 jtmp = fkNN[ireplace]; gtmp = fkNNdist[ireplace];
136                                 fkNN[ireplace] = itmp; fkNNdist[ireplace] = ftmp;
137                                 itmp = jtmp; ftmp = gtmp;
138                                 if(ftmp == 9999.) break;
139                         }
140                 }
141         }
142         timer.Stop();
143         
144         Float_t time = timer.CpuTime()/float(ntimes);
145         printf("np[%5d] nd[%d] time[%f] %6.2f\n", np, ndim, time, timer.CpuTime());
146
147         return time;
148 }
149