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