ROOT logo
#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.;
}

 pdfTest.C:1
 pdfTest.C:2
 pdfTest.C:3
 pdfTest.C:4
 pdfTest.C:5
 pdfTest.C:6
 pdfTest.C:7
 pdfTest.C:8
 pdfTest.C:9
 pdfTest.C:10
 pdfTest.C:11
 pdfTest.C:12
 pdfTest.C:13
 pdfTest.C:14
 pdfTest.C:15
 pdfTest.C:16
 pdfTest.C:17
 pdfTest.C:18
 pdfTest.C:19
 pdfTest.C:20
 pdfTest.C:21
 pdfTest.C:22
 pdfTest.C:23
 pdfTest.C:24
 pdfTest.C:25
 pdfTest.C:26
 pdfTest.C:27
 pdfTest.C:28
 pdfTest.C:29
 pdfTest.C:30
 pdfTest.C:31
 pdfTest.C:32
 pdfTest.C:33
 pdfTest.C:34
 pdfTest.C:35
 pdfTest.C:36
 pdfTest.C:37
 pdfTest.C:38
 pdfTest.C:39
 pdfTest.C:40
 pdfTest.C:41
 pdfTest.C:42
 pdfTest.C:43
 pdfTest.C:44
 pdfTest.C:45
 pdfTest.C:46
 pdfTest.C:47
 pdfTest.C:48
 pdfTest.C:49
 pdfTest.C:50
 pdfTest.C:51
 pdfTest.C:52
 pdfTest.C:53
 pdfTest.C:54
 pdfTest.C:55
 pdfTest.C:56
 pdfTest.C:57
 pdfTest.C:58
 pdfTest.C:59
 pdfTest.C:60
 pdfTest.C:61
 pdfTest.C:62
 pdfTest.C:63
 pdfTest.C:64
 pdfTest.C:65
 pdfTest.C:66
 pdfTest.C:67
 pdfTest.C:68
 pdfTest.C:69
 pdfTest.C:70
 pdfTest.C:71
 pdfTest.C:72
 pdfTest.C:73
 pdfTest.C:74
 pdfTest.C:75
 pdfTest.C:76
 pdfTest.C:77
 pdfTest.C:78
 pdfTest.C:79
 pdfTest.C:80
 pdfTest.C:81
 pdfTest.C:82
 pdfTest.C:83
 pdfTest.C:84
 pdfTest.C:85
 pdfTest.C:86
 pdfTest.C:87
 pdfTest.C:88
 pdfTest.C:89
 pdfTest.C:90
 pdfTest.C:91
 pdfTest.C:92
 pdfTest.C:93
 pdfTest.C:94
 pdfTest.C:95
 pdfTest.C:96
 pdfTest.C:97
 pdfTest.C:98
 pdfTest.C:99
 pdfTest.C:100
 pdfTest.C:101
 pdfTest.C:102
 pdfTest.C:103
 pdfTest.C:104
 pdfTest.C:105
 pdfTest.C:106
 pdfTest.C:107
 pdfTest.C:108
 pdfTest.C:109
 pdfTest.C:110
 pdfTest.C:111
 pdfTest.C:112
 pdfTest.C:113
 pdfTest.C:114
 pdfTest.C:115
 pdfTest.C:116
 pdfTest.C:117
 pdfTest.C:118
 pdfTest.C:119
 pdfTest.C:120
 pdfTest.C:121
 pdfTest.C:122
 pdfTest.C:123
 pdfTest.C:124
 pdfTest.C:125
 pdfTest.C:126
 pdfTest.C:127
 pdfTest.C:128
 pdfTest.C:129
 pdfTest.C:130
 pdfTest.C:131
 pdfTest.C:132
 pdfTest.C:133
 pdfTest.C:134
 pdfTest.C:135
 pdfTest.C:136
 pdfTest.C:137
 pdfTest.C:138
 pdfTest.C:139
 pdfTest.C:140
 pdfTest.C:141
 pdfTest.C:142
 pdfTest.C:143
 pdfTest.C:144
 pdfTest.C:145
 pdfTest.C:146
 pdfTest.C:147
 pdfTest.C:148
 pdfTest.C:149
 pdfTest.C:150
 pdfTest.C:151
 pdfTest.C:152
 pdfTest.C:153
 pdfTest.C:154
 pdfTest.C:155
 pdfTest.C:156
 pdfTest.C:157
 pdfTest.C:158
 pdfTest.C:159
 pdfTest.C:160
 pdfTest.C:161
 pdfTest.C:162
 pdfTest.C:163
 pdfTest.C:164
 pdfTest.C:165
 pdfTest.C:166
 pdfTest.C:167
 pdfTest.C:168
 pdfTest.C:169
 pdfTest.C:170
 pdfTest.C:171
 pdfTest.C:172
 pdfTest.C:173
 pdfTest.C:174
 pdfTest.C:175
 pdfTest.C:176
 pdfTest.C:177
 pdfTest.C:178
 pdfTest.C:179
 pdfTest.C:180
 pdfTest.C:181
 pdfTest.C:182
 pdfTest.C:183
 pdfTest.C:184
 pdfTest.C:185
 pdfTest.C:186
 pdfTest.C:187
 pdfTest.C:188
 pdfTest.C:189
 pdfTest.C:190
 pdfTest.C:191
 pdfTest.C:192
 pdfTest.C:193
 pdfTest.C:194
 pdfTest.C:195
 pdfTest.C:196
 pdfTest.C:197
 pdfTest.C:198
 pdfTest.C:199
 pdfTest.C:200
 pdfTest.C:201
 pdfTest.C:202
 pdfTest.C:203
 pdfTest.C:204
 pdfTest.C:205
 pdfTest.C:206
 pdfTest.C:207
 pdfTest.C:208
 pdfTest.C:209
 pdfTest.C:210
 pdfTest.C:211
 pdfTest.C:212
 pdfTest.C:213
 pdfTest.C:214
 pdfTest.C:215
 pdfTest.C:216
 pdfTest.C:217
 pdfTest.C:218
 pdfTest.C:219
 pdfTest.C:220