| The Gaudi Framework v30r0 (c919700c)

- Main Page
- Related Pages
- Modules
- Namespaces
- Classes
- Files

- File List
- File Members

NumericalDerivative.cpp
Go to the documentation of this file.
1 // ============================================================================
2 // Include files
3 // ============================================================================
4 // GaudiKernel
5 // ============================================================================
6 #include "GaudiKernel/GaudiException.h"
7 // ============================================================================
8 // local
9 // ============================================================================
10 #include "GaudiMath/Constant.h"
11 #include "GaudiMath/NumericalDerivative.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 {
35 namespace GaudiMathImplementation
36 {
37 // ========================================================================
38
39 // ========================================================================
46 // ========================================================================
47 NumericalDerivative::NumericalDerivative( const AbsFunction& function, const size_t index,
48 const NumericalDerivative::Type& type )
49 : AbsFunction()
50 , m_function( function.clone() )
51 , m_index( index )
52 , m_DIM( function.dimensionality() )
53 , m_type( type )
54 , m_argument( function.dimensionality() )
55 , m_result( GSL_NEGINF )
56 , m_error( GSL_POSINF )
57 {
58 if ( m_index >= m_DIM ) {
59 Exception( "::constructor invalid variable index " );
60 };
61 }
62 // ========================================================================
63
64 // ========================================================================
66 // ========================================================================
67 NumericalDerivative::NumericalDerivative( const NumericalDerivative& right )
68 : AbsFunction()
69 , m_function( right.m_function->clone() )
70 , m_index( right.m_index )
71 , m_DIM( right.m_DIM )
72 , m_type( right.m_type )
73 , m_argument( right.m_DIM )
74 , m_result( GSL_NEGINF )
75 , m_error( GSL_POSINF )
76 {
77 }
78 // ========================================================================
79
80 // ========================================================================
82 // ========================================================================
83 const NumericalDerivative::Type& NumericalDerivative::setType( const NumericalDerivative::Type& value )
84 {
85 m_type = value;
86 return type();
87 }
88 // ========================================================================
89
90 // ========================================================================
92 // ========================================================================
93 Genfun::Derivative NumericalDerivative::partial( unsigned int idx ) const
94 {
95 if ( idx >= m_DIM ) {
96 Exception( "::partial(i): invalid variable index" );
97 }
98 const AbsFunction& aux = NumericalDerivative( *this, idx, type() );
99 return FunctionNoop( &aux );
100 }
101 // ========================================================================
102
103 // ========================================================================
105 // ========================================================================
106 double NumericalDerivative::operator()( const Argument& argument ) const
107 {
108 // reset the result and the error
109 m_result = GSL_NEGINF;
110 m_error = GSL_POSINF;
111
112 // check the argument
113 if ( argument.dimension() != m_DIM ) {
114 Exception( "::operator():invalid argument size" );
115 };
116
118 {
119 for ( size_t i = 0; i < m_DIM; ++i ) {
120 m_argument[i] = argument[i];
121 }
122 }
123
124 // create the helper object
125 GSL_Helper helper( *m_function, m_argument, m_index );
126
128 gsl_function F;
129 F.function = &GSL_Adaptor;
130 F.params = &helper;
131
132 double x = argument[m_index];
133 int ierrno = 0;
134 switch ( type() ) {
135 case Central:
136 ierrno = gsl_diff_central( &F, x, &m_result, &m_error );
137 break;
138 case Forward:
139 ierrno = gsl_diff_forward( &F, x, &m_result, &m_error );
140 break;
141 case Backward:
142 ierrno = gsl_diff_backward( &F, x, &m_result, &m_error );
143 break;
144 default:
145 Exception( "::operator(): invalid diffrentiation type " );
146 }
147
148 if ( ierrno ) {
149 gsl_error( " NumericalDerivative:: the error from gsl_diff_XXXX", __FILE__, __LINE__, ierrno );
150 }
151
152 return result();
153 }
154 // ========================================================================
155
156 // ========================================================================
158 // ========================================================================
159 double NumericalDerivative::operator()( const double argument ) const
160 {
161 // reset the result and the error
162 m_result = GSL_NEGINF;
163 m_error = GSL_POSINF;
164
165 // check the argument
166 if ( 1 != m_DIM ) {
167 Exception( "operator(): invalid argument size " );
168 }
169
170 Argument arg( 1 );
171 arg[0] = argument;
172 return ( *this )( arg );
173 }
174 // ========================================================================
175
176 // ========================================================================
178 // ========================================================================
179 StatusCode NumericalDerivative::Exception( const std::string& message, const StatusCode& sc ) const
180 {
181 throw GaudiException( "NumericalDerivative" + message, "*GaudiMath*", sc );
182 return sc;
183 }
184 // ========================================================================
185
186 } // end of namespace GaudiMathImplementation
187
188 } // end of namespace Genfun (due to GSL)
189
190 // ============================================================================
191 // The END
192 // ============================================================================
GaudiException
Define general base for Gaudi exception.
Definition: GaudiException.h:22
Genfun::GaudiMathImplementation::NumericalDerivative::type
const Type & type() const
type
Definition: NumericalDerivative.h:118
Genfun::GaudiMathImplementation::NumericalDerivative::Type
Type
the type of numerical differentiation
Definition: NumericalDerivative.h:48
Genfun::GaudiMathImplementation::NumericalDerivative::m_argument
Argument m_argument
Definition: NumericalDerivative.h:143
Genfun::GaudiMathImplementation::NumericalDerivative::m_function
std::unique_ptr< const AbsFunction > m_function
Definition: NumericalDerivative.h:138
Genfun::GaudiMathImplementation::NumericalDerivative::NumericalDerivative
NumericalDerivative()
GaudiMath::Derivative
Genfun::GaudiMathImplementation::NumericalDerivative Derivative
Definition: GaudiMath.h:27
Genfun::GaudiMathImplementation::NumericalDerivative::operator()
double operator()(double argument) const override
Function value.
Definition: NumericalDerivative.cpp:159
Genfun::GaudiMathImplementation::NumericalDerivative::partial
Derivative partial(unsigned int index) const override
Derivatives.
Definition: NumericalDerivative.cpp:93
std::string
STL class.
Genfun::GaudiMathImplementation::NumericalDerivative::m_error
double m_error
Definition: NumericalDerivative.h:145
NumericalDerivative.h
Genfun::GaudiMathImplementation::NumericalDerivative::setType
const Type & setType(const Type &value)
change the type of the adaptive differentiation
Definition: NumericalDerivative.cpp:83
Genfun::GaudiMathImplementation::NumericalDerivative::Forward
Definition: NumericalDerivative.h:48
StatusCode
This class is used for returning status codes from appropriate routines.
Definition: StatusCode.h:26
GaudiException.h
gaudirun.type
type
Definition: gaudirun.py:133
Genfun::GaudiMathImplementation::NumericalDerivative::m_index
size_t m_index
Definition: NumericalDerivative.h:139
Genfun::GaudiMathImplementation::GSL_Adaptor
double GSL_Adaptor(double x, void *params)
Definition: Helpers.cpp:31
Genfun::GaudiMathImplementation::NumericalDerivative::m_DIM
size_t m_DIM
Definition: NumericalDerivative.h:140
Genfun::GaudiMathImplementation::NumericalDerivative
Numerical derivative (using GSL adaptive numerical differentiation)
Definition: NumericalDerivative.h:44
Genfun::GaudiMathImplementation::NumericalDerivative::m_type
Type m_type
Definition: NumericalDerivative.h:141
Genfun::GaudiMathImplementation::NumericalDerivative::Central
Definition: NumericalDerivative.h:48
Genfun::GaudiMathImplementation::NumericalDerivative::m_result
double m_result
Definition: NumericalDerivative.h:144
Helpers.h
Constant.h
Genfun::GaudiMathImplementation::NumericalDerivative::Backward
Definition: NumericalDerivative.h:48
Genfun::GaudiMathImplementation::NumericalDerivative::Exception
StatusCode Exception(const std::string &message, const StatusCode &sc=StatusCode::FAILURE) const
Exception.
Definition: NumericalDerivative.cpp:179
Genfun
CLHEP.
Definition: IEqSolver.h:13
Genfun::GaudiMathImplementation::GSL_Helper
the simple structure to be used for adaption interface Genfun::AbsFunction to gsl_function structure ...
Definition: Helpers.h:25
Genfun::GaudiMathImplementation::NumericalDerivative::result
double result() const
the result of the last call
Definition: NumericalDerivative.h:124