Gaudi Framework, version v20r4
Generated: 8 Jan 2009
|
|
|
Generated: 8 Jan 2009 |
00001 // $Id: NumericalIndefiniteIntegral.cpp,v 1.2 2007/05/24 14:39:11 hmd Exp $ 00002 // ============================================================================ 00003 // CVS tag $Name: $ 00004 // ============================================================================ 00005 // Include 00006 // ============================================================================ 00007 // STD & STL 00008 // ============================================================================ 00009 #include <vector> 00010 #include <algorithm> 00011 // ============================================================================ 00012 // GSL 00013 // ============================================================================ 00014 #include "gsl/gsl_errno.h" 00015 #include "gsl/gsl_integration.h" 00016 // ============================================================================ 00017 // GaudiMath 00018 // ============================================================================ 00019 #include "GaudiMath/NumericalIndefiniteIntegral.h" 00020 #include "GaudiMath/NumericalDerivative.h" 00021 // ============================================================================ 00022 // GaudiKernel 00023 // ============================================================================ 00024 #include "GaudiKernel/GaudiException.h" 00025 // ============================================================================ 00026 // local 00027 // ============================================================================ 00028 #include "Helpers.h" 00029 // ============================================================================ 00030 00031 // ============================================================================ 00036 // ============================================================================ 00037 00038 namespace Genfun 00039 { 00040 namespace GaudiMathImplementation 00041 { 00042 00043 struct NumericalIndefiniteIntegral::_Workspace 00044 { gsl_integration_workspace* ws ; }; 00045 struct NumericalIndefiniteIntegral::_Function 00046 { gsl_function* fn ; }; 00047 00048 // ======================================================================== 00050 // ======================================================================== 00051 FUNCTION_OBJECT_IMP( NumericalIndefiniteIntegral ); 00052 // ======================================================================== 00053 00054 // ======================================================================== 00069 // ======================================================================== 00070 NumericalIndefiniteIntegral::NumericalIndefiniteIntegral 00071 ( const AbsFunction& function , 00072 const size_t index , 00073 const double a , 00074 const GaudiMath::Integration::Limit limit , 00075 const GaudiMath::Integration::Type type , 00076 const GaudiMath::Integration::KronrodRule rule , 00077 const double epsabs , 00078 const double epsrel , 00079 const size_t size ) 00080 : AbsFunction () 00081 , m_function ( function.clone() ) 00082 , m_DIM ( function.dimensionality() ) 00083 , m_index ( index ) 00084 , m_a ( a ) 00085 , m_limit ( limit ) 00086 , m_type ( type ) 00087 , m_category ( GaudiMath::Integration::Finite ) 00088 , m_rule ( rule ) 00089 // 00090 , m_points ( ) 00091 , m_pdata ( 0 ) 00092 // 00093 , m_epsabs ( epsabs ) 00094 , m_epsrel ( epsrel ) 00095 // 00096 , m_result ( GSL_NEGINF ) 00097 , m_error ( GSL_POSINF ) 00098 // 00099 , m_size ( size ) 00100 , m_ws ( 0 ) 00101 , m_argument ( function.dimensionality() ) 00102 { 00103 if ( GaudiMath::Integration::Fixed == m_rule ) 00104 { m_rule = GaudiMath::Integration::Default ; } 00105 if ( m_index >= m_DIM ) 00106 { Exception("::constructor: invalid variable index") ; } 00107 }; 00108 // ======================================================================== 00109 00120 NumericalIndefiniteIntegral::NumericalIndefiniteIntegral 00121 ( const AbsFunction& function , 00122 const size_t index , 00123 const double a , 00124 const Points& points , 00125 const GaudiMath::Integration::Limit limit , 00126 const double epsabs , 00127 const double epsrel , 00128 const size_t size ) 00129 : AbsFunction () 00130 , m_function ( function.clone() ) 00131 , m_DIM ( function.dimensionality() ) 00132 , m_index ( index ) 00133 , m_a ( a ) 00134 , m_limit ( limit ) 00135 , m_type ( GaudiMath::Integration:: Other ) 00136 , m_category ( GaudiMath::Integration:: Singular ) 00137 , m_rule ( GaudiMath::Integration:: Fixed ) 00138 , m_points ( points ) 00139 , m_pdata ( 0 ) 00140 , m_epsabs ( epsabs ) 00141 , m_epsrel ( epsrel ) 00142 // 00143 , m_result ( GSL_NEGINF ) 00144 , m_error ( GSL_POSINF ) 00145 // 00146 , m_size ( size ) 00147 , m_ws ( 0 ) 00148 , m_argument ( function.dimensionality() ) 00149 { 00150 if ( m_index >= m_DIM ) 00151 { Exception("::constructor: invalid variable index") ; } 00152 m_pdata = new double[ 2 + m_points.size() ] ; 00153 m_points.push_back( a ) ; 00154 std::sort( m_points.begin() , m_points.end() ) ; 00155 m_points.erase ( std::unique( m_points.begin () , 00156 m_points.end () ) , m_points.end() ); 00157 }; 00158 00159 // ======================================================================== 00167 // ======================================================================== 00168 NumericalIndefiniteIntegral::NumericalIndefiniteIntegral 00169 ( const AbsFunction& function , 00170 const size_t index , 00171 const GaudiMath::Integration::Limit limit , 00172 const double epsabs , 00173 const double epsrel , 00174 const size_t size ) 00175 : AbsFunction () 00176 , m_function ( function.clone() ) 00177 , m_DIM ( function.dimensionality() ) 00178 , m_index ( index ) 00179 , m_a ( GSL_NEGINF ) // should not be used! 00180 , m_limit ( limit ) 00181 , m_type ( GaudiMath::Integration:: Other ) 00182 , m_category ( GaudiMath::Integration:: Infinite ) 00183 , m_rule ( GaudiMath::Integration:: Fixed ) 00184 , m_points ( ) 00185 , m_pdata ( 0 ) 00186 , m_epsabs ( epsabs ) 00187 , m_epsrel ( epsrel ) 00188 , m_result ( GSL_NEGINF ) 00189 , m_error ( GSL_POSINF ) 00190 , m_size ( size ) 00191 , m_ws ( 0 ) 00192 , m_argument ( function.dimensionality() ) 00193 { 00194 if ( m_index >= m_DIM ) 00195 { Exception("::constructor: invalid variable index") ; } 00196 }; 00197 // ======================================================================== 00198 00199 00200 // ======================================================================== 00202 // ======================================================================== 00203 NumericalIndefiniteIntegral:: 00204 NumericalIndefiniteIntegral( const NumericalIndefiniteIntegral& right ) 00205 : AbsFunction () 00206 , m_function ( right.m_function->clone() ) 00207 , m_DIM ( right.m_DIM ) 00208 , m_index ( right.m_index ) 00209 , m_a ( right.m_a ) 00210 , m_limit ( right.m_limit ) 00211 , m_type ( right.m_type ) 00212 , m_category ( right.m_category ) 00213 , m_rule ( right.m_rule ) 00214 , m_points ( right.m_points ) 00215 , m_pdata ( 0 ) // attention 00216 , m_epsabs ( right.m_epsabs ) 00217 , m_epsrel ( right.m_epsrel ) 00218 , m_result ( GSL_NEGINF ) 00219 , m_error ( GSL_POSINF ) 00220 , m_size ( right.m_size ) 00221 , m_ws ( 0 ) 00222 , m_argument ( right.m_argument ) 00223 { 00224 m_pdata = new double[ 2 + m_points.size() ] ; // attention! 00225 }; 00226 // ======================================================================== 00227 00228 00229 // ======================================================================== 00231 // ======================================================================== 00232 NumericalIndefiniteIntegral::~NumericalIndefiniteIntegral() 00233 { 00234 if( 0 != m_ws ) 00235 { 00236 gsl_integration_workspace_free ( m_ws->ws ) ; 00237 delete m_ws ; 00238 m_ws = 0 ; 00239 } 00240 if ( 0 != m_pdata ) { delete m_pdata ; m_pdata = 0 ; } 00241 if ( 0 != m_function ) { delete m_function ; m_function = 0 ; } 00242 }; 00243 // ======================================================================== 00244 00245 // ======================================================================== 00246 // throw the exception 00247 // ======================================================================== 00248 StatusCode NumericalIndefiniteIntegral::Exception 00249 ( const std::string& message , 00250 const StatusCode& sc ) const 00251 { 00252 throw GaudiException( "NumericalIndefiniteIntegral::" + message , 00253 "*GaudiMath*" , sc ) ; 00254 return sc ; 00255 }; 00256 // ======================================================================== 00257 00258 // ======================================================================== 00260 // ======================================================================== 00261 double NumericalIndefiniteIntegral::operator() 00262 ( const double argument ) const 00263 { 00264 // reset the result and the error 00265 m_result = GSL_NEGINF ; 00266 m_error = GSL_POSINF ; 00267 00268 // check the argument 00269 if( 1 != m_DIM ) { Exception ( "operator(): invalid argument size " ) ; }; 00270 00271 m_argument[0] = argument ; 00272 return (*this) ( m_argument ); 00273 }; 00274 // ======================================================================== 00275 00276 // ======================================================================== 00278 // ======================================================================== 00279 Genfun::Derivative 00280 NumericalIndefiniteIntegral::partial ( unsigned int idx ) const 00281 { 00282 if ( idx >= m_DIM ) 00283 { Exception ( "::partial(i): invalid variable index " ) ; }; 00284 if ( idx != m_index ) 00285 { 00286 const AbsFunction& aux = NumericalDerivative( *this , idx ) ; 00287 return Genfun::FunctionNoop( &aux ) ; 00288 } 00289 else if ( GaudiMath::Integration::VariableLowLimit == limit () ) 00290 { 00291 const AbsFunction& aux = -1 * function() ; 00292 return Genfun::FunctionNoop( &aux ) ; 00293 } 00294 const AbsFunction& aux = function() ; 00295 return Genfun::FunctionNoop( &aux ) ; 00296 }; 00297 00298 // ======================================================================== 00300 // ======================================================================== 00301 double NumericalIndefiniteIntegral::operator() 00302 ( const Argument& argument ) const 00303 { 00304 // reset the result and the error 00305 m_result = GSL_NEGINF ; 00306 m_error = GSL_POSINF ; 00307 00308 // check the argument 00309 if( argument.dimension() != m_DIM ) 00310 { Exception ( "operator(): invalid argument size " ) ; }; 00311 00312 // copy the argument 00313 {for( size_t i = 0 ; i < m_DIM ; ++i ){ m_argument[i] = argument[i];}} 00314 00315 // create the helper object 00316 GSL_Helper helper( *m_function , m_argument , m_index ); 00317 00318 // use GSL to evaluate the numerical derivative 00319 gsl_function F ; 00320 F.function = &GSL_Adaptor ; 00321 F.params = &helper ; 00322 _Function F1 ; 00323 F1.fn = &F ; 00324 00325 if ( GaudiMath::Integration::Infinite == category () ) 00326 { return QAGI ( &F1 ) ; } // RETURN 00327 else if ( GaudiMath::Integration::Singular == category () ) 00328 { return QAGP ( &F1 ) ; } // RETURN 00329 else if ( GaudiMath::Integration::Finite == category () ) 00330 if ( GaudiMath::Integration::NonAdaptive == type () ) 00331 { return QNG ( &F1 ) ; } // RETURN 00332 else if ( GaudiMath::Integration::Adaptive == type () ) 00333 { return QAG ( &F1 ) ; } // RETURN 00334 else if ( GaudiMath::Integration::AdaptiveSingular == type () ) 00335 { return QAGS ( &F1 ) ; } // RETURN 00336 else 00337 { Exception ( "::operator(): invalid type " ); } 00338 else 00339 { Exception ( "::operator(): invalid category " ); } 00340 00341 return 0 ; 00342 }; 00343 // ======================================================================== 00344 00345 // ======================================================================== 00347 // ======================================================================== 00348 NumericalIndefiniteIntegral::_Workspace* 00349 NumericalIndefiniteIntegral::allocate() const 00350 { 00351 if ( 0 != m_ws ) { return m_ws; } 00352 gsl_integration_workspace* aux = 00353 gsl_integration_workspace_alloc( size () ); 00354 if ( 0 == aux ) { Exception ( "allocate()::invalid workspace" ) ; }; 00355 m_ws = new _Workspace() ; 00356 m_ws->ws = aux ; 00357 return m_ws ; 00358 }; 00359 // ======================================================================== 00360 00361 // ======================================================================== 00362 // adaptive integration on infinite interval 00363 // ======================================================================== 00364 double NumericalIndefiniteIntegral::QAGI ( _Function* F ) const 00365 { 00366 // check the argument 00367 if( 0 == F ) { Exception("::QAGI: invalid function"); } 00368 00369 const double x = m_argument[m_index] ; 00370 00371 // allocate workspace 00372 if( 0 == ws() ) { allocate() ; } 00373 00374 int ierror = 0 ; 00375 switch ( limit() ) 00376 { 00377 case GaudiMath::Integration::VariableLowLimit : 00378 ierror = gsl_integration_qagiu ( F->fn , x , 00379 m_epsabs , m_epsrel , 00380 size () , ws()->ws , 00381 &m_result , &m_error ) ; break ; 00382 case GaudiMath::Integration::VariableHighLimit : 00383 ierror = gsl_integration_qagil ( F->fn , x , 00384 m_epsabs , m_epsrel , 00385 size () , ws()->ws , 00386 &m_result , &m_error ) ; break ; 00387 default : 00388 Exception ( "::QAGI: invalid mode" ) ; 00389 }; 00390 00391 if( ierror ) { gsl_error( "NumericalIndefiniteIntegral::QAGI" , 00392 __FILE__ , __LINE__ , ierror ) ;} 00393 00394 return m_result ; 00395 }; 00396 // ======================================================================== 00397 00398 // ======================================================================== 00399 // adaptive integration with known singular points 00400 // ======================================================================== 00401 double NumericalIndefiniteIntegral::QAGP( _Function* F ) const 00402 { 00403 if( 0 == F ) { Exception("QAGP::invalid function") ; } 00404 00405 const double x = m_argument[m_index] ; 00406 if ( m_a == x ) 00407 { 00408 m_result = 0 ; 00409 m_error = 0 ; // EXACT ! 00410 return m_result ; 00411 } 00412 00413 // no known singular points ? 00414 if( points().empty() ) { return QAGS( F ) ; } 00415 00416 // integration limits 00417 const double a = std::min ( m_a , x ) ; 00418 const double b = std::max ( m_a , x ) ; 00419 00420 // "active" singular points 00421 Points::const_iterator lower = 00422 std::lower_bound ( points().begin() , points().end() , a ) ; 00423 Points::const_iterator upper = 00424 std::upper_bound ( points().begin() , points().end() , b ) ; 00425 00426 Points pnts ( upper - lower ) ; 00427 std::copy( lower , upper , pnts.begin() ); 00428 if ( *lower != a ) { pnts.insert( pnts.begin () , a ) ; } 00429 if ( *upper != b ) { pnts.insert( pnts.end () , b ) ; } 00430 std::copy( pnts.begin() , pnts.end() , m_pdata ); 00431 const size_t npts = pnts.size() ; 00432 00433 // use GSL 00434 int ierror = 00435 gsl_integration_qagp ( F->fn , 00436 m_pdata , npts , 00437 m_epsabs , m_epsrel , 00438 size () , ws()->ws , 00439 &m_result , &m_error ) ; 00440 00441 if( ierror ) { gsl_error( "NumericalIndefiniteIntegral::QAGI " , 00442 __FILE__ , __LINE__ , ierror ) ; } 00443 00444 // sign 00445 if ( GaudiMath::Integration::VariableHighLimit == limit() 00446 && x < m_a ) { m_result *= -1 ; } 00447 else if ( GaudiMath::Integration::VariableLowLimit == limit() 00448 && x > m_a ) { m_result *= -1 ; } 00449 00450 return m_result ; 00451 }; 00452 // ======================================================================== 00453 00454 // ======================================================================== 00455 // non-adaptive integration 00456 // ======================================================================== 00457 double NumericalIndefiniteIntegral::QNG ( _Function* F ) const 00458 { 00459 if( 0 == F ) { Exception("QNG::invalid function") ; } 00460 00461 const double x = m_argument[m_index] ; 00462 if ( m_a == x ) 00463 { 00464 m_result = 0 ; 00465 m_error = 0 ; // EXACT ! 00466 return m_result ; 00467 } 00468 00469 // integration limits 00470 const double a = std::min ( m_a , x ) ; 00471 const double b = std::max ( m_a , x ) ; 00472 00473 size_t neval = 0 ; 00474 int ierror = 00475 gsl_integration_qng ( F->fn , 00476 a , b , 00477 m_epsabs , m_epsrel , 00478 &m_result , &m_error , &neval ) ; 00479 00480 if( ierror ) { gsl_error( "NumericalIndefiniteIntegral::QNG " , 00481 __FILE__ , __LINE__ , ierror ) ; } 00482 00483 // sign 00484 if ( GaudiMath::Integration::VariableHighLimit == limit() 00485 && x < m_a ) { m_result *= -1 ; } 00486 else if ( GaudiMath::Integration::VariableLowLimit == limit() 00487 && x > m_a ) { m_result *= -1 ; } 00488 00489 return m_result ; 00490 }; 00491 00492 // ======================================================================== 00493 // simple adaptive integration 00494 // ======================================================================== 00495 double NumericalIndefiniteIntegral::QAG ( _Function* F ) const 00496 { 00497 if( 0 == F ) { Exception("QAG::invalid function") ; } 00498 00499 const double x = m_argument[m_index] ; 00500 if ( m_a == x ) 00501 { 00502 m_result = 0 ; 00503 m_error = 0 ; // EXACT ! 00504 return m_result ; 00505 } 00506 00507 // allocate workspace 00508 if( 0 == ws () ) { allocate () ; } 00509 00510 // integration limits 00511 const double a = std::min ( m_a , x ) ; 00512 const double b = std::max ( m_a , x ) ; 00513 00514 int ierror = 00515 gsl_integration_qag ( F->fn , 00516 a , b , 00517 m_epsabs , m_epsrel , 00518 size () , (int) rule() , ws ()->ws , 00519 &m_result , &m_error ); 00520 00521 if( ierror ) { gsl_error( "NumericalIndefiniteIntegral::QAG " , 00522 __FILE__ , __LINE__ , ierror ) ; } 00523 00524 // sign 00525 if ( GaudiMath::Integration::VariableHighLimit == limit() 00526 && x < m_a ) { m_result *= -1 ; } 00527 else if ( GaudiMath::Integration::VariableLowLimit == limit() 00528 && x > m_a ) { m_result *= -1 ; } 00529 00530 return m_result ; 00531 }; 00532 // ======================================================================== 00533 00534 // ======================================================================== 00535 // adaptive integration with singularities 00536 // ======================================================================== 00537 double NumericalIndefiniteIntegral::QAGS ( _Function* F ) const 00538 { 00539 if( 0 == F ) { Exception("QAG::invalid function") ; } 00540 00541 const double x = m_argument[m_index] ; 00542 if ( m_a == x ) 00543 { 00544 m_result = 0 ; 00545 m_error = 0 ; // EXACT ! 00546 return m_result ; 00547 } 00548 00549 // allocate workspace 00550 if( 0 == ws () ) { allocate () ; } 00551 00552 // integration limits 00553 const double a = std::min ( m_a , x ) ; 00554 const double b = std::max ( m_a , x ) ; 00555 00556 int ierror = 00557 gsl_integration_qags ( F->fn , 00558 a , b , 00559 m_epsabs , m_epsrel , 00560 size () , ws()->ws , 00561 &m_result , &m_error ); 00562 00563 if( ierror ) { gsl_error( "NumericalIndefiniteIntegral::QAGS " , 00564 __FILE__ , __LINE__ , ierror ) ; } 00565 00566 // sign 00567 if ( GaudiMath::Integration::VariableHighLimit == limit() 00568 && x < m_a ) { m_result *= -1 ; } 00569 else if ( GaudiMath::Integration::VariableLowLimit == limit() 00570 && x > m_a ) { m_result *= -1 ; } 00571 00572 return m_result ; 00573 }; 00574 // ======================================================================== 00575 00576 }; // end of namespace GaudiMathImplementation 00577 }; // end of namespace Genfun 00578 00579 00580 // ============================================================================ 00581 // The END 00582 // ============================================================================