19 _arg1(right._arg1->clone()),
20 _wrtIndex(right._wrtIndex)
34 #define ROBUST_DERIVATIVES
35 #ifdef ROBUST_DERIVATIVES
37 double FunctionNumDeriv::f_x (
double x)
const {
42 double FunctionNumDeriv::f_Arg (
double x)
const {
43 _xArg [_wrtIndex] = x;
44 return (*_arg1)(_xArg);
50 assert (_wrtIndex==0);
51 return numericalDerivative ( &FunctionNumDeriv::f_x, x );
58 double xx = x[_wrtIndex];
59 return numericalDerivative ( &FunctionNumDeriv::f_Arg, xx );
63 double FunctionNumDeriv::numericalDerivative
66 const double h0 = 5 * pow(2.0, -17);
68 const double maxErrorA = .0012;
69 const double maxErrorB = .0000026;
71 const double maxErrorC = .0003;
75 const int nItersMax = 6;
77 double bestError = 1.0E30;
80 const double valFactor = pow(2.0, -16);
82 const double w = 5.0/8;
83 const double wi2 = 64.0/25.0;
84 const double wi4 = wi2*wi2;
86 double size = fabs((this->*f)(x));
87 if (size==0) size = pow(2.0, -53);
89 const double adjustmentFactor[nItersMax] = {
97 for ( nIters = 0; nIters < nItersMax; ++nIters ) {
99 double h = h0 * adjustmentFactor[nIters];
103 double A1 = ((this->*
f)(x+h) - (this->*
f)(x-h))/(2.0*h);
105 if (fabs(A1) > size) size = fabs(A1);
108 double A2 = ((this->*
f)(x+hh) - (this->*
f)(x-hh))/(2.0*hh);
110 if (fabs(A2) > size) size = fabs(A2);
113 double A3 = ((this->*
f)(x+hh) - (this->*
f)(x-hh))/(2.0*hh);
115 if (fabs(A3) > size) size = fabs(A3);
117 if ( (fabs(A1-A2)/size > maxErrorA) || (fabs(A1-A3)/size > maxErrorA) ) {
123 double B1 = ( A2 * wi2 - A1 ) / ( wi2 - 1 );
124 double B2 = ( A3 * wi2 - A2 ) / ( wi2 - 1 );
125 if ( fabs(B1-B2)/size > maxErrorB ) {
131 double ans = ( B2 * wi4 - B1 ) / ( wi4 - 1 );
132 double err = fabs ( ans - B1 );
133 if ( err < bestError ) {
141 double val = ((this->*
f)(x+hh) - (this->*
f)(x-hh))/(2.0*hh);
142 if ( fabs(val-ans)/size > maxErrorC ) {
153 #endif // ROBUST_DERIVATIVES
157 #ifdef SIMPLER_DERIVATIVES
160 assert (_wrtIndex==0);
161 const double h=1.0E-6;
162 return ((*_arg1)(x+h) - (*_arg1)(x-h))/(2.0*h);
167 assert (_wrtIndex<x.dimension());
168 const double h=1.0E-6;
172 return ((*_arg1)(x1) - (*_arg1)(x0))/(2.0*h);
174 #endif // SIMPLER_DERIVATIVES