Logo ROOT  
Reference Guide
multidimSampling.C File Reference

Detailed Description

View in nbviewer Open in SWAN
Example of sampling a multi-dim distribution using the DistSampler class NOTE: This tutorial must be run with ACLIC

// function (a 4d gaussian)
#include "TMath.h"
#include "TF2.h"
#include "TStopwatch.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;
// Define the GausND strcture
struct GausND {
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;
}
};
// Use the Math namespace
using namespace ROOT::Math;
void multidimSampling() {
const int N = 10000;
/*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;
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 ",
}
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];
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();
}
// 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();
}
#define f(i)
Definition: RSha256.hxx:104
void Info(const char *location, const char *msgfmt,...)
Use this function for informational messages.
Definition: TError.cxx:221
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition: TError.cxx:188
void Fatal(const char *location, const char *msgfmt,...)
Use this function in case of a fatal error. It will abort the program.
Definition: TError.cxx:247
#define N
winID h TVirtualViewer3D TVirtualGLPainter p
char name[80]
Definition: TGX11.cxx:110
float xmin
Definition: THbookFile.cxx:95
float xmax
Definition: THbookFile.cxx:95
TRObject operator()(const T1 &t1) const
static void SetDefaultSampler(const char *type)
static const std::string & DefaultSampler()
Interface class for generic sampling of a distribution, i.e.
Definition: DistSampler.h:57
const double * Sample()
Sample one event and return an array x with sample coordinates values.
Definition: DistSampler.h:193
void SetRange(double xmin, double xmax, int icoord=0)
Set the range in a given dimension.
Definition: DistSampler.cxx:40
void SetFunction(Function &func, unsigned int dim)
set the parent function distribution to use for sampling (generic case)
Definition: DistSampler.h:72
virtual bool Init(const char *="")
Initialize the sampling generator with the given algorithm.
Definition: DistSampler.h:101
The Canvas class.
Definition: TCanvas.h:23
1-Dim function class
Definition: TF1.h:213
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format.
Definition: TFile.h:54
void Print(Option_t *name="") const override
Print the matrix as a table of elements.
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...
Double_t Determinant() const override
TMatrixTSym< Element > & Invert(Double_t *det=nullptr)
Invert the matrix and calculate its determinant Notice that the LU decomposition is used instead of B...
Stopwatch class.
Definition: TStopwatch.h:28
Basic string class.
Definition: TString.h:136
A TTree represents a columnar dataset.
Definition: TTree.h:79
Int_t GetNrows() const
Definition: TVectorT.h:75
void Print(Option_t *option="") const override
Print the vector as a list of elements.
Definition: TVectorT.cxx:1363
RVec< PromoteTypes< T0, T1 > > pow(const T0 &x, const RVec< T1 > &v)
Definition: RVec.hxx:1792
RVec< PromoteType< T > > exp(const RVec< T > &v)
Definition: RVec.hxx:1783
Double_t x[n]
Definition: legend1.C:17
return c2
Definition: legend2.C:14
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
constexpr Double_t Pi()
Definition: TMath.h:37
Definition: file.py:1
auto * t1
Definition: textangle.C:20
Author
Lorenzo Moneta

Definition in file multidimSampling.C.