ROOT logo

From $ROOTSYS/tutorials/math/multidimSampling.C

//+ example of sampling a multi-dim distribution using the DistSampler class
//NOTE: This tutorial must be run with ACLIC
//Author L. Moneta Dec 2010

// function (a 4d gaussian)
#include "TMath.h"
#include "TF2.h"
#include "TStopwatch.h"
#include "Math/DistSampler.h"
#include "Math/DistSamplerOptions.h"
#include "Math/MinimizerOptions.h"
#include "Math/Factory.h"

#include "TKDTreeBinning.h"

#include "TTree.h"
#include "TFile.h"
#include "TMatrixDSym.h"
#include "TVectorD.h"
#include "TCanvas.h"
#include <cmath>

// Gauss ND function
// make a class in order to avoid constructing the 
// matrices for every call
// This however requires that  the code must be compiled with ACLIC

bool debug = false;
struct GausND { 
   
   TVectorD X; 
   TVectorD Mu; 
   TMatrixDSym CovMat; 
   
   GausND( int dim ) : 
      X(TVectorD(dim)), 
      Mu(TVectorD(dim)), 
      CovMat(TMatrixDSym(dim) )
   {}
   
   double operator() (double *x, double *p) { 
      // 4 parameters 
      int dim = X.GetNrows();
      int k = 0;
      for (int i = 0; i<dim; ++i) { X[i] = x[i] - p[k]; k++; } 
      for (int i = 0; i<dim; ++i) { 
         CovMat(i,i) = p[k]*p[k];
         k++;
      }
      for (int i = 0; i<dim; ++i) { 
         for (int j = i+1; j<dim; ++j) { 
            // p now are the correlations N(N-1)/2
               CovMat(i,j) = p[k]*sqrt(CovMat(i,i)*CovMat(j,j));            
               CovMat(j,i) = CovMat(i,j);
               k++;
         }
      }
      if (debug) { 
         X.Print();
         CovMat.Print();
      }

      double det = CovMat.Determinant(); 
      if (det <= 0) {
         Fatal("GausND","Determinant is <= 0 det = %f",det);
         CovMat.Print();
         return 0;  
      }
      double norm = std::pow( 2. * TMath::Pi(), dim/2) * sqrt(det);  
      // compute the gaussians
      CovMat.Invert();
      double fval  = std::exp( - 0.5 * CovMat.Similarity(X) )/ norm; 

      if (debug) { 
         std::cout << "det  " << det << std::endl; 
         std::cout << "norm " << norm << std::endl; 
         std::cout << "fval " << fval << std::endl; 
      }

      return fval;
   }
};

using namespace ROOT::Math; 

void multidimSampling() { 

#ifdef __CINT__
   std::cout << "DO NOT RUN WITH CINT:" << std::endl;
   std::cout << "we are using a custom function which requires" << std::endl;
   std::cout << "that this tutorial must be compiled with ACLIC" << std::endl;
  return;
#endif

   const int N = 100000;
   //const int NBin = 1000;
   const int DIM = 4; 

   double xmin[] = {-10,-10,-10, -10};
   double xmax[] = { 10, 10, 10,  10};
   double par0[] = { 1., -1., 2, 0, // the gaussian mu
                     1, 2, 1, 3, // the sigma
                     0.5,0.,0.,0.,0.,0.8 };  // the correlation
   
   const int NPAR = DIM + DIM*(DIM+1)/2; // 14 in the 4d case
   // generate the sample 
   GausND gaus4d(4);
   TF1 * f = new TF1("functionND",gaus4d,0,1,14);
   f->SetParameters(par0);

   double x0[] = {0,0,0,0};
   // for debugging
   if (debug) f->EvalPar(x0,0);
   debug = false;

   TString name; 
   for (int i = 0; i < NPAR; ++i )  {
      if (i < DIM) f->SetParName(i, name.Format("mu_%d",i+1) );
      else if (i < 2*DIM) f->SetParName(i, name.Format("sig_%d",i-DIM+1) );
      else if (i < 2*DIM) f->SetParName(i, name.Format("sig_%d",i-2*DIM+1) );
   }

   //ROOT::Math::DistSamplerOptions::SetDefaultSampler("Foam");
   DistSampler * sampler = Factory::CreateDistSampler();
   if (sampler == 0) { 
      Info("multidimSampling","Default sampler %s is not available try with Foam ",
           ROOT::Math::DistSamplerOptions::DefaultSampler().c_str() );           
      ROOT::Math::DistSamplerOptions::SetDefaultSampler("Foam");
   }
   sampler = Factory::CreateDistSampler();
   if (sampler == 0) { 
      Error("multidimSampling","Foam sampler is not available - exit ");
      return;
   }
            

   sampler->SetFunction(*f,DIM);
   sampler->SetRange(xmin,xmax);
   bool ret = sampler->Init();

   std::vector<double> data1(DIM*N); 
   double v[DIM];
   TStopwatch w; 

   if (!ret) { 
      Error("Sampler::Init","Error initializing unuran sampler");
      return; 
   }
   // generate the data 
   w.Start(); 
   for (int i = 0; i < N; ++i) { 
      sampler->Sample(v);
      for (int j = 0; j < DIM; ++j) 
         data1[N*j + i]     = v[j];
   }
   w.Stop(); 
   w.Print(); 

   // fill tree with data 
   TFile * file = new TFile("multiDimSampling.root","RECREATE");
   double x[DIM];
   TTree * t1 = new TTree("t1","Tree from Unuran"); 
   t1->Branch("x",x,"x[4]/D");
   for (int i = 0; i < N; ++i) { 
      for (int j = 0; j < DIM; ++j) {
         x[j] = data1[i+N*j];
      }
      t1->Fill();
   }


   // try to fit the 2d data; 
   // GausND gaus2(2); 
   // TF2 * f2 = new TF2("f2",gaus2,-3,3,-3,3,5,"GausND");
   
   // f2->SetParameters(0,0,1,1,0);
   // f2->SetParLimits(4,-1,1);
   // t1->UnbinnedFit("f2","x[0]:x[1]");

   // plot the data
   
   t1->Draw("x[0]:x[1]:x[2]:x[3]","","candle");
   TCanvas * c2 = new TCanvas();
   c2->Divide(3,2);
   int ic=1;
   c2->cd(ic++);
   t1->Draw("x[0]:x[1]");
   c2->cd(ic++);
   t1->Draw("x[0]:x[2]");
   c2->cd(ic++);
   t1->Draw("x[0]:x[3]");
   c2->cd(ic++);
   t1->Draw("x[1]:x[2]");
   c2->cd(ic++);
   t1->Draw("x[1]:x[3]");
   c2->cd(ic++);
   t1->Draw("x[2]:x[3]");
   

   t1->Write();
   file->Close();


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