ROOT  6.06/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"
6 #include "Math/WrappedParamFunction.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 
37 int testDerivation() {
38 
39  int status = 0;
40 
41 
42  // Derivative of an IGenFunction
43  // Works when compiled c++, compiled ACLiC, interpreted by CINT
44  ROOT::Math::Polynomial *f1 = new ROOT::Math::Polynomial(2);
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 
52  ROOT::Math::Derivator *der = new ROOT::Math::Derivator(*f1);
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
88  ROOT::Math::Functor1D *f3 = new ROOT::Math::Functor1D (&myfunc2);
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 
104  ROOT::Math::Derivator der2;
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 
122  return status;
123 }
124 
125 
126 #ifdef HAVE_ROOTLIBS
127 
128 void testDerivPerf() {
129 
130 
131  std::cout << "\n\n***************************************************************\n";
132  std::cout << "Test derivation performances....\n\n";
133 
134  ROOT::Math::Polynomial f1(2);
135  double p[3] = {2,3,4};
136  f1.SetParameters(p);
137 
138  TStopwatch timer;
139  int n = 1000000;
140  double x1 = 0; double x2 = 10;
141  double dx = (x2-x1)/double(n);
142 
143  timer.Start();
144  double s1 = 0;
145  ROOT::Math::Derivator der(f1);
146  for (int i = 0; i < n; ++i) {
147  double x = x1 + dx*i;
148  s1+= der.EvalCentral(x);
149  }
150  timer.Stop();
151  std::cout << "Time using ROOT::Math::Derivator :\t" << timer.RealTime() << std::endl;
152  int pr = std::cout.precision(18); std::cout << s1 << std::endl; std::cout.precision(pr);
153 
154  timer.Start();
155  s1 = 0;
156  for (int i = 0; i < n; ++i) {
157  ROOT::Math::Derivator der2(f1);
158  double x = x1 + dx*i;
159  s1+= der2.EvalForward(x);
160  }
161  timer.Stop();
162  std::cout << "Time using ROOT::Math::Derivator(2):\t" << timer.RealTime() << std::endl;
163  pr = std::cout.precision(18); std::cout << s1 << std::endl; std::cout.precision(pr);
164 
165  timer.Start();
166  s1 = 0;
167  for (int i = 0; i < n; ++i) {
168  double x = x1 + dx*i;
169  s1+= ROOT::Math::Derivator::Eval(f1,x);
170  }
171  timer.Stop();
172  std::cout << "Time using ROOT::Math::Derivator(3):\t" << timer.RealTime() << std::endl;
173  pr = std::cout.precision(18); std::cout << s1 << std::endl; std::cout.precision(pr);
174 
175 
176  TF1 f2("pol","pol2",0,10);
177  f2.SetParameters(p);
178 
179  timer.Start();
180  double s2 = 0;
181  for (int i = 0; i < n; ++i) {
182  double x = x1 + dx*i;
183  s2+= f2.Derivative(x);
184  }
185  timer.Stop();
186  std::cout << "Time using TF1::Derivative :\t\t" << timer.RealTime() << std::endl;
187  pr = std::cout.precision(18);
188  std::cout << s2 << std::endl;
189  std::cout.precision(pr);
190 
191 
192 
193 }
194 
195 double userFunc(const double *x, const double *) {
196  return std::exp(-x[0]);
197 }
198 double userFunc1(double x) { return userFunc(&x, 0); }
199 
200 double userFunc2(const double * x) { return userFunc(x, 0); }
201 
202 void testDerivPerfUser() {
203 
204 
205  std::cout << "\n\n***************************************************************\n";
206  std::cout << "Test derivation performances - using a User function\n\n";
207 
208  ROOT::Math::WrappedFunction<> f1(userFunc1);
209 
210  TStopwatch timer;
211  int n = 1000000;
212  double x1 = 0; double x2 = 10;
213  double dx = (x2-x1)/double(n);
214 
215  timer.Start();
216  double s1 = 0;
217  ROOT::Math::Derivator der(f1);
218  for (int i = 0; i < n; ++i) {
219  double x = x1 + dx*i;
220  s1+= der.EvalCentral(x);
221  }
222  timer.Stop();
223  std::cout << "Time using ROOT::Math::Derivator :\t" << timer.RealTime() << std::endl;
224  int pr = std::cout.precision(18); std::cout << s1 << std::endl; std::cout.precision(pr);
225 
226  timer.Start();
227  s1 = 0;
228  for (int i = 0; i < n; ++i) {
229  ROOT::Math::Derivator der2(f1);
230  double x = x1 + dx*i;
231  s1+= der2.EvalForward(x);
232  }
233  timer.Stop();
234  std::cout << "Time using ROOT::Math::Derivator(2):\t" << timer.RealTime() << std::endl;
235  pr = std::cout.precision(18); std::cout << s1 << std::endl; std::cout.precision(pr);
236 
237  timer.Start();
238  s1 = 0;
239  for (int i = 0; i < n; ++i) {
240  double x = x1 + dx*i;
241  s1+= ROOT::Math::Derivator::Eval(f1,x);
242  }
243  timer.Stop();
244  std::cout << "Time using ROOT::Math::Derivator(3):\t" << timer.RealTime() << std::endl;
245  pr = std::cout.precision(18); std::cout << s1 << std::endl; std::cout.precision(pr);
246 
247 
248  TF1 f2("uf",userFunc,0,10,0);
249 
250  timer.Start();
251  double s2 = 0;
252  for (int i = 0; i < n; ++i) {
253  double x = x1 + dx*i;
254  s2+= f2.Derivative(x);
255  }
256  timer.Stop();
257  std::cout << "Time using TF1::Derivative :\t\t" << timer.RealTime() << std::endl;
258  pr = std::cout.precision(18);
259  std::cout << s2 << std::endl;
260  std::cout.precision(pr);
261 
262  //typedef double( * FN ) (const double *, const double * );
263  ROOT::Math::WrappedMultiFunction<> f3(userFunc2,1);
264  timer.Start();
265  s1 = 0;
266  double xx[1];
267  for (int i = 0; i < n; ++i) {
268  xx[0] = x1 + dx*i;
269  s1+= ROOT::Math::Derivator::Eval(f3,xx,0);
270  }
271  timer.Stop();
272  std::cout << "Time using ROOT::Math::Derivator Multi:\t" << timer.RealTime() << std::endl;
273  pr = std::cout.precision(18); std::cout << s1 << std::endl; std::cout.precision(pr);
274 
275 
276 
277 }
278 
279 
280 double gausFunc( const double * x, const double * p) {
281  return p[0] * ROOT::Math::normal_pdf(x[0], p[2], p[1] );
282 }
283 
284 
285 void testDerivPerfParam() {
286 
287 
288  std::cout << "\n\n***************************************************************\n";
289  std::cout << "Test derivation performances - using a Gaussian Param function\n\n";
290 
291  //TF1 gaus("gaus","gaus",-10,10);
292  TF1 gaus("gaus",gausFunc,-10,10,3);
293  double params[3] = {10,1.,1.};
294  gaus.SetParameters(params);
295 
296  ROOT::Math::WrappedTF1 f1(gaus);
297 
298  TStopwatch timer;
299  int n = 300000;
300  double x1 = 0; double x2 = 10;
301  double dx = (x2-x1)/double(n);
302 
303  timer.Start();
304  double s1 = 0;
305  for (int i = 0; i < n; ++i) {
306  double x = x1 + dx*i;
307  // param derivatives
308  s1 += ROOT::Math::Derivator::Eval(f1,x,params,0);
309  s1 += ROOT::Math::Derivator::Eval(f1,x,params,1);
310  s1 += ROOT::Math::Derivator::Eval(f1,x,params,2);
311  }
312  timer.Stop();
313  std::cout << "Time using ROOT::Math::Derivator (1D) :\t" << timer.RealTime() << std::endl;
314  int pr = std::cout.precision(18); std::cout << s1 << std::endl; std::cout.precision(pr);
315 
316  ROOT::Math::WrappedParamFunction<> f2(&gausFunc,1,params,params+3);
317  double xx[1];
318 
319  timer.Start();
320  s1 = 0;
321  for (int i = 0; i < n; ++i) {
322  xx[0] = x1 + dx*i;
323  s1 += ROOT::Math::Derivator::Eval(f2,xx,params,0);
324  s1 += ROOT::Math::Derivator::Eval(f2,xx,params,1);
325  s1 += ROOT::Math::Derivator::Eval(f2,xx,params,2);
326  }
327  timer.Stop();
328  std::cout << "Time using ROOT::Math::Derivator(ND):\t" << timer.RealTime() << std::endl;
329  pr = std::cout.precision(18); std::cout << s1 << std::endl; std::cout.precision(pr);
330 
331  // test that func parameters have not been changed
332  assert( std::fabs(params[0] - gaus.GetParameter(0)) < 1.E-15);
333  assert( std::fabs(params[1] - gaus.GetParameter(1)) < 1.E-15);
334  assert( std::fabs(params[2] - gaus.GetParameter(2)) < 1.E-15);
335 
336  timer.Start();
337  s1 = 0;
338  double g[3];
339  for (int i = 0; i < n; ++i) {
340  xx[0] = x1 + dx*i;
341  gaus.GradientPar(xx,g,1E-8);
342  s1 += g[0];
343  s1 += g[1];
344  s1 += g[2];
345  }
346  timer.Stop();
347  std::cout << "Time using TF1::ParamGradient:\t\t" << timer.RealTime() << std::endl;
348  pr = std::cout.precision(18); std::cout << s1 << std::endl; std::cout.precision(pr);
349 
350 }
351 
352 #endif
353 
354 int main() {
355 
356  int status = 0;
357 
358  status += testDerivation();
359 
360 #ifdef HAVE_ROOTLIBS
361  testDerivPerf();
362  testDerivPerfUser();
363  testDerivPerfParam();
364 #endif
365 
366  return status;
367 
368 }
TF1::SetParameters
virtual void SetParameters(const Double_t *params)
Definition: TF1.h:439
TStopwatch::RealTime
Double_t RealTime()
Stop the stopwatch (if it is running) and return the realtime (in seconds) passed between the start a...
Definition: TStopwatch.cxx:108
TStopwatch::Start
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:56
ROOT::Math::Derivator
Class for computing numerical derivative of a function.
Definition: Derivator.h:65
FP
double(* FP)(double, void *)
Definition: testDerivation.cxx:24
assert
#define assert(cond)
Definition: unittest.h:542
ROOT::Math::Derivator::Eval
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
ROOT::Math::WrappedTF1
Class to Wrap a ROOT Function class (like TF1) in a IParamFunction interface of one dimensions to be ...
Definition: WrappedTF1.h:41
status
TAlienJobStatus * status
Definition: TAlienJob.cxx:51
ROOT::Math::Derivator::EvalCentral
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
ROOT::Math::WrappedFunction
Template class to wrap any C++ callable object which takes one argument i.e.
Definition: WrappedFunction.h:58
ROOT::Math::normal_pdf
double normal_pdf(double x, double sigma=1, double x0=0)
Probability density function of the normal (Gaussian) distribution.
Definition: PdfFuncMathCore.cxx:274
sqrt
double sqrt(double)
timer
TStopwatch timer
Definition: pirndm.C:37
TStopwatch::Stop
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:75
FP2
double(* FP2)(double)
Definition: testDerivation.cxx:25
pow
double pow(double, double)
ROOT::Math::Derivator::Error
double Error() const
return the estimate of the absolute error of the last derivative calculation
Definition: Derivator.cxx:154
ROOT::Math::fabs
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
Definition: UnaryOperators.h:133
ROOT::Math::ParamFunction::SetParameters
virtual void SetParameters(const double *p)
Set the parameter values.
Definition: ParamFunction.h:104
ROOT::Math::Derivator::SetFunction
void SetFunction(const IGenFunction &f)
Set the function for calculating the derivatives.
Definition: Derivator.cxx:84
myfunc
double myfunc(double x, void *)
Definition: testDerivation.cxx:28
ROOT::Math::Derivator::EvalForward
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
TMath::E
Double_t E()
Definition: TMath.h:54
ROOT::Math::IGradientOneDim::Derivative
double Derivative(double x) const
Return the derivative of the function at a point x Use the private method DoDerivative.
Definition: IFunction.h:258
main
int main()
Definition: testDerivation.cxx:354
ROOT::Math::Polynomial
Parametric Function class describing polynomials of order n.
Definition: Polynomial.h:63
myfunc2
double myfunc2(double x)
Definition: testDerivation.cxx:33
ROOT::Math::Derivator::EvalBackward
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
x1
static const double x1[5]
Definition: RooGaussKronrodIntegrator1D.cxx:326
ROOT::Math::WrappedMultiFunction
Template class to wrap any C++ callable object implementing operator() (const double * x) in a multi-...
Definition: WrappedFunction.h:156
f2
double f2(const double *x)
Definition: testIntegration.cxx:16
TF1
1-Dim function class
Definition: TF1.h:149
f1
TF1 * f1
Definition: legend1.C:11
result
double result[121]
Definition: testSampleQuantiles.cxx:17
ROOT::Math::Functor1D
Functor1D class for one-dimensional functions.
Definition: Functor.h:492
exp
double exp(double)
testDerivation
int testDerivation()
Definition: testDerivation.cxx:37
n
const Int_t n
Definition: legend1.C:16
ROOT::Math::WrappedParamFunction
WrappedParamFunction class to wrap any multi-dimensional function pbject implementing the operator()(...
Definition: WrappedParamFunction.h:43
ROOT::Math::Derivator::Status
int Status() const
return the error status of the last derivative calculation
Definition: Derivator.cxx:156
ERRORLIMIT
const double ERRORLIMIT
Definition: testDerivation.cxx:22
TStopwatch
Stopwatch class.
Definition: TStopwatch.h:30