CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

FunctionNumDeriv.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 // $Id: FunctionNumDeriv.cc,v 1.4 2003/10/10 17:40:39 garren Exp $
3 // ---------------------------------------------------------------------------
4 
6 #include <assert.h>
7 #include <cmath> // for pow()
8 
9 namespace Genfun {
10 FUNCTION_OBJECT_IMP(FunctionNumDeriv)
11 
12 FunctionNumDeriv::FunctionNumDeriv(const AbsFunction *arg1, unsigned int index):
13  _arg1(arg1->clone()),
14  _wrtIndex(index)
15 {
16 }
17 
19  AbsFunction(right),
20  _arg1(right._arg1->clone()),
21  _wrtIndex(right._wrtIndex)
22 {
23 }
24 
25 
27 {
28  delete _arg1;
29 }
30 
31 unsigned int FunctionNumDeriv::dimensionality() const {
32  return _arg1->dimensionality();
33 }
34 
35 #define ROBUST_DERIVATIVES
36 #ifdef ROBUST_DERIVATIVES
37 
38 double FunctionNumDeriv::f_x (double x) const {
39  return (*_arg1)(x);
40 }
41 
42 
43 double FunctionNumDeriv::f_Arg (double x) const {
44  _xArg [_wrtIndex] = x;
45  return (*_arg1)(_xArg);
46 }
47 
48 
49 double FunctionNumDeriv::operator ()(double x) const
50 {
51  assert (_wrtIndex==0);
52  return numericalDerivative ( &FunctionNumDeriv::f_x, x );
53 }
54 
55 double FunctionNumDeriv::operator ()(const Argument & x) const
56 {
57  assert (_wrtIndex<x.dimension());
58  _xArg = x;
59  double xx = x[_wrtIndex];
60  return numericalDerivative ( &FunctionNumDeriv::f_Arg, xx );
61 }
62 
63 
64 double FunctionNumDeriv::numericalDerivative
65  ( double (FunctionNumDeriv::*f)(double)const, double x ) const {
66 
67  const double h0 = 5 * std::pow(2.0, -17);
68 
69  const double maxErrorA = .0012; // These are the largest errors in steps
70  const double maxErrorB = .0000026; // A, B consistent with 8-digit accuracy.
71 
72  const double maxErrorC = .0003; // Largest acceptable validation discrepancy.
73 
74  // This value of gives 8-digit accuracy for 1250 > curvature scale < 1/1250.
75 
76  const int nItersMax = 6;
77  int nIters;
78  double bestError = 1.0E30;
79  double bestAns = 0;
80 
81  const double valFactor = std::pow(2.0, -16);
82 
83  const double w = 5.0/8;
84  const double wi2 = 64.0/25.0;
85  const double wi4 = wi2*wi2;
86 
87  double size = fabs((this->*f)(x));
88  if (size==0) size = std::pow(2.0, -53);
89 
90  const double adjustmentFactor[nItersMax] = {
91  1.0,
92  std::pow(2.0, -17),
93  std::pow(2.0, +17),
94  std::pow(2.0, -34),
95  std::pow(2.0, +34),
96  std::pow(2.0, -51) };
97 
98  for ( nIters = 0; nIters < nItersMax; ++nIters ) {
99 
100  double h = h0 * adjustmentFactor[nIters];
101 
102  // Step A: Three estimates based on h and two smaller values:
103 
104  double A1 = ((this->*f)(x+h) - (this->*f)(x-h))/(2.0*h);
105 // size = max(fabs(A1), size);
106  if (fabs(A1) > size) size = fabs(A1);
107 
108  double hh = w*h;
109  double A2 = ((this->*f)(x+hh) - (this->*f)(x-hh))/(2.0*hh);
110 // size = max(fabs(A2), size);
111  if (fabs(A2) > size) size = fabs(A2);
112 
113  hh *= w;
114  double A3 = ((this->*f)(x+hh) - (this->*f)(x-hh))/(2.0*hh);
115 // size = max(fabs(A3), size);
116  if (fabs(A3) > size) size = fabs(A3);
117 
118  if ( (fabs(A1-A2)/size > maxErrorA) || (fabs(A1-A3)/size > maxErrorA) ) {
119  continue;
120  }
121 
122  // Step B: Two second-order estimates based on h h*w, from A estimates
123 
124  double B1 = ( A2 * wi2 - A1 ) / ( wi2 - 1 );
125  double B2 = ( A3 * wi2 - A2 ) / ( wi2 - 1 );
126  if ( fabs(B1-B2)/size > maxErrorB ) {
127  continue;
128  }
129 
130  // Step C: Third-order estimate, from B estimates:
131 
132  double ans = ( B2 * wi4 - B1 ) / ( wi4 - 1 );
133  double err = fabs ( ans - B1 );
134  if ( err < bestError ) {
135  bestError = err;
136  bestAns = ans;
137  }
138 
139  // Validation estimate based on much smaller h value:
140 
141  hh = h * valFactor;
142  double val = ((this->*f)(x+hh) - (this->*f)(x-hh))/(2.0*hh);
143  if ( fabs(val-ans)/size > maxErrorC ) {
144  continue;
145  }
146 
147  // Having passed both apparent accuracy and validation, we are finished:
148  break;
149  }
150 
151  return bestAns;
152 
153 }
154 #endif // ROBUST_DERIVATIVES
155 
156 
157 
158 #ifdef SIMPLER_DERIVATIVES
159 double FunctionNumDeriv::operator ()(double x) const
160 {
161  assert (_wrtIndex==0);
162  const double h=1.0E-6;
163  return ((*_arg1)(x+h) - (*_arg1)(x-h))/(2.0*h);
164 }
165 
166 double FunctionNumDeriv::operator ()(const Argument & x) const
167 {
168  assert (_wrtIndex<x.dimension());
169  const double h=1.0E-6;
170  Argument x1=x, x0=x;
171  x1[_wrtIndex] +=h;
172  x0[_wrtIndex] -=h;
173  return ((*_arg1)(x1) - (*_arg1)(x0))/(2.0*h);
174 }
175 #endif // SIMPLER_DERIVATIVES
176 
177 } // namespace Genfun
virtual unsigned int dimensionality() const
FunctionNumDeriv(const AbsFunction *arg1, unsigned int index=0)
virtual unsigned int dimensionality() const
Definition: AbsFunction.cc:79
void f(void g())
Definition: excDblThrow.cc:38
unsigned int dimension() const
#define FUNCTION_OBJECT_IMP(classname)
virtual double operator()(double argument) const