Loading [MathJax]/extensions/tex2jax.js
The Gaudi Framework  v31r0 (aeb156f0)
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
NumericalDefiniteIntegral.cpp
Go to the documentation of this file.
1 // ============================================================================
2 // Include
3 // ============================================================================
4 // STD & STL
5 // ============================================================================
6 #include <algorithm>
7 #include <vector>
8 // ============================================================================
9 // GSL
10 // ============================================================================
11 #include "gsl/gsl_errno.h"
12 #include "gsl/gsl_integration.h"
13 // ============================================================================
14 // GaudiKernel
15 // ============================================================================
16 #include "GaudiKernel/GaudiException.h"
17 // ============================================================================
18 // GaudiMath
19 // ============================================================================
20 #include "GaudiMath/NumericalDefiniteIntegral.h"
21 #include "GaudiMath/NumericalDerivative.h"
22 // ============================================================================
23 // Local
24 // ============================================================================
25 #include "Helpers.h"
26 // ============================================================================
27 
28 #ifdef __ICC
29 // disable icc remark #1572: floating-point equality and inequality comparisons are unreliable
30 // The comparison are meant
31 # pragma warning( disable : 1572 )
32 #endif
33 #ifdef WIN32
34 // Disable the warning
35 // C4996: 'std::copy': Function call with parameters that may be unsafe
36 // The parameters are checked
37 # pragma warning( disable : 4996 )
38 #endif
39 
40 namespace Genfun {
41  namespace GaudiMathImplementation {
42 
43  struct NumericalDefiniteIntegral::_Function {
44  gsl_function* fn;
45  };
46 
47  // ========================================================================
48 
100  NumericalDefiniteIntegral::NumericalDefiniteIntegral( const AbsFunction& function, const size_t index,
101  const double a, const double b,
102  const GaudiMath::Integration::Type type,
103  const GaudiMath::Integration::KronrodRule rule,
104  const double epsabs, const double epsrel, const size_t size )
105  : m_function( function.clone() )
106  , m_index( index )
107  , m_a( a )
108  , m_b( b )
109  , m_type( type )
110  , m_category( GaudiMath::Integration::Finite )
111  , m_rule( rule )
112  , m_epsabs( epsabs )
113  , m_epsrel( epsrel )
114  , m_size( size ) {
115  if ( GaudiMath::Integration::Fixed == m_rule ) { m_rule = GaudiMath::Integration::Default; }
116  if ( function.dimensionality() < 2 ) { Exception( "::constructor: invalid dimensionality " ); }
117  if ( m_index >= function.dimensionality() ) { Exception( "::constructor: invalid variable index" ); }
118 
119  m_DIM = function.dimensionality() - 1;
120  m_argument = Argument( m_DIM );
121  m_argF = Argument( m_DIM + 1 );
122  }
123 
133  NumericalDefiniteIntegral::NumericalDefiniteIntegral( const AbsFunction& function, const size_t index,
134  const double a, const double b,
135  const NumericalDefiniteIntegral::Points& points,
136  const double epsabs, const double epsrel, const size_t size )
137  : m_function( function.clone() )
138  , m_index( index )
139  , m_a( a )
140  , m_b( b )
141  , m_type( GaudiMath::Integration::Other )
142  , m_category( GaudiMath::Integration::Singular )
143  , m_rule( GaudiMath::Integration::Fixed )
144  , m_points( points )
145  , m_epsabs( epsabs )
146  , m_epsrel( epsrel )
147  , m_size( size ) {
148  if ( function.dimensionality() < 2 ) { Exception( "::constructor: invalid dimensionality " ); }
149  if ( m_index >= function.dimensionality() ) { Exception( "::constructor: invalid variable index" ); }
150 
151  m_DIM = function.dimensionality() - 1;
152  m_argument = Argument( m_DIM );
153  m_argF = Argument( m_DIM + 1 );
154 
155  const auto lim = std::minmax( a, b );
156  m_points.push_back( lim.first );
157  m_points.push_back( lim.second );
158  auto end = std::remove_if( m_points.begin(), m_points.end(),
159  [&]( double x ) { return x < lim.first || lim.second < x; } );
160  std::sort( m_points.begin(), end );
161  m_points.erase( std::unique( m_points.begin(), end ), m_points.end() );
162  }
163 
189  NumericalDefiniteIntegral::NumericalDefiniteIntegral( const AbsFunction& function, const size_t index,
190  const double a, const GaudiMath::Integration::Inf /* b */,
191  const double epsabs, const double epsrel, const size_t size )
192  : m_function( function.clone() )
193  , m_DIM( 0 )
194  , m_index( index )
195  , m_a( a )
196  , m_type( GaudiMath::Integration::Other )
197  , m_category( GaudiMath::Integration::Infinite )
198  , m_rule( GaudiMath::Integration::Fixed )
199  , m_epsabs( epsabs )
200  , m_epsrel( epsrel )
201  , m_size( size ) {
202  if ( function.dimensionality() < 2 ) { Exception( "::constructor: invalid dimensionality " ); }
203  if ( m_index >= function.dimensionality() ) { Exception( "::constructor: invalid variable index" ); }
204 
205  m_DIM = function.dimensionality() - 1;
206  m_argument = Argument( m_DIM );
207  m_argF = Argument( m_DIM + 1 );
208  }
209 
235  NumericalDefiniteIntegral::NumericalDefiniteIntegral( const AbsFunction& function, const size_t index,
236  const GaudiMath::Integration::Inf /* a */, const double b,
237  const double epsabs, const double epsrel, const size_t size )
238  : m_function( function.clone() )
239  , m_DIM( 0 )
240  , m_index( index )
241  , m_b( b )
242  , m_type( GaudiMath::Integration::Other )
243  , m_category( GaudiMath::Integration::Infinite )
244  , m_rule( GaudiMath::Integration::Fixed )
245  , m_epsabs( epsabs )
246  , m_epsrel( epsrel )
247  //
248  , m_size( size ) {
249  if ( function.dimensionality() < 2 ) { Exception( "::constructor: invalid dimensionality " ); }
250  if ( m_index >= function.dimensionality() ) { Exception( "::constructor: invalid variable index" ); }
251 
252  m_DIM = function.dimensionality() - 1;
253  m_argument = Argument( m_DIM );
254  m_argF = Argument( m_DIM + 1 );
255  }
256 
280  NumericalDefiniteIntegral::NumericalDefiniteIntegral( const AbsFunction& function, const size_t index,
281  const float epsabs, const float epsrel, const size_t size )
282  : m_function( function.clone() )
283  , m_DIM( 0 )
284  , m_index( index )
285  , m_type( GaudiMath::Integration::Other )
286  , m_category( GaudiMath::Integration::Infinite )
287  , m_rule( GaudiMath::Integration::Fixed )
288  , m_epsabs( epsabs )
289  , m_epsrel( epsrel )
290  , m_size( size ) {
291  if ( function.dimensionality() < 2 ) { Exception( "::constructor: invalid dimensionality " ); }
292  if ( m_index >= function.dimensionality() ) { Exception( "::constructor: invalid variable index" ); }
293 
294  m_DIM = function.dimensionality() - 1;
295  m_argument = Argument( m_DIM );
296  m_argF = Argument( m_DIM + 1 );
297  }
298 
300  NumericalDefiniteIntegral::NumericalDefiniteIntegral( const NumericalDefiniteIntegral& right )
301  : AbsFunction()
302  , m_function( right.m_function->clone() )
303  , m_DIM( right.m_DIM )
304  , m_index( right.m_index )
305  , m_a( right.m_a )
306  , m_b( right.m_b )
307  , m_type( right.m_type )
308  , m_category( right.m_category )
309  , m_rule( right.m_rule )
310  , m_points( right.m_points )
311  , m_epsabs( right.m_epsabs )
312  , m_epsrel( right.m_epsrel )
313  , m_size( right.m_size )
314  , m_argument( right.m_argument )
315  , m_argF( right.m_argF ) {}
316 
317  // ========================================================================
318  // throw the exception
319  // ========================================================================
320  StatusCode NumericalDefiniteIntegral::Exception( const std::string& message, const StatusCode& sc ) const {
321  throw GaudiException( "NumericalDefiniteIntegral::" + message, "*GaudiMath*", sc );
322  return sc;
323  }
324  // ========================================================================
325 
326  // ========================================================================
328  // ========================================================================
329  double NumericalDefiniteIntegral::operator()( const double argument ) const {
330  // reset the result and the error
331  m_result = -std::numeric_limits<double>::infinity();
332  m_error = std::numeric_limits<double>::infinity();
333 
334  // check the argument
335  if ( 1 != m_DIM ) { Exception( "operator(): invalid argument size " ); };
336 
337  m_argument[0] = argument;
338  return ( *this )( m_argument );
339  }
340  // ========================================================================
341 
342  // ========================================================================
344  // ========================================================================
345  Genfun::Derivative NumericalDefiniteIntegral::partial( unsigned int idx ) const {
346  if ( idx >= m_DIM ) { Exception( "::partial(i): invalid variable index " ); };
347  //
348  const AbsFunction& aux = NumericalDerivative( *this, idx );
349  return Genfun::FunctionNoop( &aux );
350  }
351  // ========================================================================
352 
353  // ========================================================================
355  // ========================================================================
356  double NumericalDefiniteIntegral::operator()( const Argument& argument ) const {
357  // reset the result and the error
358  m_result = -std::numeric_limits<double>::infinity();
359  m_error = std::numeric_limits<double>::infinity();
360 
361  // check the argument
362  if ( argument.dimension() != m_DIM ) { Exception( "operator(): invalid argument size " ); };
363 
364  // copy the argument
365 
366  for ( size_t i = 0; i < m_DIM; ++i ) {
367  m_argument[i] = argument[i];
368  const size_t j = i < m_index ? i : i + 1;
369  m_argF[j] = argument[i];
370  };
371 
372  // create the helper object
373  GSL_Helper helper( *m_function, m_argF, m_index );
374 
375  // use GSL to evaluate the numerical derivative
376  gsl_function F;
377  F.function = &GSL_Adaptor;
378  F.params = &helper;
379  _Function F1;
380  F1.fn = &F;
381 
382  if ( GaudiMath::Integration::Infinite == category() ) { return QAGI( &F1 ); } // RETURN
383 
384  if ( m_a == m_b ) {
385  m_result = 0;
386  m_error = 0; // EXACT
387  return m_result; // RETURN
388  }
389 
390  switch ( category() ) {
391  case GaudiMath::Integration::Singular:
392  return QAGP( &F1 );
393  case GaudiMath::Integration::Finite:
394  switch ( type() ) {
395  case GaudiMath::Integration::NonAdaptive:
396  return QNG( &F1 );
397  case GaudiMath::Integration::Adaptive:
398  return QAG( &F1 );
399  case GaudiMath::Integration::AdaptiveSingular:
400  return QAGS( &F1 );
401  default:
402  Exception( "::operator(): invalid type " );
403  }
404  break;
405  default:
406  Exception( "::operator(): invalid category " );
407  }
408 
409  return 0;
410  }
411  // ========================================================================
412 
413  // ========================================================================
415  // ========================================================================
416  NumericalDefiniteIntegral::_Workspace* NumericalDefiniteIntegral::allocate() const {
417  if ( !m_ws ) {
418  m_ws.reset( new _Workspace() );
419  m_ws->ws = gsl_integration_workspace_alloc( size() );
420  if ( !m_ws->ws ) { Exception( "allocate()::invalid workspace" ); };
421  }
422  return m_ws.get();
423  }
424  // ========================================================================
425 
426  // ========================================================================
427  // adaptive integration on infinite interval
428  // ========================================================================
429  double NumericalDefiniteIntegral::QAGI( _Function* F ) const {
430  // check the argument
431  if ( !F ) { Exception( "::QAGI: invalid function" ); }
432 
433  // allocate workspace
434  allocate();
435 
436  int ierror = 0;
437 
438  if ( m_a == -std::numeric_limits<double>::infinity() && m_b == std::numeric_limits<double>::infinity() ) {
439  ierror = gsl_integration_qagi( F->fn, m_epsabs, m_epsrel, size(), ws()->ws, &m_result, &m_error );
440  } else if ( m_a == -std::numeric_limits<double>::infinity() ) {
441  ierror = gsl_integration_qagil( F->fn, m_b, m_epsabs, m_epsrel, size(), ws()->ws, &m_result, &m_error );
442  } else if ( m_b == std::numeric_limits<double>::infinity() ) {
443  ierror = gsl_integration_qagiu( F->fn, m_a, m_epsabs, m_epsrel, size(), ws()->ws, &m_result, &m_error );
444  } else {
445  Exception( "::QAGI: invalid mode" );
446  };
447 
448  if ( ierror ) { gsl_error( "NumericalDefiniteIntegral::QAGI", __FILE__, __LINE__, ierror ); }
449 
450  return m_result;
451  }
452  // ========================================================================
453 
454  // ========================================================================
455  // adaptive integration with known singular points
456  // ========================================================================
457  double NumericalDefiniteIntegral::QAGP( _Function* F ) const {
458  if ( !F ) { Exception( "QAGP::invalid function" ); }
459 
460  // no known singular points ?
461  if ( points().empty() ) { return QAGS( F ); }
462 
463  // use GSL
464  int ierror = gsl_integration_qagp( F->fn, const_cast<double*>( points().data() ), points().size(), m_epsabs,
465  m_epsrel, size(), ws()->ws, &m_result, &m_error );
466 
467  if ( ierror ) { gsl_error( "NumericalDefiniteIntegral::QAGI ", __FILE__, __LINE__, ierror ); }
468 
469  // the sign
470  if ( m_b < m_a ) m_result = -m_result;
471 
472  return m_result;
473  }
474  // ========================================================================
475 
476  // ========================================================================
477  // non-adaptive integration
478  // ========================================================================
479  double NumericalDefiniteIntegral::QNG( _Function* F ) const {
480  if ( !F ) { Exception( "QNG::invalid function" ); }
481 
482  // integration limits
483  const auto lim = std::minmax( m_a, m_b );
484 
485  size_t neval = 0;
486  int ierror = gsl_integration_qng( F->fn, lim.first, lim.second, m_epsabs, m_epsrel, &m_result, &m_error, &neval );
487 
488  if ( ierror ) { gsl_error( "NumericalIndefiniteIntegral::QNG ", __FILE__, __LINE__, ierror ); }
489  // sign
490  if ( m_b < m_a ) m_result = -m_result;
491 
492  return m_result;
493  }
494  // ========================================================================
495 
496  // ========================================================================
497  // simple adaptive integration
498  // ========================================================================
499  double NumericalDefiniteIntegral::QAG( _Function* F ) const {
500  if ( !F ) { Exception( "QAG::invalid function" ); }
501 
502  // allocate workspace
503  allocate();
504 
505  // integration limits
506  const auto lim = std::minmax( m_a, m_b );
507 
508  int ierror = gsl_integration_qag( F->fn, lim.first, lim.second, m_epsabs, m_epsrel, size(), (int)rule(), ws()->ws,
509  &m_result, &m_error );
510 
511  if ( ierror ) { gsl_error( "NumericalDefiniteIntegral::QAG ", __FILE__, __LINE__, ierror ); }
512  // sign
513  if ( m_b < m_a ) m_result = -m_result;
514 
515  return m_result;
516  }
517  // ========================================================================
518 
519  // ========================================================================
520  // adaptive integration with singularities
521  // ========================================================================
522  double NumericalDefiniteIntegral::QAGS( _Function* F ) const {
523  if ( !F ) { Exception( "QAG::invalid function" ); }
524 
525  if ( m_a == m_b ) {
526  m_result = 0;
527  m_error = 0; // EXACT !
528  return m_result;
529  }
530 
531  // allocate workspace
532  allocate();
533 
534  // integration limits
535  const auto lim = std::minmax( m_a, m_b );
536 
537  int ierror = gsl_integration_qags( F->fn, lim.first, lim.second, m_epsabs, m_epsrel, size(), ws()->ws, &m_result,
538  &m_error );
539 
540  if ( ierror ) { gsl_error( "NumericalDefiniteIntegral::QAGS ", __FILE__, __LINE__, ierror ); }
541 
542  // sign
543  if ( m_b < m_a ) m_result = -m_result;
544 
545  return m_result;
546  }
547  // ========================================================================
548 
549  } // end of namespace GaudiMathImplementation
550 } // end of namespace Genfun
551 
552 // ============================================================================
553 // The END
554 // ============================================================================
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::category
GaudiMath::Integration::Category category() const
integration category
Definition: NumericalDefiniteIntegral.h:311
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::Exception
StatusCode Exception(const std::string &message, const StatusCode &sc=StatusCode::FAILURE) const
Definition: NumericalDefiniteIntegral.cpp:320
GaudiException
Define general base for Gaudi exception.
Definition: GaudiException.h:21
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::_Workspace
Definition: NumericalDefiniteIntegral.h:80
std::minmax
T minmax(T...args)
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::operator()
double operator()(double argument) const override
Function value.
Definition: NumericalDefiniteIntegral.cpp:329
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::partial
Genfun::Derivative partial(unsigned int index) const override
Derivatives.
Definition: NumericalDefiniteIntegral.cpp:345
std::vector::end
T end(T...args)
GaudiMath::Derivative
Genfun::GaudiMathImplementation::NumericalDerivative Derivative
Definition: GaudiMath.h:26
std::remove_if
T remove_if(T...args)
std::unique
T unique(T...args)
std::string
STL class.
NumericalDerivative.h
std::vector::push_back
T push_back(T...args)
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::dimensionality
unsigned int dimensionality() const override
dimensionality of the problem
Definition: NumericalDefiniteIntegral.h:274
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::_Function::fn
gsl_function * fn
Definition: NumericalDefiniteIntegral.cpp:44
NumericalDefiniteIntegral.h
StatusCode
This class is used for returning status codes from appropriate routines.
Definition: StatusCode.h:50
GaudiMath::Integration::Type
Type
type of integration (for finite limits)
Definition: Integration.h:22
GaudiException.h
std::vector::erase
T erase(T...args)
Genfun::GaudiMathImplementation::GSL_Adaptor
double GSL_Adaptor(double x, void *params)
Definition: Helpers.cpp:25
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::_Function
Definition: NumericalDefiniteIntegral.cpp:43
std::numeric_limits::infinity
T infinity(T...args)
Genfun::GaudiMathImplementation::NumericalDerivative
Numerical derivative (using GSL adaptive numerical differentiation)
Definition: NumericalDerivative.h:42
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::allocate
_Workspace * allocate() const
allocate the integration workspace
Definition: NumericalDefiniteIntegral.cpp:416
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::type
GaudiMath::Integration::Type type() const
integration type
Definition: NumericalDefiniteIntegral.h:309
std::vector::begin
T begin(T...args)
Helpers.h
GaudiMath
GaudiMath.h GaudiMath/GaudiMath.h.
Definition: Adapters.h:13
GaudiMath::Integration::KronrodRule
KronrodRule
integration rule
Definition: Integration.h:26
std::sort
T sort(T...args)
Genfun
CLHEP.
Definition: IEqSolver.h:13
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::function
const AbsFunction & function() const
accessor to the function itself
Definition: NumericalDefiniteIntegral.h:289
Genfun::GaudiMathImplementation::GSL_Helper
the simple structure to be used for adaption interface Genfun::AbsFunction to gsl_function structure ...
Definition: Helpers.h:23
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::rule
GaudiMath::Integration::KronrodRule rule() const
integration rule
Definition: NumericalDefiniteIntegral.h:313
Integration
collection of common types for classes NumericalIndefiniteIntegral and NumericalDefiniteIntegral ...
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral
This class allows the numerical evaluation of the following functions:
Definition: NumericalDefiniteIntegral.h:78