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

#include <malloc.h>

#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/TKDPDF.h"
#include "../src/TKDTree.h"

void kNNTest(const int np = 10000, const int ndim = 2);
void kNNDraw(const Float_t *p, const int kNN=20);
void build(const Int_t ndim = 2, const Int_t nstat = 1000000);
Float_t Mem();


//______________________________________________________________
void kNNTest(const int np, const int ndim)
{
// Test speed and quality of nearest neighbors search.
// The results are compared with a simple loop search through the data.
// The macro should run in compiled mode.

	const int ntimes = 1000; // times to repeate kNN search
	const Int_t kNN = 1; // number of neighbors

	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);
	
	// 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->Gaus();
	}


	Int_t *index;
	Float_t *d;
	TStopwatch timer;
	timer.Start(kTRUE);
	for(int itime=0; itime<ntimes; itime++) nnFinder.FindNearestNeighbors(p[itime], kNN, index, d);
	timer.Stop();
	printf("kDTree NN calculation.\n");
	printf("cpu = %f [ms] real = %f [ms]\n", 1.E3*timer.CpuTime()/ntimes, 1.E3*timer.RealTime()/ntimes);
	printf("Points indexes in increasing order of their distance to the target point.\n");
	for(int i=0; i<kNN; i++){
		printf("%5d [%7.4f] ", index[i], d[i]);
		if(i%5==4) printf("\n");
	}
	printf("\n");
	
	// draw kNN
	TLegend *leg = new TLegend(.7, .7, .9, .9);
	leg->SetBorderSize(1);
	leg->SetHeader("NN finders");
	TH2 *h2 = new TH2I("h2", "", 100, -2., 2., 100, -2., 2.);
	h2->Draw(); h2->SetStats(kFALSE);
	TMarker *mp = new TMarker(p[ntimes-1][0], p[ntimes-1][1], 3);
	mp->SetMarkerColor(4);
	mp->Draw(); leg->AddEntry(mp, "Target", "p");
	TGraph *gKD = new TGraph(kNN);
	gKD->SetMarkerStyle(24);
	gKD->SetMarkerColor(2);
	gKD->SetMarkerSize(.5);
	for(int ip=0; ip<kNN; ip++) gKD->SetPoint(ip, x[0][index[ip]], x[1][index[ip]]);
	gKD->Draw("p"); leg->AddEntry(gKD, "kDTree", "p");
	

	// STAND ALONE
	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.;

	// calculate
	timer.Start(kTRUE);
	for(int idx=0; idx<np; idx++){
		// calculate distance in the L1 metric
		dist = 0.;
		for(int idim=0; idim<ndim; idim++) dist += TMath::Abs(p[ntimes-1][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();
	printf("\nStand Alone NN calculation.\n");
	printf("cpu = %f [ms] real = %f [ms]\n", 1.E3*timer.CpuTime(), 1.E3*timer.RealTime());
	printf("Points indexes in increasing order of their distance to the target point.\n");
	for(int i=0; i<kNN; i++){
		printf("%5d [%7.4f] ", fkNN[i], fkNNdist[i]);
		if(i%5==4) printf("\n");
	}
	printf("\n");	

	TGraph *gSA = new TGraph(kNN);
	gSA->SetMarkerStyle(24);
	for(int ip=0; ip<kNN; ip++) gSA->SetPoint(ip, x[0][fkNN[ip]], x[1][fkNN[ip]]);

	gSA->Draw("p"); leg->AddEntry(gSA, "Stand Alone", "p");
	leg->Draw();
	
	for(int ip=0; ip<ntimes; ip++) delete [] p[ip];
	delete [] p;
	for(int idim=0; idim<ndim; idim++) delete [] x[idim];
	delete [] x;
}


//______________________________________________________________
Float_t p[]={1.4, -.6}; //}{1.7, -.4};
void kNNDraw(const Float_t *p, const int kNN)
{
// Test memory consumption and draw "kNN" nearest neighbours of point "p".
// The distance (in the L1 metric) is encoded in the color code.

	const Int_t npoints = 10000;
	Float_t *data0 =  new Float_t[npoints*2];
  Float_t *data[2];
  data[0] = &data0[0];
  data[1] = &data0[npoints];
  for (Int_t i=0;i<npoints;i++) {
    data[1][i]= gRandom->Gaus();
    data[0][i]= gRandom->Gaus();
  }

	TKDPDF pdf(npoints, 2, 100, data);
	pdf.DrawNode(pdf.GetNodeIndex(p));
	
	TMarker *mp = new TMarker(p[0], p[1], 3);
	mp->SetMarkerColor(4);
	mp->Draw();
		
	Int_t *index, color;
	Float_t *d, d0, pknn[2];
	Float_t start = Mem();
	pdf.FindNearestNeighbors(p, kNN, index, d);
	Float_t end = Mem();
	printf("FindNearestNeighbors memory usage %fKB\n", end-start);
	d0 = d[kNN-1];
	TMarker *ma = new TMarker[kNN];
	for(int ip=0; ip<kNN; ip++){
		pdf.GetDataPoint(index[ip], pknn);
		color = 101 - Int_t(50. * d[ip]/d0);
		ma[ip].SetMarkerStyle(4);
		ma[ip].SetMarkerColor(color);
		ma[ip].DrawMarker(pknn[0], pknn[1]);
	}
}


//______________________________________________________________
Float_t Mem()
{
  // size in KB
  static struct mallinfo memdebug; 
  memdebug = mallinfo();
  return memdebug.uordblks/1024.;
}

//______________________________________________________________
void build(const Int_t ndim, const Int_t nstat)
{
// Build "nstat" data points in "ndim" dimensions taken from an
// uncorrelated 2D Gauss distribution.
	

	printf("build data ... \n");
	Double_t pntTrain[ndim];
	TFile *f = TFile::Open(Form("%dD_Gauss.root", ndim), "RECREATE");
	TTree *t = new TTree("db", "gauss database");
	for(int idim=0; idim<ndim; idim++) t->Branch(Form("x%d", idim), &pntTrain[idim], Form("x%d/D", idim));
	
	for(int istat=0; istat<nstat; istat++){
		for(int idim=0; idim<ndim; idim++) pntTrain[idim] = gRandom->Gaus();
		t->Fill();
	}
	f->cd();
	t->Write();
	f->Close();
	delete f;
}
Commit	Line	Data
	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
	16	void kNNTest(const int np = 10000, const int ndim = 2);
	17	void kNNDraw(const Float_t *p, const int kNN=20);
	18	void build(const Int_t ndim = 2, const Int_t nstat = 1000000);
	19	Float_t Mem();
	20
	21
	22	//______________________________________________________________
	23	void 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.E3timer.CpuTime()/ntimes, 1.E3timer.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.E3timer.CpuTime(), 1.E3timer.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	//______________________________________________________________
	133	Float_t p[]={1.4, -.6}; //}{1.7, -.4};
	134	void 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[npoints2];
	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	//______________________________________________________________
	175	Float_t Mem()
	176	{
	177	// size in KB
	178	static struct mallinfo memdebug;
	179	memdebug = mallinfo();
	180	return memdebug.uordblks/1024.;
	181	}
	182
	183	//______________________________________________________________
	184	void 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