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