Loading [MathJax]/extensions/tex2jax.js
The Gaudi Framework  v31r0 (aeb156f0)
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
NumericalDerivative.cpp
Go to the documentation of this file.
1 // ============================================================================
2 // Include files
3 // ============================================================================
4 // GaudiKernel
5 // ============================================================================
7 // ============================================================================
8 // local
9 // ============================================================================
10 #include "GaudiMath/Constant.h"
12 // ============================================================================
13 // GSL
14 // ============================================================================
15 #include "gsl/gsl_diff.h"
16 #include "gsl/gsl_errno.h"
17 #include "gsl/gsl_math.h"
18 // ============================================================================
19 // local
20 // ============================================================================
21 #include "Helpers.h"
22 // ============================================================================
23 
24 // ============================================================================
31 // ============================================================================
32 
33 namespace Genfun {
34  namespace GaudiMathImplementation {
35  // ========================================================================
36 
37  // ========================================================================
44  // ========================================================================
45  NumericalDerivative::NumericalDerivative( const AbsFunction& function, const size_t index,
47  : AbsFunction()
48  , m_function( function.clone() )
49  , m_index( index )
50  , m_DIM( function.dimensionality() )
51  , m_type( type )
52  , m_argument( function.dimensionality() )
53  , m_result( GSL_NEGINF )
54  , m_error( GSL_POSINF ) {
55  if ( m_index >= m_DIM ) { Exception( "::constructor invalid variable index " ); };
56  }
57  // ========================================================================
58 
59  // ========================================================================
61  // ========================================================================
63  : AbsFunction()
64  , m_function( right.m_function->clone() )
65  , m_index( right.m_index )
66  , m_DIM( right.m_DIM )
67  , m_type( right.m_type )
68  , m_argument( right.m_DIM )
69  , m_result( GSL_NEGINF )
70  , m_error( GSL_POSINF ) {}
71  // ========================================================================
72 
73  // ========================================================================
75  // ========================================================================
77  m_type = value;
78  return type();
79  }
80  // ========================================================================
81 
82  // ========================================================================
84  // ========================================================================
86  if ( idx >= m_DIM ) { Exception( "::partial(i): invalid variable index" ); }
87  const AbsFunction& aux = NumericalDerivative( *this, idx, type() );
88  return FunctionNoop( &aux );
89  }
90  // ========================================================================
91 
92  // ========================================================================
94  // ========================================================================
95  double NumericalDerivative::operator()( const Argument& argument ) const {
96  // reset the result and the error
97  m_result = GSL_NEGINF;
98  m_error = GSL_POSINF;
99 
100  // check the argument
101  if ( argument.dimension() != m_DIM ) { Exception( "::operator():invalid argument size" ); };
102 
104  {
105  for ( size_t i = 0; i < m_DIM; ++i ) { m_argument[i] = argument[i]; }
106  }
107 
108  // create the helper object
110 
112  gsl_function F;
113  F.function = &GSL_Adaptor;
114  F.params = &helper;
115 
116  double x = argument[m_index];
117  int ierrno = 0;
118  switch ( type() ) {
119  case Central:
120  ierrno = gsl_diff_central( &F, x, &m_result, &m_error );
121  break;
122  case Forward:
123  ierrno = gsl_diff_forward( &F, x, &m_result, &m_error );
124  break;
125  case Backward:
126  ierrno = gsl_diff_backward( &F, x, &m_result, &m_error );
127  break;
128  default:
129  Exception( "::operator(): invalid diffrentiation type " );
130  }
131 
132  if ( ierrno ) { gsl_error( " NumericalDerivative:: the error from gsl_diff_XXXX", __FILE__, __LINE__, ierrno ); }
133 
134  return result();
135  }
136  // ========================================================================
137 
138  // ========================================================================
140  // ========================================================================
141  double NumericalDerivative::operator()( const double argument ) const {
142  // reset the result and the error
143  m_result = GSL_NEGINF;
144  m_error = GSL_POSINF;
145 
146  // check the argument
147  if ( 1 != m_DIM ) { Exception( "operator(): invalid argument size " ); }
148 
149  Argument arg( 1 );
150  arg[0] = argument;
151  return ( *this )( arg );
152  }
153  // ========================================================================
154 
155  // ========================================================================
157  // ========================================================================
159  throw GaudiException( "NumericalDerivative" + message, "*GaudiMath*", sc );
160  return sc;
161  }
162  // ========================================================================
163 
164  } // end of namespace GaudiMathImplementation
165 
166 } // namespace Genfun
167 
168 // ============================================================================
169 // The END
170 // ============================================================================
Define general base for Gaudi exception.
Genfun::GaudiMathImplementation::NumericalDerivative Derivative
Definition: GaudiMath.h:26
double operator()(double argument) const override
Function value.
Derivative partial(unsigned int index) const override
Derivatives.
STL class.
const Type & setType(const Type &value)
change the type of the adaptive differentiation
This class is used for returning status codes from appropriate routines.
Definition: StatusCode.h:50
double GSL_Adaptor(double x, void *params)
Definition: Helpers.cpp:25
Numerical derivative (using GSL adaptive numerical differentiation)
StatusCode Exception(const std::string &message, const StatusCode &sc=StatusCode::FAILURE) const
Exception.
CLHEP.
Definition: IEqSolver.h:13
the simple structure to be used for adaption interface Genfun::AbsFunction to gsl_function structure ...
Definition: Helpers.h:23
double result() const
the result of the last call