new macros for testing the STAT package

[u/mrichter/AliRoot.git] / STAT / Macros / kNNSpeed.C

#include "TSystem.h"
#include "TFile.h"
#include "TTree.h"
#include "TH2I.h"
#include "TLegend.h"
#include "TRandom.h"
#include "TString.h"
#include "TGraph.h"
#include "TMarker.h"
#include "TStopwatch.h"
#include "../src/TKDInterpolator.h"
#include "../src/TKDTree.h"

const int ntimes = 500000; // times to repeate kNN search
TStopwatch timer;

Float_t kdTreeNN(Float_t **p, const int ndim = 5, const int np = 10000);
Float_t standAloneNN(Float_t **p, const int ndim = 5, const int np = 10000);

//______________________________________________________________
TTree *db = 0x0;
void kNNSpeed(const Int_t method = 0)
{
// Estimate kNN algorithm speed:
// method = 0 : run kDTree kNN search
// method = 1 : run kNN loop search

	const Int_t ndim  = 5;
	const Int_t nstat = 9;


	// generate test sample
	Float_t **p = new Float_t*[ntimes];
	for(int ip =0 ; ip<ntimes; ip++){
		p[ip] = new Float_t[ndim];
		for(int idim=0; idim<ndim; idim++) p[ip][idim] = gRandom->Uniform(-.5, .5);
	}

	Float_t time;
	Int_t stat;
	TLegend *leg = new TLegend(.7, .7, .9, .9);
	TH1 *h1 = new TH1I("h1", "", 100, 9.5, 20.5);
	h1->SetMaximum(1.E2);
	h1->SetMinimum(1.E-1);
	h1->GetXaxis()->SetTitle("Statistics");
	h1->GetYaxis()->SetTitle("CPU [#mus]");
	h1->SetStats(0);
	h1->Draw();
	TGraph **g=new TGraph*[ndim];
	for(int idim=2; idim<3/*ndim*/; idim++){
		g[idim] = new TGraph(nstat);
		g[idim]->SetMarkerStyle(20+idim);
		g[idim]->SetMarkerColor(idim+1);
		leg->AddEntry(g[idim], Form("%dD", idim+1), "pl");
		stat = 16384;//1024;
		for(int istat = 0; istat<nstat; istat++){
			time = kdTreeNN(p, idim+1, stat);
			//time = standAloneNN(p, idim+1, stat);
			g[idim]->SetPoint(istat, TMath::Log2(float(stat)), time);
			stat *= 2;
		}
		g[idim]->Draw("pl");
	}
	leg->Draw();
	
	for(int ip=0; ip<ntimes; ip++) delete [] p[ip];
	delete [] p;
}

//______________________________________________________________
Float_t kdTreeNN(Float_t **p, const int ndim, const int np)
{
	Float_t **x = new Float_t*[ndim];
	for(int idim =0 ; idim<ndim; idim++){
		x[idim] = new Float_t[np];
		for(int ip=0; ip<np; ip++) x[idim][ip] = gRandom->Gaus();
	}
	TKDTreeIF nnFinder(np, ndim, 1, x);
	nnFinder.MakeBoundaries();
	
	Int_t *index, fNN = 1;
	Float_t *d;
	timer.Start(kTRUE);
	for(int itime=0; itime<ntimes; itime++){
		nnFinder.FindNearestNeighbors(p[itime], fNN, index, d);
	}
	timer.Stop();
	Float_t time = 1.E6*timer.CpuTime()/float(ntimes);

	printf("np[%7d] nd[%d] time = %5.2f[mus]\n", np, ndim, time);
	
	// clean
	for(int idim=0; idim<ndim; idim++) delete [] x[idim];
	delete [] x;
	

	return time;
}

//______________________________________________________________
Float_t standAloneNN(Float_t **p, const int ndim, const int np)
{	
	// STAND ALONE
	const Int_t kNN = 1;
	Float_t ftmp, gtmp, dist;
	Int_t itmp, jtmp;
	Int_t   fkNN[kNN];
	Float_t fkNNdist[kNN];
	for(int i=0; i<kNN; i++) fkNNdist[i] = 9999.;

	Float_t **x = new Float_t*[ndim];
	for(int idim =0 ; idim<ndim; idim++){
		x[idim] = new Float_t[np];
		for(int ip=0; ip<np; ip++) x[idim][ip] = gRandom->Gaus();
	}
	Int_t npoints = np;
	
	// calculate
	Int_t ntimes = 100;
	timer.Start(kTRUE);
	for(int it=0; it<ntimes; it++){
		for(int idx=0; idx<npoints; idx++){
			// calculate distance in the L1 metric
			dist = 0.;
			for(int idim=0; idim<ndim; idim++) dist += TMath::Abs(p[it*10][idim] - x[idim][idx]);
			if(dist >= fkNNdist[kNN-1]) continue;
	
			//insert neighbor
			int iNN=0;
			while(dist >= fkNNdist[iNN]) iNN++;
			itmp = fkNN[iNN]; ftmp = fkNNdist[iNN];
			fkNN[iNN]     = idx;
			fkNNdist[iNN] = dist;
			for(int ireplace=iNN+1; ireplace<kNN; ireplace++){
				jtmp = fkNN[ireplace]; gtmp = fkNNdist[ireplace];
				fkNN[ireplace] = itmp; fkNNdist[ireplace] = ftmp;
				itmp = jtmp; ftmp = gtmp;
				if(ftmp == 9999.) break;
			}
		}
	}
	timer.Stop();
	
	Float_t time = timer.CpuTime()/float(ntimes);
	printf("np[%5d] nd[%d] time[%f] %6.2f\n", np, ndim, time, timer.CpuTime());

	return time;
}