ROOT  6.06/09
Reference Guide
testMCIntegration.cxx
Go to the documentation of this file.
1 // test of multidimentional Integration
2 // Calculates an intergal of a function
3 // in 2,3,..., 8 dimensions
4 // by using adaptive Genz Malik cubature
5 // and MonteCarlo methods:
6 // --PLAIN
7 // --VEGAS
8 // --MISER
9 //
10 // from
11 // IntegratorMultiDim class
12 // and GSLMCIntegrator class
13 //
14 // Compares time performance
15 // for different dimensions
16 // draws a graph
17 //
18 // Author: Magdalena Slawinska
19 //
20 
21 #include "TMath.h"
22 #include "TStopwatch.h"
23 
24 #include <cmath>
25 #include <iostream>
26 #include <iomanip>
27 
28 #include "Math/Integrator.h"
29 #include "Math/Functor.h"
30 #include "Math/IFunction.h"
33 #include "Math/IFunctionfwd.h"
34 #include "Math/GSLMCIntegrator.h"
35 
36 
37 // for graphical comparison of performance
38 #include "TGraph.h"
39 #include "TAxis.h"
40 #include "TCanvas.h"
41 #include "TApplication.h"
42 #include "TPaveLabel.h"
43 #include "TLegend.h"
44 #include "TH1.h"
45 #include "TCanvas.h"
46 //#include "TLegend.h"
47 
48 bool showGraphics = false;
49 bool verbose = false;
50 int NMAX = 6; // maximum dimension used
51 
52 
53 using std::cout;
54 using std::endl;
55 using std::cerr;
56 
57 //const int n = 3; //default dimensionality
58 
59 Double_t Sum( const double* x, const double *p)
60 {
61  double sum = 0.;
62  for(int i = 0; i < p[0]; i++)
63  sum += x[i];
64  return sum;
65 
66 }
67 
68 //multidim function to integrate
69 Double_t SimpleFun( const double* x, const double *p)
70 {
71  double prod = 1.;
72  for(int i = 0; i < p[0]; i++)
73  prod *= TMath::Power(TMath::E(), -x[i]*x[i]);
74  return prod*Sum(x, p)*TMath::Sin(x[0]);
75 
76 }
77 
78 //singularity at (0,0, ..., 0)
79 Double_t SingularFun( const double* x, const double *p)
80 {
81 
82  double prod = 1.;
83  for(int i = 0; i < p[0]; i++)
84  prod *= TMath::Cos(x[i]);
85  return 1./((1-prod)*8*TMath::Power(TMath::Pi(), 3));
86 
87 }
88 
89 
90  // ################################################################
91  //
92  // testing IntegratorMultiDim class
93  //
94  // ################################################################
95 
96 double integral_num(unsigned int dim, double* a, double* b, double* p, double & value, double & error)
97 {
98  if (verbose) std::cout << "\nTesting IntegratorMultiDim class.." << std::endl;
99  //std::cout << p[0] << std::endl;dimensionality
100 
102  timer.Start();
103  ROOT::Math::WrappedParamFunction<> funptr1(&SimpleFun, dim, p, p+1);
104  unsigned int nmax = (unsigned int) 1.E8; // apply cut off to avoid routine to explode
105  //unsigned int size = (unsigned int) 1.E8; // apply cut off to avoid routine to explode
106  unsigned int size = 0; // use default value defined by nmax
107  ROOT::Math::AdaptiveIntegratorMultiDim ig1(funptr1, 1.E-5, 1.E-5, nmax,size);
108  // std::cout << "1. integral= " << std::endl;
109  ig1.SetFunction(funptr1);
110  ig1.Integral(a, b);
111  timer.Stop();
112  if (verbose) {
113  std::cout.precision(12);
114  std::cout << "result: \t";
115  std::cout << ig1.Result() << "\t" << "error: \t" << ig1.Error() << std::endl;
116  std::cout << "Number of function evaluations: " << ig1.NEval() << std::endl;
117  std::cout << "Time using IntegratorMultiDim: \t" << timer.RealTime() << std::endl;
118  std::cout << "------------------------------------" << std::endl;
119  }
120  else { std::cout << " . "; }
121  value = ig1.Result();
122  error = ig1.Error();
123  return timer.RealTime();
124 }
125 
126  // ################################################################
127  //
128  // testing MCIntegrator class
129  //
130  // ################################################################
131 std::vector<double> integral_MC(unsigned int dim, double* a, double* b, double* p, double * value, double * error)
132 {
133 
134  double timeVegas;
135  double timeMiser;
136  double timePlain;
137  if (verbose) {
138  std::cout << "\nTesting GSLMCIntegrator class.." << std::endl;
139  std::cout << "\t VEGAS.. " << std::endl;
140  std::cout << "" << std::endl;
141  }
142  else { std::cout << "."; }
143 
145  //timer.Start();
146  ROOT::Math::WrappedParamFunction<> funptr(&SimpleFun, dim, p, p+1);
147 
149  ig1.SetFunction(funptr);
150  //ig1.SetMode(ROOT::Math::MCIntegration::kIMPORTANCE_ONLY);
151 
152  /*
153  VegasParameters param;
154  param.iterations = 2;
155  ig1.SetParameters(param);
156  */
157 
158  ig1.Integral(a, b);
159  timer.Stop();
160  timeVegas = timer.RealTime();
161  value[0] = ig1.Result();
162  error[0] = ig1.Error();
163 
164  if (verbose) {
165  std::cout << "result: \t";
166  std::cout << ig1.Result() << "\t" << "error: \t" << ig1.Error() << std::endl;
167  std::cout << "sigma: \t" << ig1.Sigma();
168  std::cout << "\t" << "chi2: \t" << ig1.ChiSqr() << std::endl;
169  std::cout << std::endl;
170  std::cout << "Time using GSLMCIntegrator::VEGAS :\t" << timer.RealTime() << std::endl;
171  //std::cout << "Number of function evaluations: " << ig1.Eval() << std::endl;
172  //ig2_param.iterations = 1000;
173  std::cout << "" <<std::endl;
174  std::cout << "------------------------------------" << std::endl;
175  std::cout << "\t MISER.. " << std::endl;
176  std::cout << "" << std::endl;
177  }
178  else { std::cout << "."; }
179 
180 
181  timer.Start();
183 
184  ig2.SetFunction(funptr);
185 
186 
187  //par.min_calls = 4*dim;
188  //par.min_calls_per_bisection = 8*par.min_calls;
189 
190 
191  //MiserParameters par(dim);
192  //ig2.SetParameters(par);
193  ig2.Integral(a, b);
194  timer.Stop();
195  timeMiser = timer.RealTime();
196  value[1] = ig2.Result();
197  error[1] = ig2.Error();
198  if (verbose) {
199  std::cout << "result: \t";
200  std::cout << ig2.Result() << "\t" << "error: \t" << ig2.Error() << std::endl;
201 
202  std::cout << "Time using GSLMCIntegrator::MISER :\t" << timer.RealTime() << std::endl;
203  std::cout << "" << std::endl;
204  std::cout << "------------------------------------" << std::endl;
205  std::cout << "\t PLAIN.. " << std::endl;
206  }
207  else { std::cout << "."; }
208 
209  timer.Start();
210 
212  ig3.SetFunction(funptr);
213  ig3.Integral(a, b);
214  timer.Stop();
215  timePlain = timer.RealTime();
216  value[2] = ig3.Result();
217  error[2] = ig3.Error();
218 
219  if (verbose) {
220  std::cout << "" << std::endl;
221  std::cout << "result: \t";
222  std::cout << ig3.Result() << "\t" << "error: \t" << ig3.Error() << std::endl;
223  std::cout << "Time using GSLMCIntegrator::PLAIN :\t" << timer.RealTime() << std::endl;
224  std::cout << "" << std::endl;
225  }
226  std::vector<double> result(3);
227  result[0] = timeVegas; result[1] = timeMiser; result[2] = timePlain;
228  return result;
229 }
230 
232 {
233  //dimensionality
234  unsigned int Nmax = NMAX;
235  unsigned int size = Nmax-1;
236  bool ok = true;
237 
238  TH1D *num_performance = new TH1D("cubature", "", size, 1.5, Nmax+.5);
239  TH1D *Vegas_performance = new TH1D("VegasMC", "", size, 1.5, Nmax+.5);
240  TH1D *Miser_performance = new TH1D("MiserMC", "", size, 1.5, Nmax+.5);
241  TH1D *Plain_performance = new TH1D("PlainMC", "", size, 1.5, Nmax+.5);
242 
243  num_performance->SetBinContent(1, 0.0);
244  Vegas_performance->SetBinContent(1,0.0);
245  for(unsigned int N = 2; N <=Nmax; N++)//dim
246  {
247  if (verbose) {
248  std::cout<< "*********************************************" << std::endl;
249  std::cout<< "Number of dimensions: "<< N << std::endl;
250  }
251  else {
252  std::cout << "\n\tdim="<< N << " : ";
253  }
254  //integration limits
255  double * a = new double[N];
256  double * b = new double[N];
257  double p[1];
258  p[0] = N;
259  for (unsigned int i=0; i < N; i++)
260  {
261  a[i] = -1.;//-TMath::Pi();
262  b[i] = 1;//TMath::Pi();
263  }
264  //x[N] = N;
265  double val0, err0;
266  double valMC[3], errMC[3];
267  double timeNumInt = integral_num(N, a, b, p, val0, err0);
268  std::vector<double> timeMCInt = integral_MC(N, a, b, p, valMC, errMC);
269 
270  // set the histograms
271  num_performance->SetBinContent(N-1, timeNumInt);
272  Vegas_performance->SetBinContent(N-1, timeMCInt[0]);
273  Miser_performance->SetBinContent(N-1, timeMCInt[1]);
274  Plain_performance->SetBinContent(N-1, timeMCInt[2]);
275 
276  // test the values
277  for (int j = 0; j < 3; ++j) {
278  if (TMath::Abs(val0-valMC[j] ) > 5 * std::sqrt( err0*err0 + errMC[j]*errMC[j] ) ) {
279  Error("testMCIntegration","Result is not consistent for dim %d between adaptive and MC %d ",N,j);
280  ok = false;
281  }
282  }
283 
284  }
285 
286 
287 
288 
289  if ( showGraphics )
290  {
291 
292  TCanvas * c1 = new TCanvas();
293  c1->SetFillColor(kYellow-10);
294 
295  num_performance->SetBarWidth(0.23);
296  num_performance->SetBarOffset(0.04);
297  num_performance->SetFillColor(kRed+3);
298  num_performance->SetStats(0);
299  //num_performance->GetXaxis()->SetLimits(1.5, Nmax+0.5);
300  num_performance->GetXaxis()->SetTitle("number of dimensions");
301  num_performance->GetYaxis()->SetTitle("time [s]");
302  num_performance->SetTitle("comparison of performance");
303  TH1 *h1 = num_performance->DrawCopy("bar");
304  Vegas_performance->SetBarWidth(0.23);
305  Vegas_performance->SetBarOffset(0.27);
306  Vegas_performance->SetFillColor(kRed);
307  TH1 *h2 = Vegas_performance->DrawCopy("bar,same");
308  Miser_performance->SetBarWidth(0.23);
309  Miser_performance->SetBarOffset(0.50);
310  Miser_performance->SetFillColor(kRed-7);
311  TH1 *h3 = Miser_performance->DrawCopy("bar,same");
312  Plain_performance->SetBarWidth(0.23);
313  Plain_performance->SetBarOffset(0.73);
314  Plain_performance->SetFillColor(kRed-10);
315  TH1 *h4 = Plain_performance->DrawCopy("bar,same");
316 
317  TLegend *legend = new TLegend(0.25,0.65,0.55,0.82);
318  legend->AddEntry(h1,"Cubature","f");
319  legend->AddEntry(h2,"MC Vegas","f");
320  legend->AddEntry(h3,"MC Miser","f");
321  legend->AddEntry(h4,"MC Plain","f");
322  legend->Draw();
323  gPad->SetLogy(true);
324  gPad->Update();
325  }
326 
327  std::cout << "\nTest Timing results\n";
328  std::cout << " N dim \t Adaptive MC Vegas MC Miser MC Plain \n";
329  for (unsigned int i=1; i<=size; i++) {
330  std::cout.width(8);
331  std::cout.precision(6);
332  std::cout << i+1;
333  std::cout << "\t " << std::setw(12) << num_performance->GetBinContent(i) << std::setw(12) << Vegas_performance->GetBinContent(i)
334  << std::setw(12) << Miser_performance->GetBinContent(i) << std::setw(12) << Plain_performance->GetBinContent(i) << std::endl;
335  }
336  return ok;
337 }
338 
339 
340 int main(int argc, char **argv)
341 {
342  int status = 0;
343 
344  // Parse command line arguments
345  for (Int_t i=1 ; i<argc ; i++) {
346  std::string arg = argv[i] ;
347  if (arg == "-g") {
348  showGraphics = true;
349  }
350  if (arg == "-v") {
351  showGraphics = true;
352  verbose = true;
353  }
354  if (arg == "-h") {
355  cout << "Usage: " << argv[0] << " [-g] [-v]\n";
356  cout << " where:\n";
357  cout << " -g : graphics mode\n";
358  cout << " -v : verbose mode";
359  cout << endl;
360  return -1;
361  }
362  }
363 
364 
365  TApplication* theApp = 0;
366 
367  if ( showGraphics )
368  theApp = new TApplication("App",&argc,argv);
369 
370  status = performance() ? 0 : 1;
371 
372  if ( showGraphics )
373  {
374  theApp->Run();
375  delete theApp;
376  theApp = 0;
377  }
378 
379  return status;
380 }
bool performance()
int NEval() const
return number of function evaluations in calculating the integral
virtual void SetBarOffset(Float_t offset=0.25)
Definition: TH1.h:356
Double_t RealTime()
Stop the stopwatch (if it is running) and return the realtime (in seconds) passed between the start a...
Definition: TStopwatch.cxx:108
double Result() const
return the type of the integration used
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4629
This class displays a legend box (TPaveText) containing several legend entries.
Definition: TLegend.h:35
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:56
void SetFunction(const IMultiGenFunction &f)
set the integration function (must implement multi-dim function interface: IBaseFunctionMultiDim) ...
double Integral(const GSLMonteFuncPointer &f, unsigned int dim, double *a, double *b, void *p=0)
evaluate the Integral of a function f over the defined hypercube (a,b)
TCanvas * c1
Definition: legend1.C:2
Definition: Rtypes.h:61
virtual void Draw(Option_t *option="")
Draw this legend with its current attributes.
Definition: TLegend.cxx:373
Double_t Sum(const double *x, const double *p)
double Result() const
return result of integration
#define N
int Int_t
Definition: RtypesCore.h:41
TAlienJobStatus * status
Definition: TAlienJob.cxx:51
TArc * a
Definition: textangle.C:12
Definition: Rtypes.h:61
void SetFunction(const IMultiGenFunction &f)
method to set the a generic integration function
Short_t Abs(Short_t d)
Definition: TMathBase.h:110
int main(int argc, char **argv)
TLegend * legend
Definition: pirndm.C:35
THist< 1, double > TH1D
Definition: THist.h:314
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:501
Double_t SingularFun(const double *x, const double *p)
virtual void SetBarWidth(Float_t width=0.5)
Definition: TH1.h:357
double sqrt(double)
TStopwatch timer
Definition: pirndm.C:37
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:75
Double_t x[n]
Definition: legend1.C:17
double Error() const
return the estimate of the absolute Error of the last Integral calculation
double ChiSqr()
returns chi-squared per degree of freedom for the estimate of the integral in the Vegas algorithm ...
virtual void Run(Bool_t retrn=kFALSE)
Main application eventloop. Calls system dependent eventloop via gSystem.
double integral_num(unsigned int dim, double *a, double *b, double *p, double &value, double &error)
virtual TH1 * DrawCopy(Option_t *option="", const char *name_postfix="_copy") const
Copy this histogram and Draw in the current pad.
Definition: TH1.cxx:2925
TH1F * h1
Definition: legend1.C:5
void Error(const char *location, const char *msgfmt,...)
double Integral(const double *xmin, const double *xmax)
evaluate the integral with the previously given function between xmin[] and xmax[] ...
virtual void SetFillColor(Color_t fcolor)
Definition: TAttFill.h:50
std::vector< double > integral_MC(unsigned int dim, double *a, double *b, double *p, double *value, double *error)
bool showGraphics
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content see convention for numbering bins in TH1::GetBin In case the bin number is greater th...
Definition: TH1.cxx:8543
Double_t E()
Definition: TMath.h:54
TAxis * GetYaxis()
Definition: TH1.h:320
double Error() const
return integration error
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:613
Double_t Cos(Double_t)
Definition: TMath.h:424
Double_t Pi()
Definition: TMath.h:44
The Canvas class.
Definition: TCanvas.h:48
bool verbose
double Double_t
Definition: RtypesCore.h:55
TLegendEntry * AddEntry(const TObject *obj, const char *label="", Option_t *option="lpf")
Add a new entry to this legend.
Definition: TLegend.cxx:280
double Sigma()
set parameters for PLAIN method
The TH1 histogram class.
Definition: TH1.h:80
Double_t Sin(Double_t)
Definition: TMath.h:421
Double_t SimpleFun(const double *x, const double *p)
#define gPad
Definition: TVirtualPad.h:288
double result[121]
This class creates the ROOT Application Environment that interfaces to the windowing system eventloop...
Definition: TApplication.h:45
virtual void SetTitle(const char *title)
Change (i.e.
Definition: TH1.cxx:6268
int NMAX
virtual void SetTitle(const char *title="")
Change (i.e. set) the title of the TNamed.
Definition: TNamed.cxx:152
float value
Definition: math.cpp:443
virtual void SetStats(Bool_t stats=kTRUE)
Set statistics option on/off.
Definition: TH1.cxx:8320
class for adaptive quadrature integration in multi-dimensions using rectangular regions.
TAxis * GetXaxis()
Definition: TH1.h:319
WrappedParamFunction class to wrap any multi-dimensional function pbject implementing the operator()(...
Stopwatch class.
Definition: TStopwatch.h:30