Gaudi Framework, version v20r4
Generated: 8 Jan 2009
|
|
|
Generated: 8 Jan 2009 |
00001 // $Id: NumericalDerivative.cpp,v 1.1 2003/11/19 16:56:00 mato Exp $ 00002 // ============================================================================ 00003 // Include files 00004 // ============================================================================ 00005 // GaudiKernel 00006 // ============================================================================ 00007 #include "GaudiKernel/GaudiException.h" 00008 // ============================================================================ 00009 // local 00010 // ============================================================================ 00011 #include "GaudiMath/NumericalDerivative.h" 00012 #include "GaudiMath/Constant.h" 00013 // ============================================================================ 00014 // GSL 00015 // ============================================================================ 00016 #include "gsl/gsl_errno.h" 00017 #include "gsl/gsl_diff.h" 00018 #include "gsl/gsl_math.h" 00019 // ============================================================================ 00020 // local 00021 // ============================================================================ 00022 #include "Helpers.h" 00023 // ============================================================================ 00024 00025 // ============================================================================ 00032 // ============================================================================ 00033 00034 namespace Genfun 00035 { 00036 namespace GaudiMathImplementation 00037 { 00038 // ======================================================================== 00040 // ======================================================================== 00041 FUNCTION_OBJECT_IMP( NumericalDerivative ); 00042 // ======================================================================== 00043 00044 // ======================================================================== 00051 // ======================================================================== 00052 NumericalDerivative::NumericalDerivative 00053 ( const AbsFunction& function , 00054 const size_t index , 00055 const NumericalDerivative::Type& type ) 00056 : AbsFunction () 00057 , m_function ( function.clone() ) 00058 , m_index ( index ) 00059 , m_DIM ( function.dimensionality() ) 00060 , m_type ( type ) 00061 , m_argument ( function.dimensionality() ) 00062 , m_result ( GSL_NEGINF ) 00063 , m_error ( GSL_POSINF ) 00064 { 00065 if( m_index >= m_DIM ) 00066 { Exception ( "::constructor invalid variable index " ) ; }; 00067 } ; 00068 // ======================================================================== 00069 00070 // ======================================================================== 00072 // ======================================================================== 00073 NumericalDerivative::NumericalDerivative 00074 ( const NumericalDerivative& right ) 00075 : AbsFunction () 00076 , m_function ( right.m_function->clone() ) 00077 , m_index ( right.m_index ) 00078 , m_DIM ( right.m_DIM ) 00079 , m_type ( right.m_type ) 00080 , m_argument ( right.m_DIM ) 00081 , m_result ( GSL_NEGINF ) 00082 , m_error ( GSL_POSINF ) 00083 {} ; 00084 // ======================================================================== 00085 00086 // ======================================================================== 00088 // ======================================================================== 00089 NumericalDerivative::~NumericalDerivative() 00090 { if( 0 != m_function ) { delete m_function ; m_function = 0 ; } }; 00091 // ======================================================================== 00092 00093 // ======================================================================== 00095 // ======================================================================== 00096 const NumericalDerivative::Type& 00097 NumericalDerivative::setType ( const NumericalDerivative::Type& value ) 00098 { m_type = value ; return type() ; } 00099 // ======================================================================== 00100 00101 // ======================================================================== 00103 // ======================================================================== 00104 Genfun::Derivative NumericalDerivative::partial( unsigned int idx ) const 00105 { 00106 if( idx >= m_DIM ) 00107 { Exception ( "::partial(i): invalid variable index" ) ; } 00108 const AbsFunction& aux = 00109 NumericalDerivative( *this , idx , type() ) ; 00110 return FunctionNoop( &aux ) ; 00111 }; 00112 // ======================================================================== 00113 00114 // ======================================================================== 00116 // ======================================================================== 00117 double NumericalDerivative::operator() 00118 ( const Argument& argument ) const 00119 { 00120 // reset the result and the error 00121 m_result = GSL_NEGINF ; 00122 m_error = GSL_POSINF ; 00123 00124 // check the argument 00125 if( argument.dimension() != m_DIM ) 00126 { Exception ( "::operator():invalid argument size" ) ; }; 00127 00129 {for( size_t i = 0 ; i < m_DIM ; ++i ){ m_argument[i] = argument[i];}} 00130 00131 // create the helper object 00132 GSL_Helper helper( *m_function , m_argument , m_index ); 00133 00135 gsl_function F ; 00136 F.function = &GSL_Adaptor ; 00137 F.params = &helper ; 00138 00139 double x = argument[m_index]; 00140 int ierrno = 0 ; 00141 switch ( type() ) 00142 { 00143 case Central : 00144 ierrno = gsl_diff_central ( &F , x , &m_result , &m_error ) ; break ; 00145 case Forward : 00146 ierrno = gsl_diff_forward ( &F , x , &m_result , &m_error ) ; break ; 00147 case Backward : 00148 ierrno = gsl_diff_backward ( &F , x , &m_result , &m_error ) ; break ; 00149 default: 00150 Exception ( "::operator(): invalid diffrentiation type " ) ; 00151 } 00152 00153 if( ierrno ) 00154 { gsl_error ( " NumericalDerivative:: the error from gsl_diff_XXXX" , 00155 __FILE__ , __LINE__ , ierrno ) ; } 00156 00157 return result() ; 00158 }; 00159 // ======================================================================== 00160 00161 // ======================================================================== 00163 // ======================================================================== 00164 double NumericalDerivative::operator() ( const double argument ) const 00165 { 00166 // reset the result and the error 00167 m_result = GSL_NEGINF ; 00168 m_error = GSL_POSINF ; 00169 00170 // check the argument 00171 if( 1 != m_DIM ) { Exception ( "operator(): invalid argument size " ); } 00172 00173 Argument arg(1) ; 00174 arg[0] = argument ; 00175 return (*this)( arg ) ; 00176 }; 00177 // ======================================================================== 00178 00179 // ======================================================================== 00181 // ======================================================================== 00182 StatusCode NumericalDerivative::Exception 00183 ( const std::string& message , 00184 const StatusCode& sc ) const 00185 { 00186 throw GaudiException( "NumericalDerivative" + message , 00187 "*GaudiMath*" , sc ) ; 00188 return sc ; 00189 }; 00190 // ======================================================================== 00191 00192 00193 }; // end of namespace GaudiMathImplementation 00194 00195 }; // end of namespace Genfun (due to GSL) 00196 00197 // ============================================================================ 00198 // The END 00199 // ============================================================================