+++ /dev/null
-const Int_t ndim = 2;
-Double_t testInterpolator(const Int_t nstat = 100000)
-{
-// Macro for testing the TKDInterpolator.
-//
-// The function which it is interpolated is an uncorrelated landau
-// distribution in "ndim" dimensions. The function returns the chi2 of
-// the interpolation.
-//
-// Parameters
-// nstat - number of points to be used for training
-// kBuild - on/off generate data
-// kTransform - on/off outliers compresion
-// kPCA - on/off pricipal component analysis
-
- const Bool_t kBuild = 1;
- const Bool_t kTransform = 1;
- const Bool_t kPCA = 1;
- gStyle->SetPalette(1);
-
- Double_t pntTrain[ndim], pntTest[ndim], pntRotate[ndim];
- Double_t pdf;
- TFile *f = 0x0, *fEval = 0x0;
- TTree *t = 0x0, *tEval = 0x0;
-
-
- // build data
- if(kBuild){
- printf("build data ... \n");
- f = TFile::Open(Form("%dD_LL.root", ndim), "RECREATE");
- t = new TTree("db", "Log-Log database");
- for(int idim=0; idim<ndim; idim++) t->Branch(Form("x%d", idim), &pntTrain[idim], Form("x%d/D", idim));
-
- fEval = TFile::Open(Form("%dD_Eval.root", ndim), "RECREATE");
- tEval = new TTree("db", "Evaluation database");
- for(int idim=0; idim<ndim; idim++) tEval->Branch(Form("x%d", idim), &pntTest[idim], Form("x%d/D", idim));
-
- for(int istat=0; istat<nstat; istat++){
- for(int idim=0; idim<ndim; idim++) pntTrain[idim] = gRandom->Landau(5.);
- if(!(istat%3)){ // one third of the statistics is for testing
- memcpy(pntTest, pntTrain, ndim*sizeof(Double_t));
- tEval->Fill();
- continue;
- }
- if(kTransform)
- for(int idim=0; idim<ndim; idim++)
- if(pntTrain[idim] > 0.) pntTrain[idim] = TMath::Log(pntTrain[idim]);
- else pntTrain[idim] = 0.;
- t->Fill();
- }
- f->cd();
- t->Write();
- f->Flush();
-
- fEval->cd();
- tEval->Write();
- fEval->Flush();
- } else {// link data
- printf("link data ... \n");
- f = TFile::Open(Form("%dD_LL.root", ndim));
- t = (TTree*)f->Get("db");
- for(int idim=0; idim<ndim; idim++) t->SetBranchAddress(Form("x%d", idim), &pntTrain[idim]);
-
- fEval = TFile::Open(Form("%dD_Eval.root", ndim));
- tEval = (TTree*)fEval->Get("db");
- for(int idim=0; idim<ndim; idim++) tEval->SetBranchAddress(Form("x%d", idim), &pntTest[idim]);
- }
-
-
- // do principal component analysis (PCA)
- TPrincipal princ(ndim, "N");
- if(kPCA && kBuild){
- printf("do principal component analysis (PCA) ... \n");
- f->cd();
- TTree *tt = new TTree("db1", "PCA database");
- for(int idim=0; idim<ndim; idim++) tt->Branch(Form("x%d", idim), &pntRotate[idim], Form("x%d/D", idim));
- for(int ientry=0; ientry<t->GetEntries(); ientry++){
- t->GetEntry(ientry);
- princ.AddRow(pntTrain);
- }
- princ.MakePrincipals();
- for(int ientry=0; ientry<t->GetEntries(); ientry++){
- t->GetEntry(ientry);
- princ.X2P(pntTrain, pntRotate);
- tt->Fill();
- }
- tt->Write();
- f->Flush();
- for(int idim=0; idim<ndim; idim++) tt->SetBranchAddress(Form("x%d", idim), &pntTrain[idim]);
- t = tt;
- }
- gROOT->cd();
-
- // do interpolation
- printf("do interpolation ... \n");
- Double_t pdf, pdf_estimate, chi2;
- TString vl = "x0";
- for(int idim=1; idim<ndim; idim++) vl+=Form(":x%d", idim);
- TKDInterpolator in(t, vl.Data(), "", 200.);
- chi2 = 0.;
-/* for(int ip=0; ip<tEval->GetEntries(); ip++){
- tEval->GetEntry(ip);
- printf("\nEval %d\n", ip);*/
- TH1 *h1 = new TH2F("h1", "", 50, 0., 100., 50, 0., 100.);
- TH1 *h2 = new TH2F("h2", "", 50, 0., 100., 50, 0., 100.);
- TAxis *ax = h2->GetXaxis(), *ay = h2->GetYaxis();
- for(int ix=2; ix<ax->GetNbins(); ix++){
- pntTest[0] = ax->GetBinCenter(ix);
- for(int iy=2; iy<ay->GetNbins(); iy++){
- pntTest[1] = ay->GetBinCenter(iy);
- memcpy(pntTrain, pntTest, ndim*sizeof(Double_t));
-
- if(kTransform)
- for(int idim=0; idim<ndim; idim++)
- if(pntTrain[idim] > 0.) pntTrain[idim] = TMath::Log(pntTrain[idim]);
- else pntTrain[idim] = 0.;
-
- if(kPCA){
- princ.X2P(pntTrain, pntRotate);
- memcpy(pntTrain, pntRotate, ndim*sizeof(Double_t));
- }
-
- pdf_estimate = in.Eval(pntTrain, 30);
- // calculate chi2
- if(kTransform)
- for(int idim=0; idim<ndim; idim++)
- if(pntTest[idim] > 0.) pdf_estimate /= pntTest[idim];
- else continue;
-
- h1->SetBinContent(ix, iy, pdf_estimate);
-
- pdf = 1.; for(int idim=0; idim<ndim; idim++) pdf *= TMath::Landau(pntTest[idim], 5.);
- h2->SetBinContent(ix, iy, pdf);
- pdf_estimate -= pdf;
- chi2 += pdf_estimate*pdf_estimate/pdf;
- }}
- f->Close(); delete f;
- fEval->Close(); delete fEval;
-
- // results presentation
- printf("chi2 = %f\n", chi2);
- TCanvas *c = 0x0;
- if(!(c = (TCanvas*)gROOT->FindObject("c"))){
- c = new TCanvas("c", "", 10, 10, 900, 500);
- c->Divide(2, 1);
- }
- c->cd(1);
- h1->Draw("lego2"); h1->GetZaxis()->SetRangeUser(1.e-9, 5.e-2); gPad->SetLogz(); gPad->Modified(); gPad->Update();
-
- c->cd(2);
- h2->Draw("lego2"); h2->GetZaxis()->SetRangeUser(1.e-9, 5.e-2); gPad->SetLogz(); gPad->Modified(); gPad->Update();
- return chi2;
-}
-
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff