Logo ROOT   6.10/09
Reference Guide
testIntegrationMultiDim.cxx
Go to the documentation of this file.
1 // Compare the AdaptiveIntegratorMultiDim and
2 // TF1::IntegralMultiple performance and results
3 
4 // Compares time performance
5 // for different dimensions
6 // draws a graph
7 //
8 // Author: David Gonzalez Maline
9 //
10 
11 #include "TMath.h"
12 #include "TStopwatch.h"
13 #include <cmath>
14 #include <iostream>
15 
16 #include "Math/Integrator.h"
17 #include "Math/Functor.h"
18 #include "Math/IFunction.h"
21 #include "Math/IFunctionfwd.h"
22 #include "TF1.h"
23 
24 // for graphical comparison of performance
25 #include "TGraph.h"
26 #include "TAxis.h"
27 #include "TCanvas.h"
28 #include "TApplication.h"
29 #include "TPaveLabel.h"
30 #include "TLegend.h"
31 #include "TH1.h"
32 
33 bool showGraphics = false;
34 bool verbose = false;
35 
36 using namespace std;
37 
38 int NMAX = 6; //maximum dimension
39 
40 Double_t Sum( const double* x, const double *p)
41 {
42  double sum = 0.;
43  for(int i = 0; i < p[0]; i++)
44  sum += x[i];
45  return sum;
46 
47 }
48 
49 //multidim function to integrate
50 Double_t SimpleFun( const double* x, const double *p)
51 {
52  double prod = 1.;
53  for(int i = 0; i < p[0]; i++)
54  prod *= TMath::Power(TMath::E(), -x[i]*x[i]);
55  return prod*Sum(x, p)*TMath::Sin(x[0]);
56 
57 }
58 
59  // ################################################################
60  //
61  // testing IntegratorMultiDim class
62  //
63  // ################################################################
64 
65 double integral_num(unsigned int dim, double* a, double* b, double* p)
66 {
67 
68  if (verbose) std::cout << "\nTesting IntegratorMultiDim class.." << std::endl;
69  //std::cout << p[0] << std::endl;dimensionality
70 
72  timer.Start();
73  ROOT::Math::WrappedParamFunction<> funptr1(&SimpleFun, dim, p, p+1);
74  unsigned int nmax = (unsigned int) 1.E7; // apply cut off to avoid routine to explode
75  ROOT::Math::AdaptiveIntegratorMultiDim ig1(funptr1, 1.E-5, 1.E-5, nmax);
76 
77  ig1.SetFunction(funptr1);
78  ig1.Integral(a, b);
79  timer.Stop();
80 
81  if (verbose) {
82  std::cout.precision(12);
83  std::cout << "result: \t";
84  std::cout << ig1.Result() << "\t" << "error: \t" << ig1.Error() << std::endl;
85  std::cout << "Number of function evaluations: " << ig1.NEval() << std::endl;
86  std::cout << "Time using IntegratorMultiDim: \t" << timer.RealTime() << std::endl;
87  std::cout << "------------------------------------" << std::endl;
88  }
89  return timer.RealTime();
90 }
91 
92  // ################################################################
93  //
94  // testing TF1::IntegralMultiple class
95  //
96  // ################################################################
97 double integral_TF1(unsigned int dim, double* a, double* b, double* p)
98 {
99 
100  double timeTF1;
101  if (verbose) std::cout << "\ntesting TF1::IntegralMultiple.." << std::endl;
102 
104  //timer.Start();
105  //ROOT::Math::WrappedParamFunction<> funptr(&SimpleFun, dim, p, p+1);
106 
107  TF1 function("function", &SimpleFun, 0, 0, 1);
108  function.SetParameters(p);
109  double error, result;
110  int iter, fail;
111  result = function.IntegralMultiple(dim, a, b, 0, (int) 1.E7, 1.E-5, error, iter, fail);
112 
113 // ROOT::Math::GSLMCIntegrator ig1;
114 // ig1.SetType(ROOT::Math::MCIntegration::VEGAS);
115 // ig1.SetFunction(funptr);
116 // ig1.Integral(a, b);
117 
118  timer.Stop();
119 // std::cout << "result: \t";
120 // std::cout << ig1.Result() << "\t" << "error: \t" << ig1.Error() << std::endl;
121 // std::cout << "sigma: \t" << ig1.Sigma();
122 // std::cout << "\t" << "chi2: \t" << ig1.ChiSqr() << std::endl;
123 // std::cout << "\nTime using TF1::IntegralMultiple :\t" << timer.RealTime() << std::endl;
124 // std::cout << "\n------------------------------------" << std::endl;
125 // std::cout << "\t MISER.. \n" << std::endl;
126 
127  if (verbose) {
128  std::cout.precision(12);
129  std::cout << "result: \t";
130  std::cout << result << "\t" << "error: \t" << error << std::endl;
131  std::cout << "Number of function evaluations: " << iter << std::endl;
132  std::cout << "Time using TF1::IntegralMultiple : \t" << timer.RealTime() << std::endl;
133  std::cout << "------------------------------------" << std::endl;
134  }
135 
136  timeTF1 = timer.RealTime();
137 
138  return timeTF1;
139 }
140 
142 {
143  //dimensionality
144  unsigned int Nmax = NMAX;
145  unsigned int size = Nmax-1;
146  TH1D *num_performance = new TH1D("cubature", "", size, 1.5, Nmax+.5);
147  TH1D *TF1_performance = new TH1D("montecarlo", "", size, 1.5, Nmax+.5);
148 
149  num_performance->SetBinContent(1, 0.0);
150  TF1_performance->SetBinContent(1,0.0);
151  std::cout << "Testing multidim integration performances for various dimensions\n";
152  for(unsigned int N = 2; N <=Nmax; N++)//dim
153  {
154  if (verbose) {
155  std::cout<< "*********************************************" << std::endl;
156  std::cout<< "Number of dimensions: "<< N << std::endl;
157  }
158  else {
159  std::cout << "\tdim="<< N << std::endl;
160  }
161  //integration limits
162  double * a = new double[N];
163  double * b = new double[N];
164  double p[1];
165  p[0] = N;
166  for (unsigned int i=0; i < N; i++)
167  {
168  a[i] = -1.;//-TMath::Pi();
169  b[i] = 1;//TMath::Pi();
170  }
171  num_performance->SetBinContent(N-1, integral_num(N, a, b, p));
172  TF1_performance->SetBinContent(N-1,integral_TF1(N, a, b, p));
173  }
174 
175  if (showGraphics) {
176  num_performance->SetBarWidth(0.45);
177  num_performance->SetBarOffset(0.05);
178  num_performance->SetFillColor(49);
179  num_performance->SetStats(0);
180  //num_performance->GetXaxis()->SetLimits(1.5, Nmax+0.5);
181  num_performance->GetXaxis()->SetTitle("number of dimensions");
182  num_performance->GetYaxis()->SetTitle("time [s]");
183  num_performance->SetTitle("comparison of performance");
184  TH1 *h1 = num_performance->DrawCopy("bar2");
185  TF1_performance->SetBarWidth(0.40);
186  TF1_performance->SetBarOffset(0.5);
187  TF1_performance->SetFillColor(kRed);
188  TH1 *h2 = TF1_performance->DrawCopy("bar2,same");
189 
190  TLegend *legend = new TLegend(0.25,0.65,0.55,0.82);
191  legend->AddEntry(h1,"AdaptiveIntegratorMultiDim","f");
192  legend->AddEntry(h2,"TF1::IntegralMultiple()","f");
193  legend->Draw();
194 
195  gPad->Update();
196  }
197  std::cout << "Test Timing results\n";
198  std::cout << "N dim \t Adaptive time \t\t TF1 time \n";
199  for (unsigned int i=1; i<=size; i++)
200  std::cout << i+1 << " " << num_performance->GetBinContent(i) << "\t" << TF1_performance->GetBinContent(i)<<std::endl;
201 
202 }
203 
204 int main(int argc, char **argv)
205 {
206  // Parse command line arguments
207  for (Int_t i=1 ; i<argc ; i++) {
208  std::string arg = argv[i] ;
209  if (arg == "-g") {
210  showGraphics = true;
211  }
212  if (arg == "-v") {
213  showGraphics = true;
214  verbose = true;
215  }
216  if (arg == "-h") {
217  cerr << "Usage: " << argv[0] << " [-g] [-v]\n";
218  cerr << " where:\n";
219  cerr << " -g : graphics mode\n";
220  cerr << " -v : verbose mode";
221  cerr << endl;
222  return -1;
223  }
224  }
225 
226  TApplication* theApp = 0;
227  if ( showGraphics )
228  theApp = new TApplication("App",&argc,argv);
229 
230  performance();
231 
232  if ( showGraphics )
233  {
234  theApp->Run();
235  delete theApp;
236  theApp = 0;
237  }
238 
239  return 0;
240 
241 }
virtual void SetBarOffset(Float_t offset=0.25)
Definition: TH1.h:337
static long int sum(long int i)
Definition: Factory.cxx:2162
virtual void SetParameters(const Double_t *params)
Definition: TF1.h:588
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:23
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
int main(int argc, char **argv)
void SetFunction(const IMultiGenFunction &f)
set the integration function (must implement multi-dim function interface: IBaseFunctionMultiDim) ...
double integral_TF1(unsigned int dim, double *a, double *b, double *p)
double integral_num(unsigned int dim, double *a, double *b, double *p)
Definition: Rtypes.h:56
void performance()
virtual void Draw(Option_t *option="")
Draw this legend with its current attributes.
Definition: TLegend.cxx:452
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4639
bool showGraphics
Double_t SimpleFun(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:2898
STL namespace.
double Result() const
return result of integration
TLegend * legend
Definition: pirndm.C:35
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:628
virtual void SetBarWidth(Float_t width=0.5)
Definition: TH1.h:338
TStopwatch timer
Definition: pirndm.C:37
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:77
Double_t x[n]
Definition: legend1.C:17
virtual void Run(Bool_t retrn=kFALSE)
Main application eventloop. Calls system dependent eventloop via gSystem.
int NEval() const
return number of function evaluations in calculating the integral
TH1F * h1
Definition: legend1.C:5
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)
Set the fill area color.
Definition: TAttFill.h:37
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:8325
double Error() const
return integration error
TAxis * GetYaxis()
Definition: TH1.h:301
tomato 1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:594
constexpr Double_t E()
Definition: TMath.h:74
Double_t Sum(const double *x, const double *p)
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:359
The TH1 histogram class.
Definition: TH1.h:56
1-Dim function class
Definition: TF1.h:150
Double_t Sin(Double_t)
Definition: TMath.h:548
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:310
#define gPad
Definition: TVirtualPad.h:284
double result[121]
This class creates the ROOT Application Environment that interfaces to the windowing system eventloop...
Definition: TApplication.h:39
virtual void SetTitle(const char *title)
Change (i.e.
Definition: TH1.cxx:6028
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:8103
Class for adaptive quadrature integration in multi-dimensions using rectangular regions.
TAxis * GetXaxis()
Definition: TH1.h:300
WrappedParamFunction class to wrap any multi-dimensional function pbject implementing the operator()(...
Stopwatch class.
Definition: TStopwatch.h:28