ROOT   6.08/07 Reference Guide
testGSLIntegration.cxx
Go to the documentation of this file.
1 #include "Math/Polynomial.h"
2 #include "Math/Functor.h"
3 #include "Math/Error.h"
4 #include <iostream>
5
6
7 #ifdef HAVE_ROOTLIBS
8 #include "TStopwatch.h"
9 #include "TF1.h"
10 #include "TError.h"
11 #endif
12
13
14 #include "Math/GSLIntegrator.h"
15 // temp before having new Integrator class
16 namespace ROOT {
17  namespace Math {
18  typedef GSLIntegrator Integrator;
19  }
20 }
21
22 const double ERRORLIMIT = 1E-8;
23
24 double exactIntegral ( const std::vector<double> & par, double a, double b) {
25
26  ROOT::Math::Polynomial *func = new ROOT::Math::Polynomial( par.size() +1);
27
28  std::vector<double> p = par;
29  p.push_back(0);
30  p[0] = 0;
31  for (unsigned int i = 1; i < p.size() ; ++i) {
32  p[i] = par[i-1]/double(i);
33  }
34  func->SetParameters(&p.front());
35
36  return (*func)(b)-(*func)(a);
37 }
38
39 double singularFunction(double x) {
40  if (x >= 0)
41  return 1./sqrt(x);
42  else
43  return 1./sqrt(-x);
44 }
45
46
48
49  int status = 0;
50
51
52 #ifdef HAVE_ROOTLIBS
53  gErrorIgnoreLevel = 5000;
54 #endif
55
57
58  std::vector<double> p(3);
59  p[0] = 4;
60  p[1] = 2;
61  p[2] = 6;
62  f->SetParameters(&p[0]);
64
65
66  double exactresult = exactIntegral(p, 0,3);
67  std::cout << "Exact value " << exactresult << std::endl << std::endl;
68
69
70  //ROOT::Math::Integrator ig(func, 0.001, 0.01, 100 );
71  ROOT::Math::GSLIntegrator ig(0.001, 0.01, 100 );
72  ig.SetFunction(func);
73
74
75  double value = ig.Integral( 0, 3);
76  // or ig.Integral(*f, 0, 10); if new function
77
78  std::cout.precision(20);
79
80  std::cout << "Adaptive singular integration:" << std::endl;
81  status &= ig.Status();
82  if (status) std::cout << "Error - Return code " << ig.Status() << std::endl;
83  std::cout << "Result " << value << " +/- " << ig.Error() << std::endl << std::endl;
84  status &= fabs(exactresult-value) > ERRORLIMIT;
85  if (status) MATH_ERROR_MSG("testGSLIntegration","Adaptive integration failed on a polynomial function");
86
87  // integrate again ADAPTIve, with different rule
89  ig2.SetFunction(func);
90  value = ig2.Integral(0, 3);
91  // or ig2.Integral(*f, 0, 10); if different function
92
93  std::cout << "Adaptive Gauss61 integration:" << std::endl;
94  status &= ig.Status();
95  if (status) std::cout << "Error - Return code " << ig.Status() << std::endl;
96  std::cout << "Result " << value << " +/- " << ig2.Error() << std::endl << std::endl;
97  status &= fabs(exactresult-value) > ERRORLIMIT;
98  if (status) MATH_ERROR_MSG("testGSLIntegration","Adaptive GAUSS61 integration failed on a polynomial function");
99
100
101  std::cout << "Testing SetFunction member function" << std::endl;
104
105  pol->SetParameters(&p.front());
106  ROOT::Math::IGenFunction &func2 = *pol;
107  ig3.SetFunction(func2);
108  std::cout << "Result " << ig3.Integral( 0, 3) << " +/- " << ig3.Error() << std::endl;
109  status += fabs(exactresult-ig3.Integral( 0, 3)) > ERRORLIMIT;
110  if (status) MATH_ERROR_MSG("testGSLIntegration","Default Adaptive integration failed on a polynomial function");
111
112  // test error
113  //typedef double ( * FreeFunc ) ( double);
114
115  std::cout << "\nTesting a singular function: 1/sqrt(x)" << std::endl;
116  //ROOT::Math::WrappedFunction<FreeFunc> wf(&singularFunction);
118
119  ig.SetFunction(wf);
120  double r = ig.Integral(0,1);
121  if (ig.Status() != 0)
122  std::cout << "Error integrating a singular function " << std::endl;
123  else
124  std::cout << "Result:(0,1] = " << r << " +/- " << ig.Error() << " (should be 2) " << std::endl;
125  status &= ig.Status();
126  status &= fabs(ig.Result() - 2.0) > ERRORLIMIT;
127  if (status) MATH_ERROR_MSG("testGSLIntegration","Singular Adaptive integration failed on 1./sqrt(x)");
128
129  double singularPts[3] = {-1,0,1};
130  std::vector<double> sp(singularPts, singularPts+3);
131  double r2 = ig.Integral(sp);
132  status &= ig.Status();
133  if (ig.Status() != 0)
134  std::cout << "Error integrating a singular function using vector of points" << std::endl;
135  else
136  std::cout << "Result:[-1,1] = " << r2 << " +/- " << ig.Error() << " (should be 4) " << std::endl;
137  status &= fabs(ig.Result() - 4.0) > ERRORLIMIT;
138  if (status) MATH_ERROR_MSG("testGSLIntegration","Adaptive integration with singular points failed on 1./sqrt(x)");
139
140
141  std::vector<double> sp2(2);
142  sp2[0] = -1.; sp2[1] = -0.5;
143  double r3 = ig.Integral(sp2);
144  std::cout << "Result on [-1,-0.5] = " << r3 << " +/- " << ig.Error() << " (should be 0.5857) " << std::endl;
145  status &= ig.Status();
146  status &= fabs(ig.Result() - (2.-2*sqrt(0.5))) > ERRORLIMIT;
147  if (status) MATH_ERROR_MSG("testGSLIntegration","Adaptive integration with singular points 2 failed on 1./sqrt(x)");
148
149
150  return status;
151 }
152
154
155 #ifdef HAVE_ROOTLIBS
156
157  std::cout << "\n\n***************************************************************\n";
158  std::cout << "Test integration performances....\n\n";
159
160
162  double p[3] = {2,3,4};
163  f1.SetParameters(p);
164
166  int n = 100000;
167  double x1 = 0; double x2 = 10;
168  double dx = (x2-x1)/double(n);
169  double a = -1;
170
171  timer.Start();
173  double s1 = 0;
174  for (int i = 0; i < n; ++i) {
175  double x = x1 + dx*i;
176  s1+= ig.Integral(a,x);
177  }
178  timer.Stop();
179  std::cout << "Time using ROOT::Math::Integrator :\t" << timer.RealTime() << std::endl;
180  int pr = std::cout.precision(18); std::cout << s1 << std::endl; std::cout.precision(pr);
181
182  timer.Start();
183  s1 = 0;
184  for (int i = 0; i < n; ++i) {
185  ROOT::Math::Integrator ig2; ig2.SetFunction(f1);
186  double x = x1 + dx*i;
187  s1+= ig2.Integral(a,x);
188  }
189  timer.Stop();
190  std::cout << "Time using ROOT::Math::Integrator(2):\t" << timer.RealTime() << std::endl;
191  pr = std::cout.precision(18); std::cout << s1 << std::endl; std::cout.precision(pr);
192
193
194  TF1 f2("pol","pol2",0,10);
195  f2.SetParameters(p);
196
197  timer.Start();
198  double s2 = 0;
199  for (int i = 0; i < n; ++i) {
200  double x = x1 + dx*i;
201  s2+= f2.Integral(a,x);
202  }
203  timer.Stop();
204  std::cout << "Time using TF1::Integral :\t\t" << timer.RealTime() << std::endl;
205  pr = std::cout.precision(18); std::cout << s1 << std::endl; std::cout.precision(pr);
206
207 #endif
208
209 }
210
211
212
213 int main() {
214
215  int status = 0;
216
217  status += testIntegration();
218  testIntegPerf();
219
220  return status;
221
222 }
double exactIntegral(const std::vector< double > &par, double a, double b)
double Result() const
return the Result of the last Integral calculation
double par[1]
Definition: unuranDistr.cxx:38
virtual void SetParameters(const Double_t *params)
Definition: TF1.h:439
double singularFunction(double x)
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
Interface (abstract class) for generic functions objects of one-dimension Provides a method to evalua...
Definition: IFunction.h:133
int testIntegration()
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
R__EXTERN Int_t gErrorIgnoreLevel
Definition: TError.h:107
This namespace contains pre-defined functions to be used in conjuction with TExecutor::Map and TExecu...
Definition: StringConv.hxx:21
const double ERRORLIMIT
TArc * a
Definition: textangle.C:12
int Status() const
return the Error Status of the last Integral calculation
virtual Double_t Integral(Double_t a, Double_t b, Double_t epsrel=1.e-12)
IntegralOneDim or analytical integral.
Definition: TF1.cxx:2277
double Integral(const IGenFunction &f, double a, double b)
evaluate the Integral of a function f over the defined interval (a,b)
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
int main()
static const double x2[5]
TStopwatch timer
Definition: pirndm.C:37
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:77
Double_t x[n]
Definition: legend1.C:17
unsigned int r3[N_CITIES]
Definition: simanTSP.cxx:323
Class for performing numerical integration of a function in one dimension.
Definition: GSLIntegrator.h:98
#define MATH_ERROR_MSG(loc, str)
Definition: Error.h:50
double Error() const
return the estimate of the absolute Error of the last Integral calculation
Definition: Integrator.h:420
double Integral(Function &f, double a, double b)
evaluate the Integral of a function f over the defined interval (a,b)
Definition: Integrator.h:506
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
virtual void SetParameters(const double *p)
Set the parameter values.
TRandom2 r(17)
void SetFunction(const IGenFunction &f)
method to set the a generic integration function
User Class for performing numerical integration of a function in one dimension.
Definition: Integrator.h:102
Double_t E()
Definition: TMath.h:54
Parametric Function class describing polynomials of order n.
Definition: Polynomial.h:63
double Error() const
return the estimate of the absolute Error of the last Integral calculation
static const double x1[5]
double f(double x)
double func(double *x, double *p)
Definition: stressTF1.cxx:213
void testIntegPerf()
Namespace for new Math classes and functions.
void SetFunction(Function &f)
method to set the a generic integration function
Definition: Integrator.h:499
double f2(const double *x)
1-Dim function class
Definition: TF1.h:149
TF1 * f1
Definition: legend1.C:11
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
IntegratorOneDim Integrator
Definition: Integrator.h:484
Functor1D class for one-dimensional functions.
Definition: Functor.h:494
const Int_t n
Definition: legend1.C:16
unsigned int r2[N_CITIES]
Definition: simanTSP.cxx:322
Stopwatch class.
Definition: TStopwatch.h:30