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