Logo ROOT   6.14/05
Reference Guide
multidimSampling.C
Go to the documentation of this file.
1 /// \file
2 /// \ingroup tutorial_math
3 /// \notebook
4 /// Example of sampling a multi-dim distribution using the DistSampler class
5 /// NOTE: This tutorial must be run with ACLIC
6 ///
7 /// \macro_image
8 /// \macro_code
9 ///
10 /// \author Lorenzo Moneta
11 
12 // function (a 4d gaussian)
13 #include "TMath.h"
14 #include "TF2.h"
15 #include "TStopwatch.h"
16 #include "Math/DistSampler.h"
18 #include "Math/MinimizerOptions.h"
19 #include "Math/Factory.h"
20 
21 #include "TKDTreeBinning.h"
22 
23 #include "TTree.h"
24 #include "TFile.h"
25 #include "TMatrixDSym.h"
26 #include "TVectorD.h"
27 #include "TCanvas.h"
28 #include <cmath>
29 
30 // Gauss ND function
31 // make a class in order to avoid constructing the
32 // matrices for every call
33 // This however requires that the code must be compiled with ACLIC
34 
35 bool debug = false;
36 
37 // Define the GausND strcture
38 struct GausND {
39 
40  TVectorD X;
41  TVectorD Mu;
42  TMatrixDSym CovMat;
43 
44  GausND( int dim ) :
45  X(TVectorD(dim)),
46  Mu(TVectorD(dim)),
47  CovMat(TMatrixDSym(dim) )
48  {}
49 
50  double operator() (double *x, double *p) {
51  // 4 parameters
52  int dim = X.GetNrows();
53  int k = 0;
54  for (int i = 0; i<dim; ++i) { X[i] = x[i] - p[k]; k++; }
55  for (int i = 0; i<dim; ++i) {
56  CovMat(i,i) = p[k]*p[k];
57  k++;
58  }
59  for (int i = 0; i<dim; ++i) {
60  for (int j = i+1; j<dim; ++j) {
61  // p now are the correlations N(N-1)/2
62  CovMat(i,j) = p[k]*sqrt(CovMat(i,i)*CovMat(j,j));
63  CovMat(j,i) = CovMat(i,j);
64  k++;
65  }
66  }
67  if (debug) {
68  X.Print();
69  CovMat.Print();
70  }
71 
72  double det = CovMat.Determinant();
73  if (det <= 0) {
74  Fatal("GausND","Determinant is <= 0 det = %f",det);
75  CovMat.Print();
76  return 0;
77  }
78  double norm = std::pow( 2. * TMath::Pi(), dim/2) * sqrt(det);
79  // compute the gaussians
80  CovMat.Invert();
81  double fval = std::exp( - 0.5 * CovMat.Similarity(X) )/ norm;
82 
83  if (debug) {
84  std::cout << "det " << det << std::endl;
85  std::cout << "norm " << norm << std::endl;
86  std::cout << "fval " << fval << std::endl;
87  }
88 
89  return fval;
90  }
91 };
92 
93 // Use the Math namespace
94 using namespace ROOT::Math;
95 
96 void multidimSampling() {
97 
98 
99  const int N = 10000;
100  /*const int NBin = 1000;*/
101  const int DIM = 4;
102 
103  double xmin[] = {-10,-10,-10, -10};
104  double xmax[] = { 10, 10, 10, 10};
105  double par0[] = { 1., -1., 2, 0, // the gaussian mu
106  1, 2, 1, 3, // the sigma
107  0.5,0.,0.,0.,0.,0.8 }; // the correlation
108 
109  const int NPAR = DIM + DIM*(DIM+1)/2; // 14 in the 4d case
110  // generate the sample
111  GausND gaus4d(4);
112  TF1 * f = new TF1("functionND",gaus4d,0,1,14);
113  f->SetParameters(par0);
114 
115  double x0[] = {0,0,0,0};
116  // for debugging
117  if (debug) f->EvalPar(x0,0);
118  debug = false;
119 
120  TString name;
121  for (int i = 0; i < NPAR; ++i ) {
122  if (i < DIM) f->SetParName(i, name.Format("mu_%d",i+1) );
123  else if (i < 2*DIM) f->SetParName(i, name.Format("sig_%d",i-DIM+1) );
124  else if (i < 2*DIM) f->SetParName(i, name.Format("sig_%d",i-2*DIM+1) );
125  }
126 
127  /*ROOT::Math::DistSamplerOptions::SetDefaultSampler("Foam");*/
128  DistSampler * sampler = Factory::CreateDistSampler();
129  if (sampler == 0) {
130  Info("multidimSampling","Default sampler %s is not available try with Foam ",
133  }
134  sampler = Factory::CreateDistSampler();
135  if (sampler == 0) {
136  Error("multidimSampling","Foam sampler is not available - exit ");
137  return;
138  }
139 
140  sampler->SetFunction(*f,DIM);
141  sampler->SetRange(xmin,xmax);
142  bool ret = sampler->Init();
143 
144  std::vector<double> data1(DIM*N);
145  double v[DIM];
146  TStopwatch w;
147 
148  if (!ret) {
149  Error("Sampler::Init","Error initializing unuran sampler");
150  return;
151  }
152 
153  // generate the data
154  w.Start();
155  for (int i = 0; i < N; ++i) {
156  sampler->Sample(v);
157  for (int j = 0; j < DIM; ++j)
158  data1[N*j + i] = v[j];
159  }
160  w.Stop();
161  w.Print();
162 
163  // fill tree with data
164  TFile * file = new TFile("multiDimSampling.root","RECREATE");
165  double x[DIM];
166  TTree * t1 = new TTree("t1","Tree from Unuran");
167  t1->Branch("x",x,"x[4]/D");
168  for (int i = 0; i < N; ++i) {
169  for (int j = 0; j < DIM; ++j) {
170  x[j] = data1[i+N*j];
171  }
172  t1->Fill();
173  }
174 
175  // plot the data
176  t1->Draw("x[0]:x[1]:x[2]:x[3]","","candle");
177  TCanvas * c2 = new TCanvas();
178  c2->Divide(3,2);
179  int ic=1;
180  c2->cd(ic++);
181  t1->Draw("x[0]:x[1]");
182  c2->cd(ic++);
183  t1->Draw("x[0]:x[2]");
184  c2->cd(ic++);
185  t1->Draw("x[0]:x[3]");
186  c2->cd(ic++);
187  t1->Draw("x[1]:x[2]");
188  c2->cd(ic++);
189  t1->Draw("x[1]:x[3]");
190  c2->cd(ic++);
191  t1->Draw("x[2]:x[3]");
192 
193  t1->Write();
194  file->Close();
195 
196 }
virtual void SetParameters(const Double_t *params)
Definition: TF1.h:628
float xmin
Definition: THbookFile.cxx:93
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
virtual void SetParName(Int_t ipar, const char *name)
Set name of parameter number ipar.
Definition: TF1.cxx:3356
virtual void Draw(Option_t *option="")=0
Default Draw method for all objects.
void Print(Option_t *option="") const
Print the real and cpu time passed between the start and stop events.
Definition: TStopwatch.cxx:219
void Fatal(const char *location, const char *msgfmt,...)
static void SetDefaultSampler(const char *type)
TVirtualPad * cd(Int_t subpadnumber=0)
Set current canvas & pad.
Definition: TCanvas.cxx:688
TVectorT.
Definition: TMatrixTBase.h:77
const double * Sample()
sample one event and rerturning array x with coordinates
Definition: DistSampler.h:169
#define N
Int_t GetNrows() const
Definition: TVectorT.h:75
#define f(i)
Definition: RSha256.hxx:104
void Print(Option_t *option="") const
Print the vector as a list of elements.
Definition: TVectorT.cxx:1361
TRObject operator()(const T1 &t1) const
double sqrt(double)
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:77
Double_t x[n]
Definition: legend1.C:17
void SetFunction(Function &func, unsigned int dim)
set the parent function distribution to use for sampling (generic case)
Definition: DistSampler.h:72
double pow(double, double)
void Info(const char *location, const char *msgfmt,...)
constexpr Double_t Pi()
Definition: TMath.h:38
void Error(const char *location, const char *msgfmt,...)
virtual Double_t Determinant() const
SVector< double, 2 > v
Definition: Dict.h:5
float xmax
Definition: THbookFile.cxx:93
void SetRange(double xmin, double xmax, int icoord=0)
set range in a given dimension
Definition: DistSampler.cxx:40
The Canvas class.
Definition: TCanvas.h:31
return c2
Definition: legend2.C:14
Interface class for generic sampling of a distribution, i.e.
Definition: DistSampler.h:57
auto * t1
Definition: textangle.C:20
static const std::string & DefaultSampler()
TMatrixTSym< Element > & Similarity(const TMatrixT< Element > &n)
Calculate B * (*this) * B^T , final matrix will be (nrowsb x nrowsb) This is a similarity transform w...
void Print(Option_t *name="") const
Print the matrix as a table of elements.
TMatrixTSym< Element > & Invert(Double_t *det=0)
Invert the matrix and calculate its determinant Notice that the LU decomposition is used instead of B...
virtual void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0)
Automatic pad generation by division.
Definition: TPad.cxx:1162
virtual bool Init(const char *="")
initialize the generators with the given algorithm Implemented by derived classes who needs it (like ...
Definition: DistSampler.h:101
Definition: file.py:1
1-Dim function class
Definition: TF1.h:211
double exp(double)
virtual Double_t EvalPar(const Double_t *x, const Double_t *params=0)
Evaluate function with given coordinates and parameters.
Definition: TF1.cxx:1365
char name[80]
Definition: TGX11.cxx:109
Stopwatch class.
Definition: TStopwatch.h:28