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
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
00047
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 ,
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 ,
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
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
00473 m_result = GSL_NEGINF ;
00474 m_error = GSL_POSINF ;
00475
00476
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
00506 m_result = GSL_NEGINF ;
00507 m_error = GSL_POSINF ;
00508
00509
00510 if( argument.dimension() != m_DIM )
00511 { Exception ( "operator(): invalid argument size " ) ; };
00512
00513
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
00523 GSL_Helper helper( *m_function , m_argF , m_index );
00524
00525
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 ) ; }
00534
00535 if ( m_a == m_b )
00536 {
00537 m_result = 0 ; m_error = 0 ;
00538 return m_result ;
00539 }
00540
00541 if ( GaudiMath::Integration::Singular == category () )
00542 { return QAGP ( &F1 ) ; }
00543 else if ( GaudiMath::Integration::Finite == category () )
00544 if ( GaudiMath::Integration::NonAdaptive == type () )
00545 { return QNG ( &F1 ) ; }
00546 else if ( GaudiMath::Integration::Adaptive == type () )
00547 { return QAG ( &F1 ) ; }
00548 else if ( GaudiMath::Integration::AdaptiveSingular == type () )
00549 { return QAGS ( &F1 ) ; }
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
00577
00578 double NumericalDefiniteIntegral::QAGI ( _Function* F ) const
00579 {
00580
00581 if ( 0 == F ) { Exception("::QAGI: invalid function"); }
00582
00583
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
00621
00622 double NumericalDefiniteIntegral::QAGP( _Function* F ) const
00623 {
00624 if( 0 == F ) { Exception("QAGP::invalid function") ; }
00625
00626
00627 if( points().empty() || 0 == m_pdata ) { return QAGS( F ) ; }
00628
00629 const size_t npts = points().size();
00630
00631
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
00643 if ( m_a > m_b ) { m_result *= -1 ; }
00644
00645 return m_result ;
00646 }
00647
00648
00649
00650
00651
00652 double NumericalDefiniteIntegral::QNG ( _Function* F ) const
00653 {
00654 if( 0 == F ) { Exception("QNG::invalid function") ; }
00655
00656
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
00671 if ( m_a > m_b ) { m_result *= -1 ; }
00672
00673 return m_result ;
00674 }
00675
00676
00677
00678
00679
00680
00681 double NumericalDefiniteIntegral::QAG ( _Function* F ) const
00682 {
00683 if( 0 == F ) { Exception("QAG::invalid function") ; }
00684
00685
00686 if( 0 == ws () ) { allocate () ; }
00687
00688
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
00703 if ( m_a > m_b ) { m_result *= -1 ; }
00704
00705 return m_result ;
00706 }
00707
00708
00709
00710
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 ;
00720 return m_result ;
00721 }
00722
00723
00724 if( 0 == ws () ) { allocate () ; }
00725
00726
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
00741 if ( m_a > m_b ) { m_result *= -1 ; }
00742
00743 return m_result ;
00744 }
00745
00746
00747
00748 }
00749 }
00750
00751
00752
00753
00754
00755
