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