1. Removed obsolete test functions
[u/mrichter/AliRoot.git] / STAT / Macros / kNNTest.C
CommitLineData
8a8f8274 1#include <malloc.h>
2
3#include "TSystem.h"
4#include "TFile.h"
5#include "TTree.h"
6#include "TH2I.h"
7#include "TLegend.h"
8#include "TRandom.h"
9#include "TString.h"
10#include "TGraph.h"
11#include "TMarker.h"
12#include "TStopwatch.h"
13#include "../src/TKDPDF.h"
14#include "../src/TKDTree.h"
15
16void kNNTest(const int np = 10000, const int ndim = 2);
17void kNNDraw(const Float_t *p, const int kNN=20);
18void build(const Int_t ndim = 2, const Int_t nstat = 1000000);
19Float_t Mem();
20
21
22//______________________________________________________________
23void kNNTest(const int np, const int ndim)
24{
25// Test speed and quality of nearest neighbors search.
26// The results are compared with a simple loop search through the data.
27// The macro should run in compiled mode.
28
29 const int ntimes = 1000; // times to repeate kNN search
30 const Int_t kNN = 1; // number of neighbors
31
32 Float_t **x = new Float_t*[ndim];
33 for(int idim =0 ; idim<ndim; idim++){
34 x[idim] = new Float_t[np];
35 for(int ip=0; ip<np; ip++) x[idim][ip] = gRandom->Gaus();
36 }
37 TKDTreeIF nnFinder(np, ndim, 1, x);
38
39 // generate test sample
40 Float_t **p = new Float_t*[ntimes];
41 for(int ip =0 ; ip<ntimes; ip++){
42 p[ip] = new Float_t[ndim];
43 for(int idim=0; idim<ndim; idim++) p[ip][idim] = gRandom->Gaus();
44 }
45
46
47 Int_t *index;
48 Float_t *d;
49 TStopwatch timer;
50 timer.Start(kTRUE);
51 for(int itime=0; itime<ntimes; itime++) nnFinder.FindNearestNeighbors(p[itime], kNN, index, d);
52 timer.Stop();
53 printf("kDTree NN calculation.\n");
54 printf("cpu = %f [ms] real = %f [ms]\n", 1.E3*timer.CpuTime()/ntimes, 1.E3*timer.RealTime()/ntimes);
55 printf("Points indexes in increasing order of their distance to the target point.\n");
56 for(int i=0; i<kNN; i++){
57 printf("%5d [%7.4f] ", index[i], d[i]);
58 if(i%5==4) printf("\n");
59 }
60 printf("\n");
61
62 // draw kNN
63 TLegend *leg = new TLegend(.7, .7, .9, .9);
64 leg->SetBorderSize(1);
65 leg->SetHeader("NN finders");
66 TH2 *h2 = new TH2I("h2", "", 100, -2., 2., 100, -2., 2.);
67 h2->Draw(); h2->SetStats(kFALSE);
68 TMarker *mp = new TMarker(p[ntimes-1][0], p[ntimes-1][1], 3);
69 mp->SetMarkerColor(4);
70 mp->Draw(); leg->AddEntry(mp, "Target", "p");
71 TGraph *gKD = new TGraph(kNN);
72 gKD->SetMarkerStyle(24);
73 gKD->SetMarkerColor(2);
74 gKD->SetMarkerSize(.5);
75 for(int ip=0; ip<kNN; ip++) gKD->SetPoint(ip, x[0][index[ip]], x[1][index[ip]]);
76 gKD->Draw("p"); leg->AddEntry(gKD, "kDTree", "p");
77
78
79
80 // STAND ALONE
81 Float_t ftmp, gtmp, dist;
82 Int_t itmp, jtmp;
83 Int_t fkNN[kNN];
84 Float_t fkNNdist[kNN];
85 for(int i=0; i<kNN; i++) fkNNdist[i] = 9999.;
86
87 // calculate
88 timer.Start(kTRUE);
89 for(int idx=0; idx<np; idx++){
90 // calculate distance in the L1 metric
91 dist = 0.;
92 for(int idim=0; idim<ndim; idim++) dist += TMath::Abs(p[ntimes-1][idim] - x[idim][idx]);
93 if(dist >= fkNNdist[kNN-1]) continue;
94
95 //insert neighbor
96 int iNN=0;
97 while(dist >= fkNNdist[iNN]) iNN++;
98 itmp = fkNN[iNN]; ftmp = fkNNdist[iNN];
99 fkNN[iNN] = idx;
100 fkNNdist[iNN] = dist;
101 for(int ireplace=iNN+1; ireplace<kNN; ireplace++){
102 jtmp = fkNN[ireplace]; gtmp = fkNNdist[ireplace];
103 fkNN[ireplace] = itmp; fkNNdist[ireplace] = ftmp;
104 itmp = jtmp; ftmp = gtmp;
105 if(ftmp == 9999.) break;
106 }
107 }
108 timer.Stop();
109 printf("\nStand Alone NN calculation.\n");
110 printf("cpu = %f [ms] real = %f [ms]\n", 1.E3*timer.CpuTime(), 1.E3*timer.RealTime());
111 printf("Points indexes in increasing order of their distance to the target point.\n");
112 for(int i=0; i<kNN; i++){
113 printf("%5d [%7.4f] ", fkNN[i], fkNNdist[i]);
114 if(i%5==4) printf("\n");
115 }
116 printf("\n");
117
118 TGraph *gSA = new TGraph(kNN);
119 gSA->SetMarkerStyle(24);
120 for(int ip=0; ip<kNN; ip++) gSA->SetPoint(ip, x[0][fkNN[ip]], x[1][fkNN[ip]]);
121
122 gSA->Draw("p"); leg->AddEntry(gSA, "Stand Alone", "p");
123 leg->Draw();
124
125 for(int ip=0; ip<ntimes; ip++) delete [] p[ip];
126 delete [] p;
127 for(int idim=0; idim<ndim; idim++) delete [] x[idim];
128 delete [] x;
129}
130
131
132//______________________________________________________________
133Float_t p[]={1.4, -.6}; //}{1.7, -.4};
134void kNNDraw(const Float_t *p, const int kNN)
135{
136// Test memory consumption and draw "kNN" nearest neighbours of point "p".
137// The distance (in the L1 metric) is encoded in the color code.
138
139 const Int_t npoints = 10000;
140 Float_t *data0 = new Float_t[npoints*2];
141 Float_t *data[2];
142 data[0] = &data0[0];
143 data[1] = &data0[npoints];
144 for (Int_t i=0;i<npoints;i++) {
145 data[1][i]= gRandom->Gaus();
146 data[0][i]= gRandom->Gaus();
147 }
148
149 TKDPDF pdf(npoints, 2, 100, data);
150 pdf.DrawNode(pdf.GetNodeIndex(p));
151
152 TMarker *mp = new TMarker(p[0], p[1], 3);
153 mp->SetMarkerColor(4);
154 mp->Draw();
155
156 Int_t *index, color;
157 Float_t *d, d0, pknn[2];
158 Float_t start = Mem();
159 pdf.FindNearestNeighbors(p, kNN, index, d);
160 Float_t end = Mem();
161 printf("FindNearestNeighbors memory usage %fKB\n", end-start);
162 d0 = d[kNN-1];
163 TMarker *ma = new TMarker[kNN];
164 for(int ip=0; ip<kNN; ip++){
165 pdf.GetDataPoint(index[ip], pknn);
166 color = 101 - Int_t(50. * d[ip]/d0);
167 ma[ip].SetMarkerStyle(4);
168 ma[ip].SetMarkerColor(color);
169 ma[ip].DrawMarker(pknn[0], pknn[1]);
170 }
171}
172
173
174//______________________________________________________________
175Float_t Mem()
176{
177 // size in KB
178 static struct mallinfo memdebug;
179 memdebug = mallinfo();
180 return memdebug.uordblks/1024.;
181}
182
183//______________________________________________________________
184void build(const Int_t ndim, const Int_t nstat)
185{
186// Build "nstat" data points in "ndim" dimensions taken from an
187// uncorrelated 2D Gauss distribution.
188
189
190 printf("build data ... \n");
191 Double_t pntTrain[ndim];
192 TFile *f = TFile::Open(Form("%dD_Gauss.root", ndim), "RECREATE");
193 TTree *t = new TTree("db", "gauss database");
194 for(int idim=0; idim<ndim; idim++) t->Branch(Form("x%d", idim), &pntTrain[idim], Form("x%d/D", idim));
195
196 for(int istat=0; istat<nstat; istat++){
197 for(int idim=0; idim<ndim; idim++) pntTrain[idim] = gRandom->Gaus();
198 t->Fill();
199 }
200 f->cd();
201 t->Write();
202 f->Close();
203 delete f;
204}
205