Logo ROOT   6.08/07
Reference Guide
testRootFinder.cxx
Go to the documentation of this file.
1 #include "TF1.h"
2 #include "Math/Functor.h"
3 #include "TStopwatch.h"
4 #include "RConfigure.h"
5 
6 #include "Math/RootFinder.h"
7 #include "Math/DistFunc.h"
8 
9 #include <iostream>
10 #include <iomanip>
11 #include <cmath>
12 
13 /**
14  Test of ROOT finder for various function
15 
16  case = 0 simple function (polynomial)
17  case = 1 function which fails for a bug in BrentMethod::MinimBrent fixed with r = 32544 (5.27.01)
18  */
19 
20 const double ERRORLIMIT = 1E-8;
21 const int iterTest = 10000;
22 int myfuncCalls = 0;
23 
24 const double Y0_P2 = 5.0; // Y0 for test 1 (parabola)
25 // these are value which gave problems in 5.26 for gamma_cdf
26 const double Y0_GAMMA = 0.32; // Y0 for test 2 (gamma cdf)
27 const double ALPHA_GAMMA = 16.; // alpha of gamma cdf
28 const double THETA_GAMMA = 0.4; // theta of gamma cdf
29 
30 int gTestCase = 0;
31 
32 double myfunc ( double x ) {
33  myfuncCalls ++;
34  if (gTestCase == 0) // polynomial
35  return x*x - Y0_P2;
36  if (gTestCase == 1) // function showing bug in BrentMethod::
38  return 0;
39 }
40 
41 double ExactResult() {
42  if (gTestCase == 0) {
43  return std::sqrt(Y0_P2);
44  }
45  if (gTestCase == 1)
46  //return ROOT::Math::gamma_quantile(Y0_GAMMA,ALPHA_GAMMA,THETA_GAMMA);
47  // put the value to avoid direct MathMore dependency
48  return 5.55680381022934800;
49 
50  return 0;
51 }
52 
53 
54 
55 double myfunc_p (double *x, double *) { return myfunc(x[0]); }
56 
57 int printStats(TStopwatch& timer, double root) {
58 
59  //std::cout << "Return code: " << status << std::endl;
60  double difference = root - ExactResult();
61  int pr = std::cout.precision(16);
62  std::cout << "Result: " << root << std::endl;
63  std::cout << "Exact result: " << ExactResult();
64  std::cout.precision(pr);
65  std::cout << " difference: " << difference << std::endl;
66  std::cout << "Time: " << timer.RealTime()/(double) iterTest << std::endl;
67  std::cout << "Number of calls to function: " << myfuncCalls/iterTest << std::endl;
68 
69  return difference > ERRORLIMIT;
70 }
71 
72 int runTest(int testcase = 0) {
73 
74  double xmin = 0;
75  double xmax = 10;
76 
77  std::cout << "*************************************************************\n";
78  gTestCase = testcase;
79  if (gTestCase == 0)
80  std::cout << "Test for simple function f(x) = x^2 - 5 \n" << std::endl;
81  if (gTestCase == 1) {
82  std::cout << "\nTest for f(x) = gamma_cdf \n" << std::endl;
83  xmin = 3.955687382047723;
84  xmax = 9.3423159494328623;
85  }
86 
88  double root = 0.;
89  int status = 0;
90 
91  double tol = 1.E-14;
92  int maxiter = 100;
93 
95 
96  TF1* f1 = new TF1("f1", myfunc_p, xmin, xmax);
97  timer.Reset(); timer.Start(); myfuncCalls = 0;
98  for (int i = 0; i < iterTest; ++i)
99  {
100  //brf.SetFunction( *func, 0, 10 ); // Just to make a fair comparision!
101  root = f1->GetX(0, xmin, xmax,tol,maxiter);
102  }
103  timer.Stop();
104  std::cout << "\nTF1 Stats:" << std::endl;
105  status += printStats(timer, root);
106 
108  timer.Reset(); timer.Start(); myfuncCalls = 0;
109  for (int i = 0; i < iterTest; ++i)
110  {
111  brf.SetFunction( *func, xmin, xmax );
112  bool ret = brf.Solve(maxiter,tol,tol);
113  if (!ret && i == 0) std::cout << "Error returned from RootFinder::Solve BRENT " << std::endl;
114  root = brf.Root();
115  }
116  timer.Stop();
117  std::cout << "Brent RootFinder Stats:" << std::endl;
118  status += printStats(timer, root);
119 
120 #ifdef R__HAS_MATHMORE
122  timer.Reset(); timer.Start(); myfuncCalls = 0;
123  for (int i = 0; i < iterTest; ++i)
124  {
125  grf.SetFunction( *func, xmin, xmax );
126  bool ret = grf.Solve(maxiter,tol,tol);
127  root = grf.Root();
128  if (!ret && i == 0) std::cout << "Error returned from RootFinder::Solve GSL_BRENT" << std::endl;
129  }
130  timer.Stop();
131  std::cout << "GSL Brent RootFinder Stats:" << std::endl;
132  status += printStats(timer, root);
133 #endif
134 
135  if (status) std::cout << "Test-case " << testcase << " FAILED" << std::endl;
136 
137 
138  return status;
139 }
140 
142 
143  int status = 0;
144  status |= runTest(0); // test pol function
145  if (status) std::cerr << "Test pol function FAILED" << std::endl;
146 
147  status |= runTest(1); // test gamma_cdf
148  if (status) std::cerr << "Test gamma function FAILED" << std::endl;
149 
150  std::cerr << "*************************************************************\n";
151  std::cerr << "\nTest RootFinder :\t";
152  if (status == 0)
153  std::cerr << "OK " << std::endl;
154  else
155  std::cerr << "Failed !" << std::endl;
156 
157  return status;
158 }
159 
160 int main() {
161  return testRootFinder();
162 }
User Class to find the Root of one dimensional functions.
Definition: RootFinder.h:88
int printStats(TStopwatch &timer, double root)
float xmin
Definition: THbookFile.cxx:93
double Root() const
Return the current and latest estimate of the Root.
Definition: RootFinder.h:175
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
double ExactResult()
double myfunc(double x)
const double THETA_GAMMA
const int iterTest
const double Y0_P2
double sqrt(double)
TStopwatch timer
Definition: pirndm.C:37
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:77
Double_t x[n]
Definition: legend1.C:17
int testRootFinder()
double myfunc_p(double *x, double *)
int main()
const double tol
int myfuncCalls
Double_t E()
Definition: TMath.h:54
float xmax
Definition: THbookFile.cxx:93
int runTest(int testcase=0)
bool SetFunction(const IGenFunction &f, double xlow, double xup)
Provide to the solver the function and the initial search interval [xlow, xup] for algorithms not usi...
Definition: RootFinder.h:124
double func(double *x, double *p)
Definition: stressTF1.cxx:213
virtual Double_t GetX(Double_t y, Double_t xmin=0, Double_t xmax=0, Double_t epsilon=1.E-10, Int_t maxiter=100, Bool_t logx=false) const
Returns the X value corresponding to the function value fy for (xmin<x<xmax).
Definition: TF1.cxx:1579
const double Y0_GAMMA
double gamma_cdf(double x, double alpha, double theta, double x0=0)
Cumulative distribution function of the gamma distribution (lower tail).
1-Dim function class
Definition: TF1.h:149
int gTestCase
TF1 * f1
Definition: legend1.C:11
void Reset()
Definition: TStopwatch.h:54
Functor1D class for one-dimensional functions.
Definition: Functor.h:494
const double ERRORLIMIT
Test of ROOT finder for various function.
bool Solve(Function &f, Derivative &d, double start, int maxIter=100, double absTol=1E-8, double relTol=1E-10)
Definition: RootFinder.h:233
const double ALPHA_GAMMA
Stopwatch class.
Definition: TStopwatch.h:30