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

// #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;
}
Commit	Line	Data
8a8f8274	1	// #include "TSystem.h"
	2	// #include "TFile.h"
	3	// #include "TTree.h"
	4	// #include "TProfile.h"
	5	// #include "TF1.h"
	6	// #include "TF2.h"
	7	// #include "TF3.h"
	8	// #include "TLegend.h"
	9	// #include "TRandom.h"
	10	// #include "TString.h"
	11	// #include "TGraph.h"
	12	// #include "TMarker.h"
	13	// #include "TStopwatch.h"
	14	// #include "../src/TKDPDF.h"
	15	// #include "../src/TKDInterpolator.h"
	16
	17	void interpolTest(const Int_t ndim = 1);
	18	Double_t interpolation(const int nstat, const int ndim, const int first);
	19	void build(const int ndim, const int npoints);
	20
	21	TF1 *f = 0x0;
	22
	23	//__________________________________________________________
	24	void interpolTest(const Int_t ndim)
	25	{
	26	// Test interpolator on a variable dimension (ndim<=5) uncorrelated
	27	// Gauss model. The macro displays the chi2 between interpolation and
	28	// model for various statistics of experimental data.
	29	//
	30	// Check the worker function interpolation() to fine tune the
	31	// algorithm.
	32	//
	33	// Attention:
	34	// The macro works in compiled mode !
	35
	36	// define model
	37	switch(ndim){
	38	case 1:
	39	f=new TF1("f1", "gaus(0);x;f(x)", -5., 5.);
	40	f->SetParameter(0, 1.);
	41	f->SetParameter(1, 0.);
	42	f->SetParameter(2, 1.);
	43	f->SetParameter(0, 1./f->Integral(-5., 5.));
	44	break;
	45	case 2:
	46	f=new TF2("f2", "[0]exp(-0.5((x-[1])/[2])*2)exp(-0.5((y-[3])/[4])*2)", -5., 5., -5., 5.);
	47	f->SetParameter(0, 1.);
	48	f->SetParameter(1, 0.);
	49	f->SetParameter(2, 1.);
	50	f->SetParameter(3, 0.);
	51	f->SetParameter(4, 1.);
	52	f->SetParameter(0, 1./f->Integral(-5., 5., -5., 5.));
	53	break;
	54	case 3:
	55	case 4:
	56	case 5:
	57	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.);
	58	f->SetParameter(0, 1.);
	59	f->SetParameter(1, 0.);
	60	f->SetParameter(2, 1.);
	61	f->SetParameter(3, 0.);
	62	f->SetParameter(4, 1.);
	63	f->SetParameter(5, 0.);
	64	f->SetParameter(6, 1.);
65	f->SetParameter(0, 1./f->Integral(-5., 5., -5., 5., -5., 5.));
66	if(ndim>3) printf("Warning : chi2 results are unreliable !\n");
67	break;
68	default:
69	printf("%dD not supported in this macro.\n", ndim);
70	return;
71	}
72
73	Double_t chi2;
74	const Int_t nchi2 = 18;
75	const Int_t nfirst = 1;
76	//TProfile *hp = new TProfile("hp", "", 100, 1.E4, 1.E6);
77	//hp->SetMarkerStyle(24);
78	TGraph *gChi2 = new TGraph(nchi2);
79	gChi2->SetMarkerStyle(24);
80	Int_t ip = 0, nstat, first;
81	for(int ifirst=0; ifirst<nfirst; ifirst++){
82	first = 0;//Int_t(gRandom->Uniform(1.E4));
83	for(int i=2; i<10; i++){
84	nstat = 10000*i;
85	chi2 = interpolation(nstat, ndim, first);
86	gChi2->SetPoint(ip++, TMath::Log10(float(nstat)), chi2);
87	//hp->Fill(float(nstat), chi2);
88	}
89	for(int i=1; i<10; i++){
90	nstat = 100000*i;
91	chi2 = interpolation(nstat, ndim, first);
92	gChi2->SetPoint(ip++, TMath::Log10(float(nstat)), chi2);
93	//hp->Fill(float(nstat), chi2);
94	}
95	gChi2->Draw("apl");
96	}
97	//hp->Draw("pl");
98	}
99
100	void test()
101	{
102	TClonesArray ar("TKDNodeInfo", 10);
103	}
104
105	//__________________________________________________________
106	Double_t interpolation(const int nstat, const int ndim, const int first)
107	{
108	// Steer interpolation of a nD (n<=5) uncorrelated Gauss distribution.
109	// The user is suppose to give the experimental statistics (nstat) the
110	// dimension of the experimental point (ndim). The entry from which the
111	// data base is being read is steered by "first".
112
113	const int bs = 400;
114	const int npoints = 100;
115	// switch on chi2 calculation between interpolation and model.
116	// Default calculates distance.
117	const Bool_t kCHI2 = kFALSE;
118	const Bool_t kINT = kFALSE; // switch on integral interpolator
119
120	// open data base
121	TString fn("5D_Gauss.root");
122	if(!gSystem->FindFile(".", fn)) build(5, 1000000);
123	TFile::Open("5D_Gauss.root");
124	TTree t = (TTree)gFile->Get("db");
125
126	// build interpolator
127	TString var = "x0"; for(int idim=1; idim<ndim; idim++) var += Form(":x%d", idim);
128	TKDPDF pdf(t, var.Data(), "", bs, nstat, first);
129	printf("\n\nFINISH BUILDING PDF\n\n");
130
131	Double_t fx;
132	Int_t inode = 0;
133	TKDNodeInfo *node = 0x0;
134	TKDInterpolator in(ndim, pdf.GetNTNodes());
135	/* while(node = pdf.GetNodeInfo(inode)){
136	if(node->fVal[0] > 0.){
137	fx = node->fVal[0];
138	node->fVal[0] = TMath::Log(fx);
139	node->fVal[1] = node->fVal[1]/fx;
140	}
141	in.SetNode(inode, *node);
142	inode++;
143	}*/
144	//pdf.SetInterpolationMethod(kINT);
145	pdf.SetAlpha(2.);
146	//in.SetStore();
147
148
149	Double_t x[2], r, e, chi2;
150	TH2 *h2 = new TH2F("h2", "", 50, -4., 4., 50, -4., 4.);
151	TAxis ax = h2->GetXaxis(), ay = h2->GetYaxis();
152	for(int ix=1; ix<=ax->GetNbins(); ix++){
153	x[0] = ax->GetBinCenter(ix);
154	for(int iy=1; iy<=ay->GetNbins(); iy++){
155	x[1] = ay->GetBinCenter(iy);
156
157	chi2 = pdf.Eval(x, r, e);
158	printf("x[%2d] x[%2d] r[%f] e[%f] chi2[%f]\n", ix, iy, r, e, chi2);
159	h2->SetBinContent(ix, iy, r/TMath::Exp(r)/);
160	}
161	}
162	h2->Draw("lego2");
163	return;
164
165	// define testing grid
166	Double_t x[5], r, e;
167	Double_t dx[5] = {4., 4., 4., 4., 4.};
168	Double_t chi2, theor, result, err;
169	for(int idim=0; idim<ndim; idim++){
170	dx[idim] = (ndim<<2)/float(npoints);
171	x[idim] = -2.+dx[idim]/2.;
172	}
173
174
175	// setup the integral interpolator and do memory test
176	Float_t c[5], v, ve;
177	for(int ip=0; ip<in.GetNTNodes(); ip++){
178	in.GetCOGPoint(ip, c, v, ve);
179	for(int idim=0; idim<ndim; idim++) x[idim] = (Double_t)c[idim];
180	in.Eval(x, result, err);
181	}
182	printf("\nInterpolating %d data points in %dD ...\n", nstat, ndim);
183	printf("\tbucket size %d\n", bs);
184	printf("\tstarting interpolation ...\n");
185	return;
186
187
188	// evaluate relative error
189	for(int idim=0; idim<ndim; idim++) x[idim] = 0.;
190	in.Eval(x, result, err);
191	Double_t threshold = err/result;
192
193	// starting interpolation over the testing grid
194	chi2 = 0.; Int_t nsample = 0;
195	x[4] = -2.+dx[4]/2.;
196	do{
197	x[3] = -2.+dx[3]/2.;
198	do{
199	x[2] = -2.+dx[2]/2.;
200	do{
201	x[1] = -2.+dx[1]/2.;
202	do{
203	x[0] = -2.+dx[0]/2.;
204	do{
205	in.Eval(x, result, err);
206	if(err/TMath::Abs(result) < 1.1*threshold){
207	theor = f->Eval(x[0], x[1], x[2]);
208	Double_t tmp = result - theor; tmp *= tmp;
209	chi2 += tmp/(kCHI2 ? err*err : theor);
210	nsample++;
211	}
212	} while((x[0] += dx[0]) < 2.);
213	} while((x[1] += dx[1]) < 2.);
214	} while((x[2] += dx[2]) < 2.);
215	} while((x[3] += dx[3]) < 2.);
216	} while((x[4] += dx[4]) < 2.);
217	chi2 /= float(nsample);
218	printf("\tinterpolation quality (chi2) %f\n", chi2);
219
220	return kCHI2 ? chi2 : TMath::Sqrt(chi2);
221	}
222
223
224	//___________________________________________________________
225	void build(const int ndim, const int npoints)
226	{
227	printf("Building DB ...\n");
228	Float_t data[ndim];
229	TFile::Open(Form("%dD_Gauss.root", ndim), "RECREATE");
230	TTree *t = new TTree("db", Form("%dD data base for kD statistics", ndim));
231	for(int idim=0; idim<ndim; idim++) t->Branch(Form("x%d", idim), &data[idim], Form("x%d/F", idim));
232
233	for (Int_t ip=0; ip<npoints; ip++){
234	for (Int_t id=0; id<ndim; id++) data[id]= gRandom->Gaus();
235	t->Fill();
236	}
237
238	t->Write();
239	gFile->Close();
240	delete gFile;
241	}