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 
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
Genfun::AbsFunction
Definition: CLHEP/GenericFunctions/AbsFunction.hh:48
Genfun::FunctionNumDeriv
Definition: CLHEP/GenericFunctions/FunctionNumDeriv.hh:19
Genfun::FunctionNumDeriv::FunctionNumDeriv
FunctionNumDeriv(const AbsFunction *arg1, unsigned int index=0)
Definition: FunctionNumDeriv.cc:12
Genfun::FunctionNumDeriv::operator()
virtual double operator()(double argument) const
Definition: FunctionNumDeriv.cc:49
size
user code seldom needs to call this function directly ZMerrno whether or not they are still recorded ZMerrno size() Return the(integer) number of ZMthrow 'n exceptions currently recorded. 5) ZMerrno.clear() Set an internal counter to zero. This counter is available(see next function) to user code to track ZMthrow 'n exceptions that have occurred during any arbitrary time interval. 6) ZMerrno.countSinceCleared() Return the(integer) number of ZMthrow 'n exceptions that have been recorded via ZMerrno.write()
f
void f(void g())
Definition: excDblThrow.cc:38
Genfun::AbsFunction::dimensionality
virtual unsigned int dimensionality() const
Definition: AbsFunction.cc:79
Genfun::FunctionNumDeriv::dimensionality
virtual unsigned int dimensionality() const
Definition: FunctionNumDeriv.cc:31
Genfun::FunctionNumDeriv::~FunctionNumDeriv
virtual ~FunctionNumDeriv()
Definition: FunctionNumDeriv.cc:26
Genfun::Argument
Definition: CLHEP/GenericFunctions/Argument.hh:17
FunctionNumDeriv.hh
x
any side effects of that construction would occur twice The semantics of throw x
Definition: whyZMthrowRethrows.txt:37
FUNCTION_OBJECT_IMP
#define FUNCTION_OBJECT_IMP(classname)
Definition: CLHEP/GenericFunctions/AbsFunction.hh:156
Genfun
Definition: CLHEP/GenericFunctions/Abs.hh:14