ROOT   6.10/09 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 #include <iomanip>
6
7
8 #ifdef HAVE_ROOTLIBS
9 #include "TStopwatch.h"
10 #include "TF1.h"
11 #include "TError.h"
12 #endif
13
14
15 #include "Math/GSLIntegrator.h"
16 // temp before having new Integrator class
17 namespace ROOT {
18  namespace Math {
19  typedef GSLIntegrator Integrator;
20  }
21 }
22
23 const double ERRORLIMIT = 1E-8;
24
25 double exactIntegral ( const std::vector<double> & par, double a, double b) {
26
27  ROOT::Math::Polynomial *func = new ROOT::Math::Polynomial( par.size() +1);
28
29  std::vector<double> p = par;
30  p.push_back(0);
31  p[0] = 0;
32  for (unsigned int i = 1; i < p.size() ; ++i) {
33  p[i] = par[i-1]/double(i);
34  }
35  func->SetParameters(&p.front());
36
37  return (*func)(b)-(*func)(a);
38 }
39
40 double singularFunction(double x) {
41  if (x >= 0)
42  return 1./sqrt(x);
43  else
44  return 1./sqrt(-x);
45 }
46
47
49
50  int status = 0;
51
52
53 #ifdef HAVE_ROOTLIBS
54  gErrorIgnoreLevel = 5000;
55 #endif
56
58
59  std::vector<double> p(3);
60  p[0] = 4;
61  p[1] = 2;
62  p[2] = 6;
63  f->SetParameters(&p[0]);
65
66
67  double exactresult = exactIntegral(p, 0,3);
68  std::cout << "Exact value " << exactresult << std::endl << std::endl;
69
70
71  //ROOT::Math::Integrator ig(func, 0.001, 0.01, 100 );
72  ROOT::Math::GSLIntegrator ig(0.001, 0.01, 100 );
73  ig.SetFunction(func);
74
75
76  double value = ig.Integral( 0, 3);
77  // or ig.Integral(*f, 0, 10); if new function
78
79  std::streamsize ss = std::cout.precision();
80  std::cout.precision(20);
81
82  std::cout << "Adaptive singular integration:" << std::endl;
83  status &= ig.Status();
84  if (status) std::cout << "Error - Return code " << ig.Status() << std::endl;
85  std::cout << "Result " << value << " +/- " << ig.Error() << std::endl << std::endl;
86  status &= fabs(exactresult-value) > ERRORLIMIT;
87  if (status) MATH_ERROR_MSG("testGSLIntegration","Adaptive integration failed on a polynomial function");
88
89  // integrate again ADAPTIve, with different rule
91  ig2.SetFunction(func);
92  value = ig2.Integral(0, 3);
93  // or ig2.Integral(*f, 0, 10); if different function
94
95  std::cout << "Adaptive Gauss61 integration:" << std::endl;
96  status &= ig.Status();
97  if (status) std::cout << "Error - Return code " << ig.Status() << std::endl;
98  std::cout << "Result " << value << " +/- " << ig2.Error() << std::endl << std::endl;
99  status &= fabs(exactresult-value) > ERRORLIMIT;
100  if (status) MATH_ERROR_MSG("testGSLIntegration","Adaptive GAUSS61 integration failed on a polynomial function");
101
102
103  std::cout << "Testing SetFunction member function" << std::endl;
106
107  pol->SetParameters(&p.front());
108  ROOT::Math::IGenFunction &func2 = *pol;
109  ig3.SetFunction(func2);
110  std::cout << "Result " << ig3.Integral( 0, 3) << " +/- " << ig3.Error() << std::endl;
111  status += fabs(exactresult-ig3.Integral( 0, 3)) > ERRORLIMIT;
112  if (status) MATH_ERROR_MSG("testGSLIntegration","Default Adaptive integration failed on a polynomial function");
113
114  // test error
115  //typedef double ( * FreeFunc ) ( double);
116
117  std::cout << "\nTesting a singular function: 1/sqrt(x)" << std::endl;
118  //ROOT::Math::WrappedFunction<FreeFunc> wf(&singularFunction);
120
121  ig.SetFunction(wf);
122  double r = ig.Integral(0,1);
123  if (ig.Status() != 0)
124  std::cout << "Error integrating a singular function " << std::endl;
125  else
126  std::cout << "Result:(0,1] = " << r << " +/- " << ig.Error() << " (should be 2) " << std::endl;
127  status &= ig.Status();
128  status &= fabs(ig.Result() - 2.0) > ERRORLIMIT;
129  if (status) MATH_ERROR_MSG("testGSLIntegration","Singular Adaptive integration failed on 1./sqrt(x)");
130
131  double singularPts[3] = {-1,0,1};
132  std::vector<double> sp(singularPts, singularPts+3);
133  double r2 = ig.Integral(sp);
134  status &= ig.Status();
135  if (ig.Status() != 0)
136  std::cout << "Error integrating a singular function using vector of points" << std::endl;
137  else
138  std::cout << "Result:[-1,1] = " << r2 << " +/- " << ig.Error() << " (should be 4) " << std::endl;
139  status &= fabs(ig.Result() - 4.0) > ERRORLIMIT;
140  if (status) MATH_ERROR_MSG("testGSLIntegration","Adaptive integration with singular points failed on 1./sqrt(x)");
141
142
143  std::vector<double> sp2(2);
144  sp2[0] = -1.; sp2[1] = -0.5;
145  double r3 = ig.Integral(sp2);
146  std::cout << "Result on [-1,-0.5] = " << r3 << " +/- " << ig.Error() << " (should be 0.5857) " << std::endl;
147  status &= ig.Status();
148  status &= fabs(ig.Result() - (2.-2*sqrt(0.5))) > ERRORLIMIT;
149  if (status) MATH_ERROR_MSG("testGSLIntegration","Adaptive integration with singular points 2 failed on 1./sqrt(x)");
150
151  std::cout.precision(ss);
152  return status;
153 }
154
156
157 #ifdef HAVE_ROOTLIBS
158
159  std::cout << "\n\n***************************************************************\n";
160  std::cout << "Test integration performances....\n\n";
161
162
164  double p[3] = {2,3,4};
165  f1.SetParameters(p);
166
168  int n = 100000;
169  double x1 = 0; double x2 = 10;
170  double dx = (x2-x1)/double(n);
171  double a = -1;
172
173  timer.Start();
175  double s1 = 0;
176  for (int i = 0; i < n; ++i) {
177  double x = x1 + dx*i;
178  s1+= ig.Integral(a,x);
179  }
180  timer.Stop();
181  std::cout << "Time using ROOT::Math::Integrator :\t" << timer.RealTime() << std::endl;
182  int pr = std::cout.precision(18); std::cout << s1 << std::endl; std::cout.precision(pr);
183
184  timer.Start();
185  s1 = 0;
186  for (int i = 0; i < n; ++i) {
187  ROOT::Math::Integrator ig2; ig2.SetFunction(f1);
188  double x = x1 + dx*i;
189  s1+= ig2.Integral(a,x);
190  }
191  timer.Stop();
192  std::cout << "Time using ROOT::Math::Integrator(2):\t" << timer.RealTime() << std::endl;
193  pr = std::cout.precision(18); std::cout << s1 << std::endl; std::cout.precision(pr);
194
195
196  TF1 f2("pol","pol2",0,10);
197  f2.SetParameters(p);
198
199  timer.Start();
200  double s2 = 0;
201  for (int i = 0; i < n; ++i) {
202  double x = x1 + dx*i;
203  s2+= f2.Integral(a,x);
204  }
205  timer.Stop();
206  std::cout << "Time using TF1::Integral :\t\t" << timer.RealTime() << std::endl;
207  pr = std::cout.precision(18); std::cout << s1 << std::endl; std::cout.precision(pr);
208
209 #endif
210
211 }
212
213
214
215 int main() {
216
217  int status = 0;
218
219  status += testIntegration();
220  testIntegPerf();
221
222  return status;
223
224 }
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:588
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:134
int testIntegration()
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
R__EXTERN Int_t gErrorIgnoreLevel
Definition: TError.h:105
Namespace for new ROOT classes and functions.
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:2309
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:90
#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:412
double Integral(Function &f, double a, double b)
evaluate the Integral of a function f over the defined interval (a,b)
Definition: Integrator.h:496
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:94
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
constexpr Double_t E()
Definition: TMath.h:74
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:489
double f2(const double *x)
1-Dim function class
Definition: TF1.h:150
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:476
Functor1D class for one-dimensional functions.
Definition: Functor.h:487
const Int_t n
Definition: legend1.C:16
unsigned int r2[N_CITIES]
Definition: simanTSP.cxx:322
Stopwatch class.
Definition: TStopwatch.h:28