STAT/Macros/interpolTest.C

   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 }