ROOT   6.10/09 Reference Guide
testDerivation.cxx
Go to the documentation of this file.
1 #include "Math/Polynomial.h"
2 #include "Math/Derivator.h"
3 #include "Math/IFunction.h"
4 #include "Math/Functor.h"
5 #include "Math/WrappedFunction.h"
7 #include <iostream>
8 #include <vector>
9 #include <cassert>
10 #include <cmath>
11
12 #ifdef HAVE_ROOTLIBS
13
14 #include "TStopwatch.h"
15 #include "TF1.h"
16 #include "Math/WrappedTF1.h"
17 #include "Math/WrappedMultiTF1.h"
18 #include "Math/DistFunc.h"
19
20 #endif
21
22 const double ERRORLIMIT = 1E-5;
23
24 typedef double ( * FP ) ( double, void * );
25 typedef double ( * FP2 ) ( double );
26
27
28 double myfunc ( double x, void * ) {
29
30  return std::pow( x, 1.5);
31 }
32
33 double myfunc2 ( double x) {
34  return std::pow( x, 1.5);
35 }
36
38
39  int status = 0;
40
41
42  // Derivative of an IGenFunction
43  // Works when compiled c++, compiled ACLiC, interpreted by CINT
45
46  std::vector<double> p(3);
47  p[0] = 2;
48  p[1] = 3;
49  p[2] = 4;
50  f1->SetParameters(&p[0]);
51
53
54  double step = 1E-8;
55  double x0 = 2;
56
57  der->SetFunction(*f1);
58  double result = der->Eval(x0);
59  std::cout << "Derivative of function inheriting from IGenFunction f(x) = 2 + 3x + 4x^2 at x = 2" << std::endl;
60  std::cout << "Return code: " << der->Status() << std::endl;
61  std::cout << "Result: " << result << " +/- " << der->Error() << std::endl;
62  std::cout << "Exact result: " << f1->Derivative(x0) << std::endl;
63  std::cout << "EvalForward: " << der->EvalForward(*f1, x0) << std::endl;
64  std::cout << "EvalBackward: " << der->EvalBackward(x0, step) << std::endl << std::endl;;
65  status += fabs(result-f1->Derivative(x0)) > ERRORLIMIT;
66
67
68  // Derivative of a free function
69  // Works when compiled c++, compiled ACLiC, does not work when interpreted by CINT
70  FP f2 = &myfunc;
71  der->SetFunction(f2);
72
73  std::cout << "Derivative of a free function f(x) = x^(3/2) at x = 2" << std::endl;
74  std::cout << "EvalCentral: " << der->EvalCentral(x0) << std::endl;
75  std::cout << "EvalForward: " << der->EvalForward(x0) << std::endl;
76  std::cout << "EvalBackward: " << der->EvalBackward(x0) << std::endl;
77
78  std::cout << "Exact result: " << 1.5*sqrt(x0) << std::endl << std::endl;
79
80  status += fabs(der->EvalCentral(x0)-1.5*sqrt(x0)) > ERRORLIMIT;
81  status += fabs(der->EvalForward(x0)-1.5*sqrt(x0)) > ERRORLIMIT;
82  status += fabs(der->EvalBackward(x0)-1.5*sqrt(x0)) > ERRORLIMIT;
83
84
85
86  // Derivative of a free function wrapped in an IGenFunction
87  // Works when compiled c++, compiled ACLiC, does not work when interpreted by CINT
89
90  std::cout << "Derivative of a free function wrapped in a Functor f(x) = x^(3/2) at x = 2" << std::endl;
91  std::cout << "EvalCentral: " << der->Eval( *f3, x0) << std::endl;
92  der->SetFunction(*f3);
93  std::cout << "EvalForward: " << der->EvalForward(x0) << std::endl;
94  std::cout << "EvalBackward: " << der->EvalBackward(x0) << std::endl;
95  std::cout << "Exact result: " << 1.5*sqrt(x0) << std::endl << std::endl;
96
97  status += fabs(der->Eval( *f3, x0)-1.5*sqrt(x0)) > ERRORLIMIT;
98  status += fabs(der->EvalForward(x0)-1.5*sqrt(x0)) > ERRORLIMIT;
99  status += fabs(der->EvalBackward(x0)-1.5*sqrt(x0)) > ERRORLIMIT;
100
101
102  // tets case when an empty Derivator is used
103
105  std::cout << "Tes a derivator without a function" << std::endl;
106  std::cout << der2.Eval(1.0) << std::endl;
107
108  // Derivative of a multidim TF1 function
109
110 // #ifdef LATER
111 // TF2 * f2d = new TF2("f2d","x*x + y*y",-10,10,-10,10);
112 // // find gradient at x={1,1}
113 // double vx[2] = {1.,2.};
114 // ROOT::Math::WrappedTF1 fx(*f2d);
115
116 // std::cout << "Derivative of a f(x,y) = x^2 + y^2 at x = 1,y=2" << std::endl;
117 // std::cout << "df/dx = " << der->EvalCentral(fx,1.) << std::endl;
118 // WrappedFunc fy(*f2d,0,vx);
119 // std::cout << "df/dy = " << der->EvalCentral(fy,2.) << std::endl;
120 // #endif
121  delete der; //fix memory leak
122  der = nullptr;
123  return status;
124 }
125
126
127 #ifdef HAVE_ROOTLIBS
128
129 void testDerivPerf() {
130
131
132  std::cout << "\n\n***************************************************************\n";
133  std::cout << "Test derivation performances....\n\n";
134
136  double p[3] = {2,3,4};
137  f1.SetParameters(p);
138
140  int n = 1000000;
141  double x1 = 0; double x2 = 10;
142  double dx = (x2-x1)/double(n);
143
144  timer.Start();
145  double s1 = 0;
146  ROOT::Math::Derivator der(f1);
147  for (int i = 0; i < n; ++i) {
148  double x = x1 + dx*i;
149  s1+= der.EvalCentral(x);
150  }
151  timer.Stop();
152  std::cout << "Time using ROOT::Math::Derivator :\t" << timer.RealTime() << std::endl;
153  int pr = std::cout.precision(18); std::cout << s1 << std::endl; std::cout.precision(pr);
154
155  timer.Start();
156  s1 = 0;
157  for (int i = 0; i < n; ++i) {
158  ROOT::Math::Derivator der2(f1);
159  double x = x1 + dx*i;
160  s1+= der2.EvalForward(x);
161  }
162  timer.Stop();
163  std::cout << "Time using ROOT::Math::Derivator(2):\t" << timer.RealTime() << std::endl;
164  pr = std::cout.precision(18); std::cout << s1 << std::endl; std::cout.precision(pr);
165
166  timer.Start();
167  s1 = 0;
168  for (int i = 0; i < n; ++i) {
169  double x = x1 + dx*i;
170  s1+= ROOT::Math::Derivator::Eval(f1,x);
171  }
172  timer.Stop();
173  std::cout << "Time using ROOT::Math::Derivator(3):\t" << timer.RealTime() << std::endl;
174  pr = std::cout.precision(18); std::cout << s1 << std::endl; std::cout.precision(pr);
175
176
177  TF1 f2("pol","pol2",0,10);
178  f2.SetParameters(p);
179
180  timer.Start();
181  double s2 = 0;
182  for (int i = 0; i < n; ++i) {
183  double x = x1 + dx*i;
184  s2+= f2.Derivative(x);
185  }
186  timer.Stop();
187  std::cout << "Time using TF1::Derivative :\t\t" << timer.RealTime() << std::endl;
188  pr = std::cout.precision(18);
189  std::cout << s2 << std::endl;
190  std::cout.precision(pr);
191
192
193
194 }
195
196 double userFunc(const double *x, const double *) {
197  return std::exp(-x[0]);
198 }
199 double userFunc1(double x) { return userFunc(&x, 0); }
200
201 double userFunc2(const double * x) { return userFunc(x, 0); }
202
203 void testDerivPerfUser() {
204
205
206  std::cout << "\n\n***************************************************************\n";
207  std::cout << "Test derivation performances - using a User function\n\n";
208
210
212  int n = 1000000;
213  double x1 = 0; double x2 = 10;
214  double dx = (x2-x1)/double(n);
215
216  timer.Start();
217  double s1 = 0;
218  ROOT::Math::Derivator der(f1);
219  for (int i = 0; i < n; ++i) {
220  double x = x1 + dx*i;
221  s1+= der.EvalCentral(x);
222  }
223  timer.Stop();
224  std::cout << "Time using ROOT::Math::Derivator :\t" << timer.RealTime() << std::endl;
225  int pr = std::cout.precision(18); std::cout << s1 << std::endl; std::cout.precision(pr);
226
227  timer.Start();
228  s1 = 0;
229  for (int i = 0; i < n; ++i) {
230  ROOT::Math::Derivator der2(f1);
231  double x = x1 + dx*i;
232  s1+= der2.EvalForward(x);
233  }
234  timer.Stop();
235  std::cout << "Time using ROOT::Math::Derivator(2):\t" << timer.RealTime() << std::endl;
236  pr = std::cout.precision(18); std::cout << s1 << std::endl; std::cout.precision(pr);
237
238  timer.Start();
239  s1 = 0;
240  for (int i = 0; i < n; ++i) {
241  double x = x1 + dx*i;
242  s1+= ROOT::Math::Derivator::Eval(f1,x);
243  }
244  timer.Stop();
245  std::cout << "Time using ROOT::Math::Derivator(3):\t" << timer.RealTime() << std::endl;
246  pr = std::cout.precision(18); std::cout << s1 << std::endl; std::cout.precision(pr);
247
248
249  TF1 f2("uf",userFunc,0,10,0);
250
251  timer.Start();
252  double s2 = 0;
253  for (int i = 0; i < n; ++i) {
254  double x = x1 + dx*i;
255  s2+= f2.Derivative(x);
256  }
257  timer.Stop();
258  std::cout << "Time using TF1::Derivative :\t\t" << timer.RealTime() << std::endl;
259  pr = std::cout.precision(18);
260  std::cout << s2 << std::endl;
261  std::cout.precision(pr);
262
263  //typedef double( * FN ) (const double *, const double * );
264  ROOT::Math::WrappedMultiFunction<> f3(userFunc2,1);
265  timer.Start();
266  s1 = 0;
267  double xx[1];
268  for (int i = 0; i < n; ++i) {
269  xx[0] = x1 + dx*i;
270  s1+= ROOT::Math::Derivator::Eval(f3,xx,0);
271  }
272  timer.Stop();
273  std::cout << "Time using ROOT::Math::Derivator Multi:\t" << timer.RealTime() << std::endl;
274  pr = std::cout.precision(18); std::cout << s1 << std::endl; std::cout.precision(pr);
275
276
277
278 }
279
280
281 double gausFunc( const double * x, const double * p) {
282  return p[0] * ROOT::Math::normal_pdf(x[0], p[2], p[1] );
283 }
284
285
286 void testDerivPerfParam() {
287
288
289  std::cout << "\n\n***************************************************************\n";
290  std::cout << "Test derivation performances - using a Gaussian Param function\n\n";
291
292  //TF1 gaus("gaus","gaus",-10,10);
293  TF1 gaus("gaus",gausFunc,-10,10,3);
294  double params[3] = {10,1.,1.};
295  gaus.SetParameters(params);
296
298
300  int n = 300000;
301  double x1 = 0; double x2 = 10;
302  double dx = (x2-x1)/double(n);
303
304  timer.Start();
305  double s1 = 0;
306  for (int i = 0; i < n; ++i) {
307  double x = x1 + dx*i;
308  // param derivatives
309  s1 += ROOT::Math::Derivator::Eval(f1,x,params,0);
310  s1 += ROOT::Math::Derivator::Eval(f1,x,params,1);
311  s1 += ROOT::Math::Derivator::Eval(f1,x,params,2);
312  }
313  timer.Stop();
314  std::cout << "Time using ROOT::Math::Derivator (1D) :\t" << timer.RealTime() << std::endl;
315  int pr = std::cout.precision(18); std::cout << s1 << std::endl; std::cout.precision(pr);
316
317  ROOT::Math::WrappedParamFunction<> f2(&gausFunc,1,params,params+3);
318  double xx[1];
319
320  timer.Start();
321  s1 = 0;
322  for (int i = 0; i < n; ++i) {
323  xx[0] = x1 + dx*i;
324  s1 += ROOT::Math::Derivator::Eval(f2,xx,params,0);
325  s1 += ROOT::Math::Derivator::Eval(f2,xx,params,1);
326  s1 += ROOT::Math::Derivator::Eval(f2,xx,params,2);
327  }
328  timer.Stop();
329  std::cout << "Time using ROOT::Math::Derivator(ND):\t" << timer.RealTime() << std::endl;
330  pr = std::cout.precision(18); std::cout << s1 << std::endl; std::cout.precision(pr);
331
332  // test that func parameters have not been changed
333  assert( std::fabs(params[0] - gaus.GetParameter(0)) < 1.E-15);
334  assert( std::fabs(params[1] - gaus.GetParameter(1)) < 1.E-15);
335  assert( std::fabs(params[2] - gaus.GetParameter(2)) < 1.E-15);
336
337  timer.Start();
338  s1 = 0;
339  double g[3];
340  for (int i = 0; i < n; ++i) {
341  xx[0] = x1 + dx*i;
343  s1 += g[0];
344  s1 += g[1];
345  s1 += g[2];
346  }
347  timer.Stop();
348  std::cout << "Time using TF1::ParamGradient:\t\t" << timer.RealTime() << std::endl;
349  pr = std::cout.precision(18); std::cout << s1 << std::endl; std::cout.precision(pr);
350
351 }
352
353 #endif
354
355 int main() {
356
357  int status = 0;
358
359  status += testDerivation();
360
361 #ifdef HAVE_ROOTLIBS
362  testDerivPerf();
363  testDerivPerfUser();
364  testDerivPerfParam();
365 #endif
366
367  return status;
368
369 }
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
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
Class for computing numerical derivative of a function.
Definition: Derivator.h:61
double(* FP)(double, void *)
Class to Wrap a ROOT Function class (like TF1) in a IParamFunction interface of one dimensions to be ...
Definition: WrappedTF1.h:37
virtual Double_t Derivative(Double_t x, Double_t *params=0, Double_t epsilon=0.001) const
Returns the first derivative of the function at point x, computed by Richardson&#39;s extrapolation metho...
Definition: TF1.cxx:867
Template class to wrap any C++ callable object which takes one argument i.e.
double normal_pdf(double x, double sigma=1, double x0=0)
Probability density function of the normal (Gaussian) distribution.
double sqrt(double)
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
double(* FP2)(double)
double pow(double, double)
int Status() const
return the error status of the last derivative calculation
Definition: Derivator.cxx:156
virtual Double_t GradientPar(Int_t ipar, const Double_t *x, Double_t eps=0.01)
Compute the gradient (derivative) wrt a parameter ipar.
Definition: TF1.cxx:2176
double EvalCentral(double x, double h=1E-8) const
Computes the numerical derivative at a point x using an adaptive central difference algorithm with a ...
Definition: Derivator.cxx:97
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.
void SetFunction(const IGenFunction &f)
Set the function for calculating the derivatives.
Definition: Derivator.cxx:84
double myfunc(double x, void *)
double Eval(double x, double h=1E-8) const
Computes the numerical derivative of a function f at a point x.
Definition: Derivator.cxx:93
int main()
Parametric Function class describing polynomials of order n.
Definition: Polynomial.h:63
constexpr Double_t E()
Definition: TMath.h:74
double myfunc2(double x)
static const double x1[5]
Template class to wrap any C++ callable object implementing operator() (const double * x) in a multi-...
double f2(const double *x)
virtual Double_t GetParameter(Int_t ipar) const
Definition: TF1.h:466
1-Dim function class
Definition: TF1.h:150
TF1 * f1
Definition: legend1.C:11
double result[121]
Functor1D class for one-dimensional functions.
Definition: Functor.h:487
double exp(double)
double Derivative(double x) const
Return the derivative of the function at a point x Use the private method DoDerivative.
Definition: IFunction.h:262
int testDerivation()
const Int_t n
Definition: legend1.C:16
double EvalForward(double x, double h=1E-8) const
Computes the numerical derivative at a point x using an adaptive forward difference algorithm with a ...
Definition: Derivator.cxx:101
WrappedParamFunction class to wrap any multi-dimensional function pbject implementing the operator()(...
const double ERRORLIMIT
Stopwatch class.
Definition: TStopwatch.h:28
double EvalBackward(double x, double h=1E-8) const
Computes the numerical derivative at a point x using an adaptive backward difference algorithm with a...
Definition: Derivator.cxx:105
double Error() const
return the estimate of the absolute error of the last derivative calculation
Definition: Derivator.cxx:154