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