changes by Astrid commited
[u/mrichter/AliRoot.git] / STAT / Macros / kNNSpeed.C
CommitLineData
8a8f8274 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
14const int ntimes = 500000; // times to repeate kNN search
15TStopwatch timer;
16
17Float_t kdTreeNN(Float_t **p, const int ndim = 5, const int np = 10000);
18Float_t standAloneNN(Float_t **p, const int ndim = 5, const int np = 10000);
19
20//______________________________________________________________
21TTree *db = 0x0;
22void 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//______________________________________________________________
71Float_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//______________________________________________________________
101Float_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