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