update testing macros
authorabercuci <abercuci@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 25 Oct 2007 13:52:33 +0000 (13:52 +0000)
committerabercuci <abercuci@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 25 Oct 2007 13:52:33 +0000 (13:52 +0000)
STAT/Macros/ioInterpolator.C [deleted file]
STAT/Macros/testInterpolator.C [deleted file]
STAT/Macros/testKNN.C [deleted file]
STAT/Macros/testkdtreeIF.C [deleted file]

diff --git a/STAT/Macros/ioInterpolator.C b/STAT/Macros/ioInterpolator.C
deleted file mode 100644 (file)
index 3013645..0000000
+++ /dev/null
@@ -1,60 +0,0 @@
-const int ndim = 2;
-
-//___________________________________________________________
-void writeInterpolator()
-{
-       if(!gSystem->FindFile(".", "6Ddb.root")) build(6, 1000000);
-       TFile::Open("6Ddb.root");
-       TKDInterpolator in(db, "x0:x1", "", 400, 100000);
-       in.SetIntInterpolation();
-       in.SetSetStore();
-
-       Float_t c[ndim], val, v_err;
-       Double_t p[ndim], res, r_err;
-       TGraph2DErrors *g= new TGraph2DErrors(in.GetNTerminalNodes());g->SetMarkerStyle(7);
-       for(int inode=0; inode<in.GetNTerminalNodes(); inode++){
-               in.GetCOGPoint(inode, c, val, v_err);
-               for(int idim=0; idim<ndim; idim++){
-                       p[idim] = (Double_t)c[idim];
-                       //printf("%f ", p[idim]);
-               }
-               //printf("\n");
-               in.Eval(p, res, r_err);
-               g->SetPoint(inode, p[0], p[1], res);
-               g->SetPointError(inode, 0., 0., r_err);
-       }
-       g->Draw("ap");
-
-       TFile *fi = TFile::Open(Form("%dD_interpolator.root", ndim), "RECREATE");
-       in.Write(Form("%dDgauss", ndim));
-       fi->Close();
-       delete fi;
-}
-
-//___________________________________________________________
-void readInterpolator()
-{
-       TFile::Open(Form("%dD_interpolator.root", ndim));
-       TKDInterpolator *in = (TKDInterpolator*)gFile->Get(Form("%dDgauss", ndim));
-       in->Dump();
-       //in->GetStatus();
-}
-
-//___________________________________________________________
-void build(const int ndim, const int npoints)
-{
-       printf("building db ...\n");
-       Float_t data[ndim];
-       TFile *f = new TFile(Form("%dDdb.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();
-       f->Close();
-       delete f;
-}
diff --git a/STAT/Macros/testInterpolator.C b/STAT/Macros/testInterpolator.C
deleted file mode 100644 (file)
index 8de2cd7..0000000
+++ /dev/null
@@ -1,154 +0,0 @@
-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;
-}
-
diff --git a/STAT/Macros/testKNN.C b/STAT/Macros/testKNN.C
deleted file mode 100644 (file)
index 2e131b5..0000000
+++ /dev/null
@@ -1,148 +0,0 @@
-Float_t p[]={1.4, -.6}; //}{1.7, -.4};
-//______________________________________________________________
-void testKNN(const Float_t *p, const int kNN=20)
-{
-// Draw "kNN" nearest neighbours of point "p". The distance (in the L1
-// metric) is encoded in the color code.
-// To build the data refere to function build().
-
-       TFile::Open("2D_Gauss.root");
-       TKDInterpolator in(db, "x0:x1", "x0>-1.5&&x0<2.&&x1>-2.&&x1<2.", 300);
-       in.DrawNode(in.FindNode(p)-in.GetNNodes());
-       
-       TMarker *mp = new TMarker(p[0], p[1], 3);
-       mp->SetMarkerColor(4);
-       mp->Draw();
-               
-       Int_t *index, color;
-       Float_t d, d0, pknn[2];
-       in.FindNearestNeighbors(p, kNN, index, d0);
-       TMarker *ma = new TMarker[kNN];
-       for(int ip=0; ip<kNN; ip++){
-               in.GetDataPoint(index[ip], pknn);
-               d = TMath::Abs(p[0]-pknn[0]) + TMath::Abs(p[1]-pknn[1]);
-               ma[ip].SetMarkerStyle(4);
-               color = 101 - Int_t(50. * d/d0);
-               ma[ip].SetMarkerColor(color);
-               ma[ip].DrawMarker(pknn[0], pknn[1]);
-       }
-}
-
-void performanceKNN(const int np = 10000)
-{
-       const int ntimes = 1000; // times to repeate kNN search
-       const int ndim = 5;      // no of dimensions for data
-
-       if(!gSystem->FindFile(".", Form("%dD_Gauss.root", ndim))) build(ndim);
-       TFile::Open(Form("%dD_Gauss.root", ndim));
-       TString var="x0";
-       for(int i=1; i<ndim; i++) var+=Form(":x%d", i);
-       TKDInterpolator in(db, var.Data(), "", 30, np);
-       
-       Int_t *index, fNN = 20;
-       Float_t *d, p[ndim];
-       for(int idim=0; idim<ndim; idim++) p[idim] = gRandom->Gaus();
-       TStopwatch time;
-       time.Start();
-       for(int itime=0; itime<ntimes; itime++) in.FindNearestNeighbors(p, fNN, index, d);
-       time.Stop();
-       printf("cpu = %f [mus] real = %f [mus]\n", 1.E6*time.CpuTime()/ntimes, 1.E6*time.RealTime()/ntimes);
-       for(int i=0; i<fNN; i++){
-               printf("%5d ", index[i]);
-               if(i%5==4) printf("\n");
-       }
-       
-       // draw kNN
-       TH2 *h2 = new TH2I("h2", "", 100, -2., 2., 100, -2., 2.);
-       h2->Draw();
-       TGraph *gKD = new TGraph(fNN);
-       gKD->SetMarkerStyle(24);
-       gKD->SetMarkerColor(2);
-       gKD->SetMarkerSize(.5);
-       Float_t pknn[ndim];
-       for(int ip=0; ip<fNN; ip++){
-               in.GetDataPoint(index[ip], pknn);
-               gKD->SetPoint(ip, pknn[0], pknn[1]);
-       }
-       gKD->Draw("p");
-       TMarker *mp = new TMarker(p[0], p[1], 3);
-       mp->SetMarkerColor(4);
-       mp->Draw();
-
-       
-       // STAND ALONE
-       const Int_t kNN = fNN;
-       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.;
-
-       // read data in memory
-       Int_t npoints = db->Draw("x0", "", "goff", np);
-       Float_t **x = new Float_t*[ndim];
-       Double_t *v;
-       for(int idim=0; idim<ndim; idim++){
-               x[idim] = new Float_t[npoints];
-               db->Draw(Form("x%d", idim), "", "goff", np); v = db->GetV1();
-               for(int i=0; i<npoints; i++) x[idim][i] = v[i];
-       }
-       // calculate
-       printf("stand alone calculation for %d neighbors in %d points\n", kNN, npoints);
-       time.Start(kTRUE);
-       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[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;
-               }
-       }
-       time.Stop();
-       printf("cpu = %f [s] real = %f [s]\n", time.CpuTime(), time.RealTime());
-       
-       for(int i=0; i<kNN; i++){
-               printf("%5d ", fkNN[i]);
-               if(i%5==4) 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");
-}
-
-//______________________________________________________________
-void build(const Int_t ndim = 2, const Int_t nstat = 1000000)
-{
-// Build "nstat" data points in "ndim" dimensions taken from an
-// uncorrelated 2D Gauss distribution.
-       
-
-       printf("build data ... \n");
-       Double_t pntTrain[ndim];
-       f = TFile::Open(Form("%dD_Gauss.root", ndim), "RECREATE");
-       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;
-}
\ No newline at end of file
diff --git a/STAT/Macros/testkdtreeIF.C b/STAT/Macros/testkdtreeIF.C
deleted file mode 100644 (file)
index 55627c2..0000000
+++ /dev/null
@@ -1,189 +0,0 @@
-// #include <malloc.h>
-// #include "TMatrixD.h"
-// #include "TRandom.h"
-// #include "TGraph.h"
-// #include "TStopwatch.h"
-// 
-
-
-Int_t Mem()
-{
-  // size in 1000 bytes
-  static struct mallinfo memdebug; 
-  memdebug = mallinfo();
-  return memdebug.uordblks/1000;
-}
-
-
-
-void testindexes(Int_t nloop, Int_t npoints);
-void  testkdtreeIF(Int_t npoints=1000, Int_t nloop=1000, Int_t mode = 1, Int_t bsize=9);
-void TestSizeIF(Int_t npoints=1000, Int_t nrows = 150, Int_t nsec=18, Int_t mode = 1, Int_t bsize=9);
-
-
-
-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 TestBuild(Int_t npoints, 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();
-    //data[1][i]= 0;
-    //data[0][i]= -i;
-  }
-  Float_t dataf[2];
-  TKDTreeIF *kdtree = new TKDTreeIF(npoints, 2, bsize, data);
-  //kdtree->Build();
-  for (Int_t i=0; i<npoints;i++){
-    Int_t index = kdtree->fIndPoints[i];
-    //printf("%d\t%d\t%f\n",i,index,data[0][index]);
-  }
-  for (Int_t i=0; i<kdtree->fNnodes;i++){
-    //printf("%d\t%f\n",i,kdtree->fNodes[i].fValue);
-  }
-  Float_t sumiter = 0;
-  for (Int_t i=0;i<npoints;i++){
-    dataf[0] = data[0][i];
-    dataf[1] = data[1][i];
-    Int_t index=-1;
-    Int_t iter =0;
-    kdtree->FindPoint(dataf,index, iter);
-    if (i!=index){
-      printf("%d\t%d\t%f\t%f\n",i,index,dataf[0],data[0][index]);
-    }
-    sumiter+=iter;
-  }
-  printf("Mean iter = %f\n",float(sumiter)/float(npoints));
-}
-
-
-void TestSizeIF(Int_t npoints, Int_t nrows, Int_t nsec, Int_t mode, Int_t bsize)
-{
-  //
-  // test size to build kdtree
-  //
-  Int_t before =Mem();
-  for (Int_t isec=0; isec<nsec;isec++)
-    for (Int_t irow=0;irow<nrows;irow++){
-      testkdtreeIF(npoints,0,mode,bsize);
-    }
-  Int_t after = Mem();
-  printf("Memory usage %d\n",after-before);
-}
-
-
-void  testkdtreeIF(Int_t npoints, Int_t nloop, Int_t mode, Int_t bsize)
-{
-  //
-  // test speed and functionality of kdtree
-  //
-  Float_t rangey  = 100;
-  Float_t rangez  = 100;
-  Float_t drangey = 0.1;
-  Float_t drangez = 0.1;
-//   Float_t rangey  = 20;
-//   Float_t rangez  = 250;
-//   Float_t drangey = 1;
-//   Float_t drangez = 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);
-    //     data[i+npoints]  = TMath::Nint(gRandom->Uniform(-rangez, rangez)/10.)*10.;
-    //printf("%d %f  %f\n", i, data[i], data[i+npoints]);     
-  }
-  TStopwatch timerbuild;
-  TKDTree<Int_t, Float_t> *kdtree = new TKDTree<Int_t, Float_t>(npoints,2,bsize,data);
-   kdtree->Build();
-   timerbuild.Stop();
-   timerbuild.Print();
-  TStopwatch timer;
-
-   Float_t countern=0;
-   Float_t counteriter  = 0;
-   Float_t counterfound = 0;
-
-
-   if (mode ==2){ 
-     if (nloop) timer.Start();
-     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); 
-}