Gaudi Framework, version v20r4
Generated: 8 Jan 2009
|
|
|
Generated: 8 Jan 2009 |
00001 // $Id: NumericalDefiniteIntegral.cpp,v 1.4 2007/11/20 13:00:17 marcocle Exp $ 00002 // ============================================================================ 00003 // Include 00004 // ============================================================================ 00005 // STD & STL 00006 // ============================================================================ 00007 #include <vector> 00008 #include <algorithm> 00009 // ============================================================================ 00010 // GSL 00011 // ============================================================================ 00012 #include "gsl/gsl_errno.h" 00013 #include "gsl/gsl_integration.h" 00014 // ============================================================================ 00015 // GaudiKernel 00016 // ============================================================================ 00017 #include "GaudiKernel/GaudiException.h" 00018 // ============================================================================ 00019 // GaudiMath 00020 // ============================================================================ 00021 #include "GaudiMath/NumericalDefiniteIntegral.h" 00022 #include "GaudiMath/NumericalDerivative.h" 00023 // ============================================================================ 00024 // Local 00025 // ============================================================================ 00026 #include "Helpers.h" 00027 // ============================================================================ 00028 00029 00030 00031 namespace Genfun 00032 { 00033 namespace GaudiMathImplementation 00034 { 00035 00036 struct NumericalDefiniteIntegral::_Workspace 00037 { gsl_integration_workspace* ws ; }; 00038 struct NumericalDefiniteIntegral::_Function 00039 { gsl_function* fn ; }; 00040 00041 // ======================================================================== 00043 // ======================================================================== 00044 FUNCTION_OBJECT_IMP( NumericalDefiniteIntegral ); 00045 // ======================================================================== 00046 00098 NumericalDefiniteIntegral::NumericalDefiniteIntegral 00099 ( const AbsFunction& function , 00100 const size_t index , 00101 const double a , 00102 const double b , 00103 const GaudiMath::Integration::Type type , 00104 const GaudiMath::Integration::KronrodRule rule , 00105 const double epsabs , 00106 const double epsrel , 00107 const size_t size ) 00108 : AbsFunction () 00109 , m_function ( function.clone () ) 00110 , m_DIM ( 0 ) 00111 , m_index ( index ) 00112 , m_a ( a ) 00113 , m_b ( b ) 00114 , m_ia ( false ) 00115 , m_ib ( false ) 00116 , m_type ( type ) 00117 , m_category ( GaudiMath::Integration::Finite ) 00118 , m_rule ( rule ) 00119 , m_points ( ) 00120 , m_pdata ( 0 ) 00121 , m_epsabs ( epsabs ) 00122 , m_epsrel ( epsrel ) 00123 , m_result ( GSL_NEGINF ) 00124 , m_error ( GSL_POSINF ) 00125 , m_size ( size ) 00126 , m_ws () 00127 , m_argument () 00128 , m_argF () 00129 { 00130 if ( GaudiMath::Integration::Fixed == m_rule ) 00131 { m_rule = GaudiMath::Integration::Default ; } 00132 if ( function.dimensionality() < 2 ) 00133 { Exception("::constructor: invalid dimensionality ") ; } 00134 if ( m_index >= function.dimensionality() ) 00135 { Exception("::constructor: invalid variable index") ; } 00136 00137 m_DIM = function.dimensionality() - 1 ; 00138 m_argument = Argument( m_DIM ) ; 00139 m_argF = Argument( m_DIM + 1 ) ; 00140 00141 }; 00142 00143 00153 NumericalDefiniteIntegral::NumericalDefiniteIntegral 00154 ( const AbsFunction& function , 00155 const size_t index , 00156 const double a , 00157 const double b , 00158 const NumericalDefiniteIntegral::Points& points , 00159 const double epsabs , 00160 const double epsrel , 00161 const size_t size ) 00162 : AbsFunction () 00163 , m_function ( function.clone() ) 00164 , m_DIM ( 0 ) 00165 , m_index ( index ) 00166 , m_a ( a ) 00167 , m_b ( b ) 00168 , m_ia ( false ) 00169 , m_ib ( false ) 00170 , m_type ( GaudiMath::Integration:: Other ) 00171 , m_category ( GaudiMath::Integration:: Singular ) 00172 , m_rule ( GaudiMath::Integration:: Fixed ) 00173 , m_points ( points ) 00174 , m_pdata ( 0 ) 00175 , m_epsabs ( epsabs ) 00176 , m_epsrel ( epsrel ) 00177 // 00178 , m_result ( GSL_NEGINF ) 00179 , m_error ( GSL_POSINF ) 00180 // 00181 , m_size ( size ) 00182 , m_ws ( 0 ) 00183 , m_argument () 00184 , m_argF () 00185 { 00186 if ( function.dimensionality() < 2 ) 00187 { Exception("::constructor: invalid dimensionality ") ; } 00188 if ( m_index >= function.dimensionality() ) 00189 { Exception("::constructor: invalid variable index") ; } 00190 00191 m_DIM = function.dimensionality() - 1 ; 00192 m_argument = Argument( m_DIM ) ; 00193 m_argF = Argument( m_DIM + 1 ) ; 00194 00195 const double l1 = std::min ( a , b ) ; 00196 const double l2 = std::max ( a , b ) ; 00197 m_points.push_back ( l1 ) ; 00198 m_points.push_back ( l2 ) ; 00199 std::sort ( m_points.begin() , m_points.end() ) ; 00200 m_points.erase( std::unique( m_points.begin () , 00201 m_points.end () ) , 00202 m_points.end() ); 00203 00204 Points::iterator lower = 00205 std::lower_bound ( m_points.begin () , m_points.end () , l1 ) ; 00206 m_points.erase ( m_points.begin () , lower ) ; 00207 Points::iterator upper = 00208 std::upper_bound ( m_points.begin () , m_points.end () , l2 ) ; 00209 m_points.erase ( upper , m_points.end () ) ; 00210 00211 m_pdata = new double[ m_points.size() ] ; 00212 std::copy( m_points.begin() , m_points.end() , m_pdata ); 00213 }; 00214 00240 NumericalDefiniteIntegral::NumericalDefiniteIntegral 00241 ( const AbsFunction& function , 00242 const size_t index , 00243 const double a , 00244 const GaudiMath::Integration::Inf /* b */ , 00245 const double epsabs , 00246 const double epsrel , 00247 const size_t size ) 00248 : AbsFunction() 00249 , m_function ( function.clone() ) 00250 , m_DIM ( 0 ) 00251 , m_index ( index ) 00252 , m_a ( a ) 00253 , m_b ( GSL_POSINF ) 00254 , m_ia ( false ) 00255 , m_ib ( true ) 00256 , m_type ( GaudiMath::Integration:: Other ) 00257 , m_category ( GaudiMath::Integration:: Infinite ) 00258 , m_rule ( GaudiMath::Integration:: Fixed ) 00259 , m_points ( ) 00260 , m_pdata ( 0 ) 00261 , m_epsabs ( epsabs ) 00262 , m_epsrel ( epsrel ) 00263 // 00264 , m_result ( GSL_NEGINF ) 00265 , m_error ( GSL_POSINF ) 00266 // 00267 , m_size ( size ) 00268 , m_ws ( 0 ) 00269 , m_argument () 00270 , m_argF () 00271 { 00272 if ( function.dimensionality() < 2 ) 00273 { Exception("::constructor: invalid dimensionality ") ; } 00274 if ( m_index >= function.dimensionality() ) 00275 { Exception("::constructor: invalid variable index") ; } 00276 00277 m_DIM = function.dimensionality() - 1 ; 00278 m_argument = Argument( m_DIM ) ; 00279 m_argF = Argument( m_DIM + 1 ) ; 00280 00281 }; 00282 00283 00309 NumericalDefiniteIntegral::NumericalDefiniteIntegral 00310 ( const AbsFunction& function , 00311 const size_t index , 00312 const GaudiMath::Integration::Inf /* a */ , 00313 const double b , 00314 const double epsabs , 00315 const double epsrel , 00316 const size_t size ) 00317 : AbsFunction() 00318 , m_function ( function.clone() ) 00319 , m_DIM ( 0 ) 00320 , m_index ( index ) 00321 , m_a ( GSL_NEGINF ) 00322 , m_b ( b ) 00323 , m_ia ( true ) 00324 , m_ib ( false ) 00325 , m_type ( GaudiMath::Integration:: Other ) 00326 , m_category ( GaudiMath::Integration:: Infinite ) 00327 , m_rule ( GaudiMath::Integration:: Fixed ) 00328 , m_points ( ) 00329 , m_pdata ( 0 ) 00330 , m_epsabs ( epsabs ) 00331 , m_epsrel ( epsrel ) 00332 // 00333 , m_result ( GSL_NEGINF ) 00334 , m_error ( GSL_POSINF ) 00335 // 00336 , m_size ( size ) 00337 , m_ws ( 0 ) 00338 , m_argument () 00339 , m_argF () 00340 { 00341 if ( function.dimensionality() < 2 ) 00342 { Exception("::constructor: invalid dimensionality ") ; } 00343 if ( m_index >= function.dimensionality() ) 00344 { Exception("::constructor: invalid variable index") ; } 00345 00346 m_DIM = function.dimensionality() - 1 ; 00347 m_argument = Argument( m_DIM ) ; 00348 m_argF = Argument( m_DIM + 1 ) ; 00349 }; 00350 00374 NumericalDefiniteIntegral::NumericalDefiniteIntegral 00375 ( const AbsFunction& function , 00376 const size_t index , 00377 const float epsabs , 00378 const float epsrel , 00379 const size_t size ) 00380 : AbsFunction() 00381 , m_function ( function.clone() ) 00382 , m_DIM ( 0 ) 00383 , m_index ( index ) 00384 , m_a ( GSL_NEGINF ) 00385 , m_b ( GSL_POSINF ) 00386 , m_ia ( true ) 00387 , m_ib ( true ) 00388 , m_type ( GaudiMath::Integration:: Other ) 00389 , m_category ( GaudiMath::Integration:: Infinite ) 00390 , m_rule ( GaudiMath::Integration:: Fixed ) 00391 , m_points ( ) 00392 , m_pdata ( 0 ) 00393 , m_epsabs ( epsabs ) 00394 , m_epsrel ( epsrel ) 00395 // 00396 , m_result ( GSL_NEGINF ) 00397 , m_error ( GSL_POSINF ) 00398 // 00399 , m_size ( size ) 00400 , m_ws ( 0 ) 00401 , m_argument () 00402 , m_argF () 00403 { 00404 if ( function.dimensionality() < 2 ) 00405 { Exception("::constructor: invalid dimensionality ") ; } 00406 if ( m_index >= function.dimensionality() ) 00407 { Exception("::constructor: invalid variable index") ; } 00408 00409 m_DIM = function.dimensionality() - 1 ; 00410 m_argument = Argument( m_DIM ) ; 00411 m_argF = Argument( m_DIM + 1 ) ; 00412 00413 }; 00414 00416 NumericalDefiniteIntegral::NumericalDefiniteIntegral 00417 ( const NumericalDefiniteIntegral& right ) 00418 : AbsFunction () 00419 , m_function ( right.m_function->clone() ) 00420 , m_DIM ( right.m_DIM ) 00421 , m_index ( right.m_index ) 00422 , m_a ( right.m_a ) 00423 , m_b ( right.m_b ) 00424 , m_ia ( right.m_ia ) 00425 , m_ib ( right.m_ib ) 00426 , m_type ( right.m_type ) 00427 , m_category ( right.m_category ) 00428 , m_rule ( right.m_rule ) 00429 , m_points ( right.m_points ) 00430 , m_pdata ( 0 ) 00431 , m_epsabs ( right.m_epsabs ) 00432 , m_epsrel ( right.m_epsrel ) 00433 , m_result ( GSL_NEGINF ) 00434 , m_error ( GSL_POSINF ) 00435 , m_size ( right.m_size ) 00436 , m_ws ( 0 ) 00437 , m_argument ( right.m_argument ) 00438 , m_argF ( right.m_argF ) 00439 { 00440 if( 0 != right.m_pdata ) 00441 { 00442 m_pdata = new double[m_points.size()] ; 00443 std::copy( m_points.begin() , m_points.end() , m_pdata ); 00444 } 00445 }; 00446 00447 NumericalDefiniteIntegral::~NumericalDefiniteIntegral() 00448 { 00449 if( 0 != m_function ) { delete m_function ; m_function = 0 ; } 00450 if( 0 != m_pdata ) { delete m_pdata ; m_pdata = 0 ; } 00451 }; 00452 00453 // ======================================================================== 00454 // throw the exception 00455 // ======================================================================== 00456 StatusCode NumericalDefiniteIntegral::Exception 00457 ( const std::string& message , 00458 const StatusCode& sc ) const 00459 { 00460 throw GaudiException( "NumericalDefiniteIntegral::" + message , 00461 "*GaudiMath*" , sc ) ; 00462 return sc ; 00463 }; 00464 // ======================================================================== 00465 00466 // ======================================================================== 00468 // ======================================================================== 00469 double NumericalDefiniteIntegral::operator() 00470 ( const double argument ) const 00471 { 00472 // reset the result and the error 00473 m_result = GSL_NEGINF ; 00474 m_error = GSL_POSINF ; 00475 00476 // check the argument 00477 if( 1 != m_DIM ) { Exception ( "operator(): invalid argument size " ) ; }; 00478 00479 m_argument[0] = argument ; 00480 return (*this) ( m_argument ); 00481 }; 00482 // ======================================================================== 00483 00484 // ======================================================================== 00486 // ======================================================================== 00487 Genfun::Derivative 00488 NumericalDefiniteIntegral::partial ( unsigned int idx ) const 00489 { 00490 if ( idx >= m_DIM ) 00491 { Exception ( "::partial(i): invalid variable index " ) ; }; 00492 // 00493 const AbsFunction& aux = NumericalDerivative( *this , idx ) ; 00494 return Genfun::FunctionNoop( &aux ) ; 00495 }; 00496 // ======================================================================== 00497 00498 00499 // ======================================================================== 00501 // ======================================================================== 00502 double NumericalDefiniteIntegral::operator() 00503 ( const Argument& argument ) const 00504 { 00505 // reset the result and the error 00506 m_result = GSL_NEGINF ; 00507 m_error = GSL_POSINF ; 00508 00509 // check the argument 00510 if( argument.dimension() != m_DIM ) 00511 { Exception ( "operator(): invalid argument size " ) ; }; 00512 00513 // copy the argument 00514 00515 for( size_t i = 0 ; i < m_DIM ; ++i ) 00516 { 00517 m_argument [i] = argument [i] ; 00518 const size_t j = i < m_index ? i : i + 1 ; 00519 m_argF [j] = argument [i] ; 00520 }; 00521 00522 // create the helper object 00523 GSL_Helper helper( *m_function , m_argF , m_index ); 00524 00525 // use GSL to evaluate the numerical derivative 00526 gsl_function F ; 00527 F.function = &GSL_Adaptor ; 00528 F.params = &helper ; 00529 _Function F1 ; 00530 F1.fn = &F ; 00531 00532 if ( GaudiMath::Integration::Infinite == category () ) 00533 { return QAGI ( &F1 ) ; } // RETURN 00534 00535 if ( m_a == m_b ) 00536 { 00537 m_result = 0 ; m_error = 0 ; // EXACT 00538 return m_result ; // RETURN 00539 } 00540 00541 if ( GaudiMath::Integration::Singular == category () ) 00542 { return QAGP ( &F1 ) ; } // RETURN 00543 else if ( GaudiMath::Integration::Finite == category () ) 00544 if ( GaudiMath::Integration::NonAdaptive == type () ) 00545 { return QNG ( &F1 ) ; } // RETURN 00546 else if ( GaudiMath::Integration::Adaptive == type () ) 00547 { return QAG ( &F1 ) ; } // RETURN 00548 else if ( GaudiMath::Integration::AdaptiveSingular == type () ) 00549 { return QAGS ( &F1 ) ; } // RETURN 00550 else 00551 { Exception ( "::operator(): invalid type " ); } 00552 else 00553 { Exception ( "::operator(): invalid category " ); } 00554 00555 return 0 ; 00556 }; 00557 // ======================================================================== 00558 00559 // ======================================================================== 00561 // ======================================================================== 00562 NumericalDefiniteIntegral::_Workspace* 00563 NumericalDefiniteIntegral::allocate() const 00564 { 00565 if ( 0 != m_ws ) { return m_ws; } 00566 gsl_integration_workspace* aux = 00567 gsl_integration_workspace_alloc( size () ); 00568 if ( 0 == aux ) { Exception ( "allocate()::invalid workspace" ) ; }; 00569 m_ws = new _Workspace() ; 00570 m_ws->ws = aux ; 00571 return m_ws ; 00572 }; 00573 // ======================================================================== 00574 00575 // ======================================================================== 00576 // adaptive integration on infinite interval 00577 // ======================================================================== 00578 double NumericalDefiniteIntegral::QAGI ( _Function* F ) const 00579 { 00580 // check the argument 00581 if ( 0 == F ) { Exception("::QAGI: invalid function"); } 00582 00583 // allocate workspace 00584 if ( 0 == ws() ) { allocate() ; } 00585 00586 int ierror = 0 ; 00587 00588 if ( m_ia && m_ib ) 00589 { 00590 ierror = gsl_integration_qagi ( F->fn , 00591 m_epsabs , m_epsrel , 00592 size () , ws()->ws , 00593 &m_result , &m_error ) ; 00594 } 00595 else if ( m_ia ) 00596 { 00597 ierror = gsl_integration_qagil ( F->fn , m_b , 00598 m_epsabs , m_epsrel , 00599 size () , ws()->ws , 00600 &m_result , &m_error ) ; 00601 } 00602 else if ( m_ib ) 00603 { 00604 ierror = gsl_integration_qagiu ( F->fn , m_a , 00605 m_epsabs , m_epsrel , 00606 size () , ws()->ws , 00607 &m_result , &m_error ) ; 00608 } 00609 else 00610 { Exception ( "::QAGI: invalid mode" ) ; }; 00611 00612 if( ierror ) { gsl_error( "NumericalDefiniteIntegral::QAGI" , 00613 __FILE__ , __LINE__ , ierror ) ;} 00614 00615 return m_result ; 00616 }; 00617 // ======================================================================== 00618 00619 // ======================================================================== 00620 // adaptive integration with known singular points 00621 // ======================================================================== 00622 double NumericalDefiniteIntegral::QAGP( _Function* F ) const 00623 { 00624 if( 0 == F ) { Exception("QAGP::invalid function") ; } 00625 00626 // no known singular points ? 00627 if( points().empty() || 0 == m_pdata ) { return QAGS( F ) ; } 00628 00629 const size_t npts = points().size(); 00630 00631 // use GSL 00632 int ierror = 00633 gsl_integration_qagp ( F->fn , 00634 m_pdata , npts , 00635 m_epsabs , m_epsrel , 00636 size () , ws()->ws , 00637 &m_result , &m_error ) ; 00638 00639 if( ierror ) { gsl_error( "NumericalDefiniteIntegral::QAGI " , 00640 __FILE__ , __LINE__ , ierror ) ; } 00641 00642 // the sign 00643 if ( m_a > m_b ) { m_result *= -1 ; } 00644 00645 return m_result ; 00646 }; 00647 // ======================================================================== 00648 00649 // ======================================================================== 00650 // non-adaptive integration 00651 // ======================================================================== 00652 double NumericalDefiniteIntegral::QNG ( _Function* F ) const 00653 { 00654 if( 0 == F ) { Exception("QNG::invalid function") ; } 00655 00656 // integration limits 00657 const double low = std::min ( m_a , m_b ) ; 00658 const double high = std::max ( m_a , m_b ) ; 00659 00660 size_t neval = 0 ; 00661 int ierror = 00662 gsl_integration_qng ( F->fn , 00663 low , high , 00664 m_epsabs , m_epsrel , 00665 &m_result , &m_error , &neval ) ; 00666 00667 if( ierror ) { gsl_error( "NumericalIndefiniteIntegral::QNG " , 00668 __FILE__ , __LINE__ , ierror ) ; } 00669 00670 // sign 00671 if ( m_a > m_b ) { m_result *= -1 ; } 00672 00673 return m_result ; 00674 }; 00675 // ======================================================================== 00676 00677 00678 // ======================================================================== 00679 // simple adaptive integration 00680 // ======================================================================== 00681 double NumericalDefiniteIntegral::QAG ( _Function* F ) const 00682 { 00683 if( 0 == F ) { Exception("QAG::invalid function") ; } 00684 00685 // allocate workspace 00686 if( 0 == ws () ) { allocate () ; } 00687 00688 // integration limits 00689 const double low = std::min ( m_a , m_b ) ; 00690 const double high = std::max ( m_a , m_b ) ; 00691 00692 int ierror = 00693 gsl_integration_qag ( F->fn , 00694 low , high , 00695 m_epsabs , m_epsrel , 00696 size () , (int) rule() , ws ()->ws , 00697 &m_result , &m_error ); 00698 00699 if( ierror ) { gsl_error( "NumericalDefiniteIntegral::QAG " , 00700 __FILE__ , __LINE__ , ierror ) ; } 00701 00702 // sign 00703 if ( m_a > m_b ) { m_result *= -1 ; } 00704 00705 return m_result ; 00706 }; 00707 // ======================================================================== 00708 00709 // ======================================================================== 00710 // adaptive integration with singularities 00711 // ======================================================================== 00712 double NumericalDefiniteIntegral::QAGS ( _Function* F ) const 00713 { 00714 if( 0 == F ) { Exception("QAG::invalid function") ; } 00715 00716 if ( m_a == m_b ) 00717 { 00718 m_result = 0 ; 00719 m_error = 0 ; // EXACT ! 00720 return m_result ; 00721 } 00722 00723 // allocate workspace 00724 if( 0 == ws () ) { allocate () ; } 00725 00726 // integration limits 00727 const double low = std::min ( m_a , m_b ) ; 00728 const double high = std::max ( m_a , m_b ) ; 00729 00730 int ierror = 00731 gsl_integration_qags ( F->fn , 00732 low , high , 00733 m_epsabs , m_epsrel , 00734 size () , ws()->ws , 00735 &m_result , &m_error ); 00736 00737 if( ierror ) { gsl_error( "NumericalDefiniteIntegral::QAGS " , 00738 __FILE__ , __LINE__ , ierror ) ; } 00739 00740 // sign 00741 if ( m_a > m_b ) { m_result *= -1 ; } 00742 00743 return m_result ; 00744 }; 00745 // ======================================================================== 00746 00747 00748 }; // end of namespace GaudiMathImplementation 00749 }; // end of namespace Genfun 00750 00751 00752 00753 // ============================================================================ 00754 // The END 00755 // ============================================================================