--- /dev/null
+void interpolCompare(const Int_t method = 0, const int nstat = 40000)
+{
+// Compare interpolation methods.
+// method = 0 : Select COG method
+// method = 1 : Select INT method
+
+ const Float_t alpha = 20.;
+ const int bucket = nstat/100; //400;
+ const int nevals = 100;
+ Float_t *data[1];
+ Float_t data1[nstat];
+ for(int istat=0; istat<nstat; istat++) data1[istat] = gRandom->Gaus();
+ data[0] = &data1[0];
+
+ TLegend *leg = new TLegend(0.7157258,0.8280255,0.9637097,0.9596603);
+ leg->SetBorderSize(1);
+ TF1 *f=new TF1("f1", "gaus(0);x;f(x)", -5., 5.);
+ f->SetLineWidth(1);
+ f->SetLineColor(2);
+ f->SetParameter(0, 1.);
+ f->SetParameter(1, 0.);
+ f->SetParameter(2, 1.);
+ f->SetParameter(0, 1./f->Integral(-10., 10.));
+ f->Draw(); leg->AddEntry(f, "model", "l");
+
+
+ Float_t cog, val, val_err;
+ Double_t x, chi2, chi2pdf, chi2int, result, err;
+ Int_t npPDF, npINT;
+ TKDPDF pdf(nstat, 1, bucket, data);
+ pdf.SetInterpolationMethod(method);
+ pdf.SetAlpha(alpha);
+
+ TKDInterpolator in(1, pdf.GetNTNodes());
+ in.SetInterpolationMethod(method);
+ in.SetAlpha(alpha);
+
+ // change to the reprezentation which is closer to the local
+ // interpolating polynom.
+ Double_t fx;
+ Int_t inode = 0;
+ TKDNodeInfo *node = 0x0;
+ while(node = pdf.GetNodeInfo(inode)){
+ if(node->fVal[0] > 0.){
+ in.SetNode(inode, *node);
+ node = in.GetNodeInfo(inode);
+ fx = node->fVal[0];
+ node->fVal[0] = TMath::Log(fx);
+ node->fVal[1] = node->fVal[1]/fx;
+ }
+ inode++;
+ }
+
+ // Do preprocessing for INT
+ Float_t *pcog;
+ if(method) for(int icog=0; icog<pdf.GetNTNodes(); icog++){
+ pdf.GetCOGPoint(icog, pcog, val, val_err);
+ x = pcog[0];
+ pdf.Eval(&x, result, err);
+ in.Eval(&x, result, err);
+ }
+
+
+
+ // prepare drawing containers
+ TGraphErrors *gPDF = new TGraphErrors(nevals);
+ gPDF->SetMarkerStyle(7);
+ TGraph *gREP = new TGraph(nevals);
+ gREP->SetMarkerStyle(7);
+ gREP->SetMarkerColor(4);
+ gREP->SetLineColor(4);
+
+ TGraphErrors *gINT = new TGraphErrors(nevals);
+ gINT->SetMarkerStyle(24);
+ gINT->SetMarkerSize(.6);
+ gINT->SetMarkerColor(3);
+ gINT->SetLineColor(3);
+ TGraph *gREI = new TGraph(nevals);
+ gREI->SetMarkerStyle(24);
+ gREI->SetMarkerColor(4);
+ gREI->SetLineColor(4);
+
+ // Do the interpolation with PDF and Interpolator
+ Double_t dx = 10./nevals;
+ x = 5.-dx/2.;
+ chi2pdf = 0.; npPDF = 0;
+ chi2int = 0.; npINT = 0;
+ for(int ip=0; ip<nevals; ip++){
+ chi2 = pdf.Eval(&x, result, err);
+ gPDF->SetPoint(ip, x, result);
+ gPDF->SetPointError(ip, 0., err);
+ if(TMath::Abs(x)<3.){
+ chi2pdf += (result - f->Eval(x))*(result - f->Eval(x))/f->Eval(x)/f->Eval(x);
+ npPDF++;
+ }
+ gREP->SetPoint(ip, x, .1*err/TMath::Abs(result));
+ //printf("%2d x[%f] r[%f +- %f] chi2[%e]\n", ip, x, result, err, chi2);
+
+ chi2 = in.Eval(&x, result, err);
+ result = TMath::Exp(result);
+ err = err*result;
+ gINT->SetPoint(ip, x, result);
+ gINT->SetPointError(ip, 0., err);
+ if(TMath::Abs(x)<3.){
+ chi2int += (result - f->Eval(x))*(result - f->Eval(x))/f->Eval(x)/f->Eval(x);
+ npINT++;
+ }
+ gREI->SetPoint(ip, x, .1*err/TMath::Abs(result));
+ //printf("%2d x[%f] r[%f +- %f] chi2[%e]\n\n", ip, x, result, err, chi2);
+
+ x -= dx;
+ }
+ gINT->Draw("p"); leg->AddEntry(gINT, Form("Interpolator [%s]", method ? "INT" : "COG"), "pl");
+ gPDF->Draw("p"); leg->AddEntry(gPDF, Form("PDF [%s]", method ? "INT" : "COG"), "pl");
+ gREI->Draw("pl"); leg->AddEntry(gREI, "1% #epsilon Interpolator", "pl");
+ gREP->Draw("pl"); leg->AddEntry(gREP, "1% #epsilon PDF", "pl");
+
+ chi2pdf /= npPDF;
+ printf("PDF : chi2(%d) = %f\n", npPDF, chi2pdf);
+ chi2int /= npINT;
+ printf("INT : chi2(%d) = %f\n", npINT, chi2int);
+
+ leg->Draw();
+}
+
+
+
--- /dev/null
+// #include "TSystem.h"
+// #include "TFile.h"
+// #include "TTree.h"
+// #include "TProfile.h"
+// #include "TF1.h"
+// #include "TF2.h"
+// #include "TF3.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/TKDInterpolator.h"
+
+void interpolTest(const Int_t ndim = 1);
+Double_t interpolation(const int nstat, const int ndim, const int first);
+void build(const int ndim, const int npoints);
+
+TF1 *f = 0x0;
+
+//__________________________________________________________
+void interpolTest(const Int_t ndim)
+{
+// Test interpolator on a variable dimension (ndim<=5) uncorrelated
+// Gauss model. The macro displays the chi2 between interpolation and
+// model for various statistics of experimental data.
+//
+// Check the worker function interpolation() to fine tune the
+// algorithm.
+//
+// Attention:
+// The macro works in compiled mode !
+
+ // define model
+ switch(ndim){
+ case 1:
+ f=new TF1("f1", "gaus(0);x;f(x)", -5., 5.);
+ f->SetParameter(0, 1.);
+ f->SetParameter(1, 0.);
+ f->SetParameter(2, 1.);
+ f->SetParameter(0, 1./f->Integral(-5., 5.));
+ break;
+ case 2:
+ f=new TF2("f2", "[0]*exp(-0.5*((x-[1])/[2])**2)*exp(-0.5*((y-[3])/[4])**2)", -5., 5., -5., 5.);
+ f->SetParameter(0, 1.);
+ f->SetParameter(1, 0.);
+ f->SetParameter(2, 1.);
+ f->SetParameter(3, 0.);
+ f->SetParameter(4, 1.);
+ f->SetParameter(0, 1./f->Integral(-5., 5., -5., 5.));
+ break;
+ case 3:
+ case 4:
+ case 5:
+ f=new TF3("f3", "[0]*exp(-0.5*((x-[1])/[2])**2)*exp(-0.5*((y-[3])/[4])**2)*exp(-0.5*((z-[5])/[6])**2)", -5., 5., -5., 5., -5., 5.);
+ f->SetParameter(0, 1.);
+ f->SetParameter(1, 0.);
+ f->SetParameter(2, 1.);
+ f->SetParameter(3, 0.);
+ f->SetParameter(4, 1.);
+ f->SetParameter(5, 0.);
+ f->SetParameter(6, 1.);
+ f->SetParameter(0, 1./f->Integral(-5., 5., -5., 5., -5., 5.));
+ if(ndim>3) printf("Warning : chi2 results are unreliable !\n");
+ break;
+ default:
+ printf("%dD not supported in this macro.\n", ndim);
+ return;
+ }
+
+ Double_t chi2;
+ const Int_t nchi2 = 18;
+ const Int_t nfirst = 1;
+ //TProfile *hp = new TProfile("hp", "", 100, 1.E4, 1.E6);
+ //hp->SetMarkerStyle(24);
+ TGraph *gChi2 = new TGraph(nchi2);
+ gChi2->SetMarkerStyle(24);
+ Int_t ip = 0, nstat, first;
+ for(int ifirst=0; ifirst<nfirst; ifirst++){
+ first = 0;//Int_t(gRandom->Uniform(1.E4));
+ for(int i=2; i<10; i++){
+ nstat = 10000*i;
+ chi2 = interpolation(nstat, ndim, first);
+ gChi2->SetPoint(ip++, TMath::Log10(float(nstat)), chi2);
+ //hp->Fill(float(nstat), chi2);
+ }
+ for(int i=1; i<10; i++){
+ nstat = 100000*i;
+ chi2 = interpolation(nstat, ndim, first);
+ gChi2->SetPoint(ip++, TMath::Log10(float(nstat)), chi2);
+ //hp->Fill(float(nstat), chi2);
+ }
+ gChi2->Draw("apl");
+ }
+ //hp->Draw("pl");
+}
+
+void test()
+{
+ TClonesArray ar("TKDNodeInfo", 10);
+}
+
+//__________________________________________________________
+Double_t interpolation(const int nstat, const int ndim, const int first)
+{
+// Steer interpolation of a nD (n<=5) uncorrelated Gauss distribution.
+// The user is suppose to give the experimental statistics (nstat) the
+// dimension of the experimental point (ndim). The entry from which the
+// data base is being read is steered by "first".
+
+ const int bs = 400;
+ const int npoints = 100;
+ // switch on chi2 calculation between interpolation and model.
+ // Default calculates distance.
+ const Bool_t kCHI2 = kFALSE;
+ const Bool_t kINT = kFALSE; // switch on integral interpolator
+
+ // open data base
+ TString fn("5D_Gauss.root");
+ if(!gSystem->FindFile(".", fn)) build(5, 1000000);
+ TFile::Open("5D_Gauss.root");
+ TTree *t = (TTree*)gFile->Get("db");
+
+ // build interpolator
+ TString var = "x0"; for(int idim=1; idim<ndim; idim++) var += Form(":x%d", idim);
+ TKDPDF pdf(t, var.Data(), "", bs, nstat, first);
+ printf("\n\nFINISH BUILDING PDF\n\n");
+
+ Double_t fx;
+ Int_t inode = 0;
+ TKDNodeInfo *node = 0x0;
+ TKDInterpolator in(ndim, pdf.GetNTNodes());
+/* while(node = pdf.GetNodeInfo(inode)){
+ if(node->fVal[0] > 0.){
+ fx = node->fVal[0];
+ node->fVal[0] = TMath::Log(fx);
+ node->fVal[1] = node->fVal[1]/fx;
+ }
+ in.SetNode(inode, *node);
+ inode++;
+ }*/
+ //pdf.SetInterpolationMethod(kINT);
+ pdf.SetAlpha(2.);
+ //in.SetStore();
+
+
+ Double_t x[2], r, e, chi2;
+ TH2 *h2 = new TH2F("h2", "", 50, -4., 4., 50, -4., 4.);
+ TAxis *ax = h2->GetXaxis(), *ay = h2->GetYaxis();
+ for(int ix=1; ix<=ax->GetNbins(); ix++){
+ x[0] = ax->GetBinCenter(ix);
+ for(int iy=1; iy<=ay->GetNbins(); iy++){
+ x[1] = ay->GetBinCenter(iy);
+
+ chi2 = pdf.Eval(x, r, e);
+ printf("x[%2d] x[%2d] r[%f] e[%f] chi2[%f]\n", ix, iy, r, e, chi2);
+ h2->SetBinContent(ix, iy, r/*TMath::Exp(r)*/);
+ }
+ }
+ h2->Draw("lego2");
+ return;
+
+ // define testing grid
+ Double_t x[5], r, e;
+ Double_t dx[5] = {4., 4., 4., 4., 4.};
+ Double_t chi2, theor, result, err;
+ for(int idim=0; idim<ndim; idim++){
+ dx[idim] = (ndim<<2)/float(npoints);
+ x[idim] = -2.+dx[idim]/2.;
+ }
+
+
+ // setup the integral interpolator and do memory test
+ Float_t c[5], v, ve;
+ for(int ip=0; ip<in.GetNTNodes(); ip++){
+ in.GetCOGPoint(ip, c, v, ve);
+ for(int idim=0; idim<ndim; idim++) x[idim] = (Double_t)c[idim];
+ in.Eval(x, result, err);
+ }
+ printf("\nInterpolating %d data points in %dD ...\n", nstat, ndim);
+ printf("\tbucket size %d\n", bs);
+ printf("\tstarting interpolation ...\n");
+ return;
+
+
+ // evaluate relative error
+ for(int idim=0; idim<ndim; idim++) x[idim] = 0.;
+ in.Eval(x, result, err);
+ Double_t threshold = err/result;
+
+ // starting interpolation over the testing grid
+ chi2 = 0.; Int_t nsample = 0;
+ x[4] = -2.+dx[4]/2.;
+ do{
+ x[3] = -2.+dx[3]/2.;
+ do{
+ x[2] = -2.+dx[2]/2.;
+ do{
+ x[1] = -2.+dx[1]/2.;
+ do{
+ x[0] = -2.+dx[0]/2.;
+ do{
+ in.Eval(x, result, err);
+ if(err/TMath::Abs(result) < 1.1*threshold){
+ theor = f->Eval(x[0], x[1], x[2]);
+ Double_t tmp = result - theor; tmp *= tmp;
+ chi2 += tmp/(kCHI2 ? err*err : theor);
+ nsample++;
+ }
+ } while((x[0] += dx[0]) < 2.);
+ } while((x[1] += dx[1]) < 2.);
+ } while((x[2] += dx[2]) < 2.);
+ } while((x[3] += dx[3]) < 2.);
+ } while((x[4] += dx[4]) < 2.);
+ chi2 /= float(nsample);
+ printf("\tinterpolation quality (chi2) %f\n", chi2);
+
+ return kCHI2 ? chi2 : TMath::Sqrt(chi2);
+}
+
+
+//___________________________________________________________
+void build(const int ndim, const int npoints)
+{
+ printf("Building DB ...\n");
+ Float_t data[ndim];
+ TFile::Open(Form("%dD_Gauss.root", ndim), "RECREATE");
+ TTree *t = new TTree("db", Form("%dD data base for kD statistics", ndim));
+ for(int idim=0; idim<ndim; idim++) t->Branch(Form("x%d", idim), &data[idim], Form("x%d/F", idim));
+
+ for (Int_t ip=0; ip<npoints; ip++){
+ for (Int_t id=0; id<ndim; id++) data[id]= gRandom->Gaus();
+ t->Fill();
+ }
+
+ t->Write();
+ gFile->Close();
+ delete gFile;
+}
--- /dev/null
+#include <malloc.h>
+#include "TMatrixD.h"
+#include "TRandom.h"
+#include "TGraph.h"
+#include "TStopwatch.h"
+#include "../src/TKDTree.h"
+
+
+Float_t Mem()
+{
+ // size in KB
+ static struct mallinfo memdebug;
+ memdebug = mallinfo();
+ return memdebug.uordblks/1024.;
+}
+
+
+void TestBuild(const Int_t npoints = 1000000, const Int_t bsize = 100);
+void TestSpeed(const Int_t npower2 = 20, const Int_t bsize = 10);
+
+void testindexes(Int_t nloop, Int_t npoints);
+void testkdtreeIF(Int_t npoints=1000, Int_t bsize=9, Int_t nloop=1000, Int_t mode = 2);
+void TestSizeIF(Int_t npoints=1000, Int_t nrows = 150, Int_t nsec=18, Int_t mode = 1, Int_t bsize=9);
+
+//______________________________________________________________________
+void kDTreeTest()
+{
+ printf("\n\tTesting kDTree memory usage ...\n");
+ TestBuild();
+
+ printf("\n\tTesting kDTree speed ...\n");
+ TestSpeed();
+}
+
+//______________________________________________________________________
+void TestBuild(const Int_t npoints, const Int_t bsize){
+ 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->Rndm();
+ data[0][i]= gRandom->Rndm();
+ }
+ Float_t before =Mem();
+ TKDTreeIF *kdtree = new TKDTreeIF(npoints, 2, bsize, data);
+ Float_t after = Mem();
+ printf("Memory usage %f KB\n",after-before);
+ delete kdtree;
+ Float_t end = Mem();
+ printf("Memory leak %f KB\n", end-before);
+ return;
+}
+
+//______________________________________________________________________
+void TestSpeed(const Int_t npower2, const Int_t bsize)
+{
+ if(npower2 < 10){
+ printf("Please specify a power of 2 greater than 10\n");
+ return;
+ }
+
+ Int_t npoints = Int_t(pow(2., npower2))*bsize;
+ 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->Rndm();
+ data[0][i]= gRandom->Rndm();
+ }
+
+ TGraph *g = new TGraph(npower2-10);
+ g->SetMarkerStyle(7);
+ TStopwatch timer;
+ Int_t tpoints;
+ TKDTreeIF *kdtree = 0x0;
+ for(int i=10; i<npower2; i++){
+ tpoints = Int_t(pow(2., i))*bsize;
+ timer.Start(kTRUE);
+ kdtree = new TKDTreeIF(tpoints, 2, bsize, data);
+ timer.Stop();
+ g->SetPoint(i-10, i, timer.CpuTime());
+ printf("npoints [%d] nodes [%d] cpu time %f [s]\n", tpoints, kdtree->GetNNodes(), timer.CpuTime());
+ //timer.Print("u");
+ delete kdtree;
+ }
+ g->Draw("apl");
+ return;
+}
+
+//______________________________________________________________________
+void TestSizeIF(Int_t npoints, Int_t nrows, Int_t nsec, Int_t mode, Int_t bsize)
+{
+ //
+ // test size to build kdtree
+ //
+ Float_t before =Mem();
+ for (Int_t isec=0; isec<nsec;isec++)
+ for (Int_t irow=0;irow<nrows;irow++){
+ testkdtreeIF(npoints,0,mode,bsize);
+ }
+ Float_t after = Mem();
+ printf("Memory usage %f\n",after-before);
+}
+
+
+void testindexes(Int_t nloop, Int_t npoints){
+ //
+ // test indexing
+ //
+ TKDTree<Int_t, Float_t> *kdtree = new TKDTree<Int_t, Float_t>(0,0,0,0);
+ Int_t row =0;
+ Int_t collumn =0;
+ TStopwatch timer;
+ timer.Start();
+ row = 10;
+ for (Int_t iloop=0;iloop<nloop;iloop++)
+ for (Int_t index=1;index<npoints;index++){
+ //row = TMath::Log2(index);
+ //row=0;
+ // if (index< (16<<row)) row=0;
+// for (; index>=(32<<row);row+=5);
+// for (; index>=(2<<row);row++);
+// collumn= index-(1<<row);
+ TKDTree<Int_t, Float_t>::GetCoord(index,row,collumn);
+ //
+ //Int_t index2=kdtree->GetIndex(row,collumn);
+ //printf("%d\t%d\t%d\t%d\n",index,index2,row,collumn);
+ if (kdtree->GetIndex(row,collumn)!=index || collumn<0) {
+ printf("Problem\n");
+ }
+ }
+ timer.Stop();
+ timer.Print();
+}
+
+
+void testkdtreeIF(Int_t npoints, Int_t bsize, Int_t nloop, Int_t mode)
+{
+//
+// Test speed and functionality of 2D kdtree.
+// Input parametrs:
+// npoints - number of data points
+// bsize - bucket size
+// nloop - number of loops
+// mode - tasks to be performed by the kdTree
+// - 0 : time building the tree
+//
+
+
+ Float_t rangey = 100;
+ Float_t rangez = 100;
+ Float_t drangey = 0.1;
+ Float_t drangez = 0.1;
+
+ //
+ Float_t *data0 = new Float_t[npoints*2];
+ Float_t *data[2];
+ data[0] = &data0[0];
+ data[1] = &data0[npoints];
+ Int_t i;
+ for (i=0; i<npoints; i++){
+ data[0][i] = gRandom->Uniform(-rangey, rangey);
+ data[1][i] = gRandom->Uniform(-rangez, rangez);
+ }
+ TStopwatch timer;
+
+ // check time build
+ printf("building kdTree ...\n");
+ timer.Start(kTRUE);
+ TKDTreeIF *kdtree = new TKDTreeIF(npoints, 2, bsize, data);
+ timer.Stop();
+ timer.Print();
+ if(mode == 0) return;
+
+
+
+ Float_t countern=0;
+ Float_t counteriter = 0;
+ Float_t counterfound = 0;
+
+
+ if (mode ==2){
+ if (nloop) timer.Start(kTRUE);
+ Int_t res[npoints];
+ Int_t nfound = 0;
+ for (Int_t kloop = 0;kloop<nloop;kloop++){
+ if (kloop==0){
+ counteriter = 0;
+ counterfound= 0;
+ countern = 0;
+ }
+ for (Int_t i=0;i<npoints;i++){
+ Float_t point[2]={data[0][i],data[1][i]};
+ Float_t delta[2]={drangey,drangez};
+ Int_t iter =0;
+ nfound =0;
+ Int_t bnode =0;
+ //kdtree->FindBNode(point,delta, bnode);
+ //continue;
+ kdtree->FindInRangeA(point,delta,res,nfound,iter,bnode);
+ if (kloop==0){
+ //Bool_t isOK = kTRUE;
+ Bool_t isOK = kFALSE;
+ for (Int_t ipoint=0;ipoint<nfound;ipoint++)
+ if (res[ipoint]==i) isOK =kTRUE;
+ counteriter+=iter;
+ counterfound+=nfound;
+ if (isOK) {
+ countern++;
+ }else{
+ printf("Bug\n");
+ }
+ }
+ }
+ }
+
+ if (nloop){
+ timer.Stop();
+ timer.Print();
+ }
+ }
+ delete [] data0;
+
+ counteriter/=npoints;
+ counterfound/=npoints;
+ if (nloop) printf("Find nearest point:\t%f\t%f\t%f\n",countern, counteriter, counterfound);
+}
--- /dev/null
+#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;
+}
+
--- /dev/null
+#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;
+}
+
--- /dev/null
+void pdfCompare(const Int_t method = 0, const int nstat = 40000)
+{
+// Compare interpolation methods.
+// method = 0 : Select COG method
+// method = 1 : Select INT method
+
+ const int bucket = nstat/100; //400;
+ const int nevals = 100;
+ Float_t *data[1];
+ Float_t data1[nstat];
+ for(int istat=0; istat<nstat; istat++) data1[istat] = gRandom->Gaus();
+ data[0] = &data1[0];
+
+ TLegend *leg = new TLegend(0.7157258,0.8280255,0.9637097,0.9596603);
+ leg->SetBorderSize(1);
+ TF1 *f=new TF1("f1", "gaus(0);x;f(x)", -5., 5.);
+ f->SetLineWidth(1);
+ f->SetLineColor(2);
+ f->SetParameter(0, 1.);
+ f->SetParameter(1, 0.);
+ f->SetParameter(2, 1.);
+ f->SetParameter(0, 1./f->Integral(-10., 10.));
+ f->Draw(); leg->AddEntry(f, "model", "l");
+
+
+ Float_t cog, val, val_err;
+ Double_t x, chi2, CHI2, result, err;
+ Int_t npoints;
+ TKDPDF pdf(nstat, 1, bucket, data);
+ pdf.SetInterpolationMethod(method);
+ pdf.SetStore();
+ pdf.SetAlpha(1.);
+ Float_t *bounds;// = in.GetBoundary(0);
+ Double_t dx = 10./nevals;
+
+ TGraph *gRE = new TGraph(npoints);
+ gRE->SetMarkerStyle(4);
+ gRE->SetMarkerColor(4);
+ gRE->SetLineColor(4);
+
+ // do a preevaluation for INT method on COG points
+ Float_t *pcog;
+ for(int icog=0; icog<pdf.GetNTNodes(); icog++){
+ pdf.GetCOGPoint(icog, pcog, val, val_err);
+ x = pcog[0];
+ chi2 = pdf.Eval(&x, result, err);
+ }
+ //pdf.GetStatus();
+
+ TGraphErrors *gINT = new TGraphErrors(npoints);
+ gINT->SetMarkerStyle(7);
+ x = 5.-dx/2.;
+ CHI2 = 0.; npoints = 0;
+ for(int ip=0; ip<nevals; ip++){
+ chi2 = pdf.Eval(&x, result, err);
+ gINT->SetPoint(ip, x, result);
+ gINT->SetPointError(ip, 0., err);
+ gRE->SetPoint(ip, x, .01*err/TMath::Abs(result));
+ if(TMath::Abs(x)<2.){
+ CHI2 += (result - f->Eval(x))*(result - f->Eval(x))/err/err;
+ npoints++;
+ }
+ x -= dx;
+ }
+ gINT->Draw("p"); leg->AddEntry(gINT, Form("interpolation [%s]", method ? "INT" : "COG"), "pl");
+
+ gRE->Draw("pl"); leg->AddEntry(gRE, "1% relative errors", "pl");
+ CHI2 /= npoints;
+ printf("chi2(%d) = %f\n", npoints, CHI2);
+
+ leg->Draw();
+}
+
+
+
--- /dev/null
+TStopwatch timer;
+const Int_t nt = 100;
+//___________________________________________________________
+void pdfIO(const Int_t method = 0)
+{
+// Test interpolator IO response for several data dimensions.
+// method = 0 : use COG interpolator
+// method = 1 : use INT interpolator
+
+ Float_t tw, tr;
+ TGraph *gw = new TGraph(5);
+ gw->SetMarkerStyle(24);gw->SetMarkerColor(2);
+ TGraph *gr = new TGraph(5);
+ gr->SetMarkerStyle(25);gr->SetMarkerColor(4);
+ for(int idim = 1; idim<6; idim++){
+ tw = interpolWrite(method, idim);
+ gw->SetPoint(idim-1, idim, tw);
+ tr = interpolRead(idim);
+ gr->SetPoint(idim-1, idim, tr);
+ }
+ gw->Draw("apl");
+ gr->Draw("pl");
+}
+
+//___________________________________________________________
+Float_t interpolWrite(const Int_t method, const Int_t ndim)
+{
+ if(!gSystem->FindFile(".", "5D_Gauss.root")) build(5, 1000000);
+ TFile::Open("5D_Gauss.root");
+ printf("\nWriting %dD interpolator ...\n", ndim);
+ TTree *db = (TTree*)gFile->Get("db");
+ TString var = "x0"; for(int idim=1; idim<ndim; idim++) var += Form(":x%d", idim);
+ printf("\tBuilding interpolator ...\n");
+ TKDPDF pdf(db, var.Data(), "", 400, 100000);
+ pdf.SetInterpolationMethod(method);
+ pdf.SetStore();
+
+
+ printf("\tSetting up the interpolator ...\n");
+ Float_t *c, val, v_err;
+ Double_t *p = new Double_t[ndim], res, r_err;
+ timer.Start(kTRUE);
+ for(int inode=0; inode<pdf.GetNTNodes(); inode++){
+ pdf.GetCOGPoint(inode, c, val, v_err);
+ for(int idim=0; idim<ndim; idim++) p[idim] = (Double_t)c[idim];
+// printf("Evaluate for node [%d] ", inode);
+// for(int idim=0; idim<ndim; idim++) printf("%5.3f ", p[idim]); printf("\n");
+ for(int it=0; it<nt; it++) pdf.Eval(p, res, r_err, kTRUE);
+// printf("R = %6.4f +- %6.4f\n", res, r_err);
+ }
+ timer.Stop();
+
+ Float_t time = 1.E3*timer.CpuTime()/float(nt);
+ printf("\tFinish interpolation in %6.2f [ms]\n", time);
+
+ printf("\tSaving interpolator ...\n");
+ TFile *fi = TFile::Open(Form("%dD_Interpolator.root", ndim), "RECREATE");
+ pdf.Write(Form("%dDgauss", ndim));
+ fi->Close();
+ delete fi;
+
+ delete [] p;
+
+ return time;
+}
+
+//___________________________________________________________
+Float_t interpolRead(const Int_t ndim)
+{
+ printf("Reading %dD interpolator ...\n", ndim);
+ TFile::Open(Form("%dD_Interpolator.root", ndim));
+ TKDPDF *pdf = (TKDPDF*)gFile->Get(Form("%dDgauss", ndim));
+ gFile->Close();
+
+ printf("\tDoing interpolation ...\n");
+ Float_t *c, v, ve;
+ Double_t *p = new Double_t[ndim], r, re;
+ timer.Start(kTRUE);
+ for(int ip=0; ip<pdf->GetNTNodes(); ip++){
+ pdf->GetCOGPoint(ip, c, v, ve);
+ for(int idim=0; idim<ndim; idim++) p[idim] = (Double_t)c[idim];
+ for(int it=0; it<nt; it++) pdf->Eval(p, r, re);
+ }
+ timer.Stop();
+ Float_t time = 1.E3*timer.CpuTime()/float(nt);
+ printf("\tFinish interpolation in %6.2f [ms]\n", time);
+ delete [] p;
+ return time;
+}
+
+//___________________________________________________________
+void build(const int ndim, const int npoints)
+{
+ printf("Building DB ...\n");
+ Float_t data[ndim];
+ TFile::Open(Form("%dD_Gauss.root", ndim), "RECREATE");
+ TTree *t = new TTree("db", Form("%dD data base for kD statistics", ndim));
+ for(int idim=0; idim<ndim; idim++) t->Branch(Form("x%d", idim), &data[idim], Form("x%d/F", idim));
+
+ for (Int_t ip=0; ip<npoints; ip++){
+ for (Int_t id=0; id<ndim; id++) data[id]= gRandom->Gaus();
+ t->Fill();
+ }
+
+ t->Write();
+ gFile->Close();
+ delete gFile;
+}
+
+
--- /dev/null
+#include <malloc.h>
+
+#include "TSystem.h"
+#include "TFile.h"
+#include "TTree.h"
+#include "TProfile.h"
+#include "TF1.h"
+#include "TF2.h"
+#include "TF3.h"
+#include "TLegend.h"
+#include "TRandom.h"
+#include "TString.h"
+#include "TGraph.h"
+#include "TMarker.h"
+#include "TStopwatch.h"
+#include "../src/TKDPDF.h"
+
+void pdfTest(const Int_t ndim = 1);
+Double_t interpolation(const int nstat, const int ndim, const int first);
+void build(const int ndim, const int npoints);
+Float_t Mem();
+
+TF1 *f = 0x0;
+
+//__________________________________________________________
+void pdfTest(const Int_t ndim)
+{
+// Test interpolator on a variable dimension (ndim<=5) uncorrelated
+// Gauss model. The macro displays the chi2 between interpolation and
+// model for various statistics of experimental data.
+//
+// Check the worker function interpolation() to fine tune the
+// algorithm.
+//
+// Attention:
+// The macro works in compiled mode !
+
+ // define model
+ switch(ndim){
+ case 1:
+ f=new TF1("f1", "gaus(0);x;f(x)", -5., 5.);
+ f->SetParameter(0, 1.);
+ f->SetParameter(1, 0.);
+ f->SetParameter(2, 1.);
+ f->SetParameter(0, 1./f->Integral(-5., 5.));
+ break;
+ case 2:
+ f=new TF2("f2", "[0]*exp(-0.5*((x-[1])/[2])**2)*exp(-0.5*((y-[3])/[4])**2)", -5., 5., -5., 5.);
+ f->SetParameter(0, 1.);
+ f->SetParameter(1, 0.);
+ f->SetParameter(2, 1.);
+ f->SetParameter(3, 0.);
+ f->SetParameter(4, 1.);
+ f->SetParameter(0, 1./f->Integral(-5., 5., -5., 5.));
+ break;
+ case 3:
+ case 4:
+ case 5:
+ f=new TF3("f3", "[0]*exp(-0.5*((x-[1])/[2])**2)*exp(-0.5*((y-[3])/[4])**2)*exp(-0.5*((z-[5])/[6])**2)", -5., 5., -5., 5., -5., 5.);
+ f->SetParameter(0, 1.);
+ f->SetParameter(1, 0.);
+ f->SetParameter(2, 1.);
+ f->SetParameter(3, 0.);
+ f->SetParameter(4, 1.);
+ f->SetParameter(5, 0.);
+ f->SetParameter(6, 1.);
+ f->SetParameter(0, 1./f->Integral(-5., 5., -5., 5., -5., 5.));
+ if(ndim>3) printf("Warning : chi2 results are unreliable !\n");
+ break;
+ default:
+ printf("%dD not supported in this macro.\n", ndim);
+ return;
+ }
+
+ Double_t chi2;
+ const Int_t nchi2 = 18;
+ const Int_t nfirst = 1;
+ //TProfile *hp = new TProfile("hp", "", 100, 1.E4, 1.E6);
+ //hp->SetMarkerStyle(24);
+ TGraph *gChi2 = new TGraph(nchi2);
+ gChi2->SetMarkerStyle(24);
+ Int_t ip = 0, nstat, first;
+ for(int ifirst=0; ifirst<nfirst; ifirst++){
+ first = 0;//Int_t(gRandom->Uniform(1.E4));
+ for(int i=2; i<10; i++){
+ nstat = 10000*i;
+ chi2 = interpolation(nstat, ndim, first);
+ gChi2->SetPoint(ip++, TMath::Log10(float(nstat)), chi2);
+ //hp->Fill(float(nstat), chi2);
+ }
+ for(int i=1; i<10; i++){
+ nstat = 100000*i;
+ chi2 = interpolation(nstat, ndim, first);
+ gChi2->SetPoint(ip++, TMath::Log10(float(nstat)), chi2);
+ //hp->Fill(float(nstat), chi2);
+ }
+ gChi2->Draw("apl");
+ }
+ //hp->Draw("pl");
+}
+
+
+//__________________________________________________________
+Double_t interpolation(const int nstat, const int ndim, const int first)
+{
+// Steer interpolation of a nD (n<=5) uncorrelated Gauss distribution.
+// The user is suppose to give the experimental statistics (nstat) the
+// dimension of the experimental point (ndim). The entry from which the
+// data base is being read is steered by "first".
+
+ const int bs = 400;
+ const int npoints = 100;
+ // switch on chi2 calculation between interpolation and model.
+ // Default calculates distance.
+ const Bool_t kCHI2 = kFALSE;
+ const Bool_t kINT = kTRUE; // switch on integral interpolator
+
+ // open data base
+ TString fn("5D_Gauss.root");
+ if(!gSystem->FindFile(".", fn)) build(5, 1000000);
+ TFile::Open("5D_Gauss.root");
+ TTree *t = (TTree*)gFile->Get("db");
+
+ // build interpolator
+ TString var = "x0"; for(int idim=1; idim<ndim; idim++) var += Form(":x%d", idim);
+ TKDPDF pdf(t, var.Data(), "", bs, nstat, first);
+ pdf.SetInterpolationMethod(kINT);
+ pdf.SetStore();
+
+ // define testing grid
+ Double_t x[5];
+ Double_t dx[5] = {4., 4., 4., 4., 4.};
+ Double_t chi2, theor, result, err;
+ for(int idim=0; idim<ndim; idim++){
+ dx[idim] = (ndim<<2)/float(npoints);
+ x[idim] = -2.+dx[idim]/2.;
+ }
+
+ // setup the integral interpolator and do memory test
+ Float_t start, stop;
+ start = Mem();
+ Float_t *c, v, ve;
+ for(int ip=0; ip<pdf.GetNTNodes(); ip++){
+ pdf.GetCOGPoint(ip, c, v, ve);
+ for(int idim=0; idim<ndim; idim++) x[idim] = (Double_t)c[idim];
+ pdf.Eval(x, result, err);
+ }
+ stop = Mem();
+ printf("\nInterpolating %d data points in %dD ...\n", nstat, ndim);
+ printf("\tbucket size %d\n", bs);
+ printf("\tmemory usage in Eval() %6.4fKB\n", stop-start);
+ printf("\tstarting interpolation ...\n");
+ //return stop-start;
+
+ // evaluate relative error
+ for(int idim=0; idim<ndim; idim++) x[idim] = 0.;
+ pdf.Eval(x, result, err);
+ Double_t threshold = err/result;
+
+ // starting interpolation over the testing grid
+ chi2 = 0.; Int_t nsample = 0;
+ x[4] = -2.+dx[4]/2.;
+ do{
+ x[3] = -2.+dx[3]/2.;
+ do{
+ x[2] = -2.+dx[2]/2.;
+ do{
+ x[1] = -2.+dx[1]/2.;
+ do{
+ x[0] = -2.+dx[0]/2.;
+ do{
+ pdf.Eval(x, result, err);
+ if(err/TMath::Abs(result) < 1.1*threshold){
+ theor = f->Eval(x[0], x[1], x[2]);
+ Double_t tmp = result - theor; tmp *= tmp;
+ chi2 += tmp/(kCHI2 ? err*err : theor);
+ nsample++;
+ }
+ } while((x[0] += dx[0]) < 2.);
+ } while((x[1] += dx[1]) < 2.);
+ } while((x[2] += dx[2]) < 2.);
+ } while((x[3] += dx[3]) < 2.);
+ } while((x[4] += dx[4]) < 2.);
+ chi2 /= float(nsample);
+ printf("\tinterpolation quality (chi2) %f\n", chi2);
+
+ return kCHI2 ? chi2 : TMath::Sqrt(chi2);
+}
+
+
+//___________________________________________________________
+void build(const int ndim, const int npoints)
+{
+ printf("Building DB ...\n");
+ Float_t data[ndim];
+ TFile::Open(Form("%dD_Gauss.root", ndim), "RECREATE");
+ TTree *t = new TTree("db", Form("%dD data base for kD statistics", ndim));
+ for(int idim=0; idim<ndim; idim++) t->Branch(Form("x%d", idim), &data[idim], Form("x%d/F", idim));
+
+ for (Int_t ip=0; ip<npoints; ip++){
+ for (Int_t id=0; id<ndim; id++) data[id]= gRandom->Gaus();
+ t->Fill();
+ }
+
+ t->Write();
+ gFile->Close();
+ delete gFile;
+}
+
+//______________________________________________________________
+Float_t Mem()
+{
+ // size in KB
+ static struct mallinfo debug;
+ debug = mallinfo();
+ //printf("arena[%d] usmblks[%d] uordblks[%d] fordblks[%d]\n", debug.arena, debug.usmblks, debug.uordblks, debug.fordblks);
+ return debug.uordblks/1024.;
+}
+
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff