The Gaudi Framework  v30r0 (c919700c)
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 {
42  namespace GaudiMathImplementation
43  {
44 
45  struct NumericalDefiniteIntegral::_Function {
46  gsl_function* fn;
47  };
48 
49  // ========================================================================
50 
102  NumericalDefiniteIntegral::NumericalDefiniteIntegral( const AbsFunction& function, const size_t index,
103  const double a, const double b,
104  const GaudiMath::Integration::Type type,
105  const GaudiMath::Integration::KronrodRule rule,
106  const double epsabs, const double epsrel, const size_t size )
107  : m_function( function.clone() )
108  , m_index( index )
109  , m_a( a )
110  , m_b( b )
111  , m_type( type )
112  , m_category( GaudiMath::Integration::Finite )
113  , m_rule( rule )
114  , m_epsabs( epsabs )
115  , m_epsrel( epsrel )
116  , m_size( size )
117  {
118  if ( GaudiMath::Integration::Fixed == m_rule ) {
119  m_rule = GaudiMath::Integration::Default;
120  }
121  if ( function.dimensionality() < 2 ) {
122  Exception( "::constructor: invalid dimensionality " );
123  }
124  if ( m_index >= function.dimensionality() ) {
125  Exception( "::constructor: invalid variable index" );
126  }
127 
128  m_DIM = function.dimensionality() - 1;
129  m_argument = Argument( m_DIM );
130  m_argF = Argument( m_DIM + 1 );
131  }
132 
142  NumericalDefiniteIntegral::NumericalDefiniteIntegral( const AbsFunction& function, const size_t index,
143  const double a, const double b,
144  const NumericalDefiniteIntegral::Points& points,
145  const double epsabs, const double epsrel, const size_t size )
146  : m_function( function.clone() )
147  , m_index( index )
148  , m_a( a )
149  , m_b( b )
150  , m_type( GaudiMath::Integration::Other )
151  , m_category( GaudiMath::Integration::Singular )
152  , m_rule( GaudiMath::Integration::Fixed )
153  , m_points( points )
154  , m_epsabs( epsabs )
155  , m_epsrel( epsrel )
156  , m_size( size )
157  {
158  if ( function.dimensionality() < 2 ) {
159  Exception( "::constructor: invalid dimensionality " );
160  }
161  if ( m_index >= function.dimensionality() ) {
162  Exception( "::constructor: invalid variable index" );
163  }
164 
165  m_DIM = function.dimensionality() - 1;
166  m_argument = Argument( m_DIM );
167  m_argF = Argument( m_DIM + 1 );
168 
169  const auto lim = std::minmax( a, b );
170  m_points.push_back( lim.first );
171  m_points.push_back( lim.second );
172  auto end = std::remove_if( m_points.begin(), m_points.end(),
173  [&]( double x ) { return x < lim.first || lim.second < x; } );
174  std::sort( m_points.begin(), end );
175  m_points.erase( std::unique( m_points.begin(), end ), m_points.end() );
176  }
177 
203  NumericalDefiniteIntegral::NumericalDefiniteIntegral( const AbsFunction& function, const size_t index,
204  const double a, const GaudiMath::Integration::Inf /* b */,
205  const double epsabs, const double epsrel, const size_t size )
206  : m_function( function.clone() )
207  , m_DIM( 0 )
208  , m_index( index )
209  , m_a( a )
210  , m_type( GaudiMath::Integration::Other )
211  , m_category( GaudiMath::Integration::Infinite )
212  , m_rule( GaudiMath::Integration::Fixed )
213  , m_epsabs( epsabs )
214  , m_epsrel( epsrel )
215  , m_size( size )
216  {
217  if ( function.dimensionality() < 2 ) {
218  Exception( "::constructor: invalid dimensionality " );
219  }
220  if ( m_index >= function.dimensionality() ) {
221  Exception( "::constructor: invalid variable index" );
222  }
223 
224  m_DIM = function.dimensionality() - 1;
225  m_argument = Argument( m_DIM );
226  m_argF = Argument( m_DIM + 1 );
227  }
228 
254  NumericalDefiniteIntegral::NumericalDefiniteIntegral( const AbsFunction& function, const size_t index,
255  const GaudiMath::Integration::Inf /* a */, const double b,
256  const double epsabs, const double epsrel, const size_t size )
257  : m_function( function.clone() )
258  , m_DIM( 0 )
259  , m_index( index )
260  , m_b( b )
261  , m_type( GaudiMath::Integration::Other )
262  , m_category( GaudiMath::Integration::Infinite )
263  , m_rule( GaudiMath::Integration::Fixed )
264  , m_epsabs( epsabs )
265  , m_epsrel( epsrel )
266  //
267  , m_size( size )
268  {
269  if ( function.dimensionality() < 2 ) {
270  Exception( "::constructor: invalid dimensionality " );
271  }
272  if ( m_index >= function.dimensionality() ) {
273  Exception( "::constructor: invalid variable index" );
274  }
275 
276  m_DIM = function.dimensionality() - 1;
277  m_argument = Argument( m_DIM );
278  m_argF = Argument( m_DIM + 1 );
279  }
280 
304  NumericalDefiniteIntegral::NumericalDefiniteIntegral( const AbsFunction& function, const size_t index,
305  const float epsabs, const float epsrel, const size_t size )
306  : m_function( function.clone() )
307  , m_DIM( 0 )
308  , m_index( index )
309  , m_type( GaudiMath::Integration::Other )
310  , m_category( GaudiMath::Integration::Infinite )
311  , m_rule( GaudiMath::Integration::Fixed )
312  , m_epsabs( epsabs )
313  , m_epsrel( epsrel )
314  , m_size( size )
315  {
316  if ( function.dimensionality() < 2 ) {
317  Exception( "::constructor: invalid dimensionality " );
318  }
319  if ( m_index >= function.dimensionality() ) {
320  Exception( "::constructor: invalid variable index" );
321  }
322 
323  m_DIM = function.dimensionality() - 1;
324  m_argument = Argument( m_DIM );
325  m_argF = Argument( m_DIM + 1 );
326  }
327 
329  NumericalDefiniteIntegral::NumericalDefiniteIntegral( const NumericalDefiniteIntegral& right )
330  : AbsFunction()
331  , m_function( right.m_function->clone() )
332  , m_DIM( right.m_DIM )
333  , m_index( right.m_index )
334  , m_a( right.m_a )
335  , m_b( right.m_b )
336  , m_type( right.m_type )
337  , m_category( right.m_category )
338  , m_rule( right.m_rule )
339  , m_points( right.m_points )
340  , m_epsabs( right.m_epsabs )
341  , m_epsrel( right.m_epsrel )
342  , m_size( right.m_size )
343  , m_argument( right.m_argument )
344  , m_argF( right.m_argF )
345  {
346  }
347 
348  // ========================================================================
349  // throw the exception
350  // ========================================================================
351  StatusCode NumericalDefiniteIntegral::Exception( const std::string& message, const StatusCode& sc ) const
352  {
353  throw GaudiException( "NumericalDefiniteIntegral::" + message, "*GaudiMath*", sc );
354  return sc;
355  }
356  // ========================================================================
357 
358  // ========================================================================
360  // ========================================================================
361  double NumericalDefiniteIntegral::operator()( const double argument ) const
362  {
363  // reset the result and the error
364  m_result = -std::numeric_limits<double>::infinity();
365  m_error = std::numeric_limits<double>::infinity();
366 
367  // check the argument
368  if ( 1 != m_DIM ) {
369  Exception( "operator(): invalid argument size " );
370  };
371 
372  m_argument[0] = argument;
373  return ( *this )( m_argument );
374  }
375  // ========================================================================
376 
377  // ========================================================================
379  // ========================================================================
380  Genfun::Derivative NumericalDefiniteIntegral::partial( unsigned int idx ) const
381  {
382  if ( idx >= m_DIM ) {
383  Exception( "::partial(i): invalid variable index " );
384  };
385  //
386  const AbsFunction& aux = NumericalDerivative( *this, idx );
387  return Genfun::FunctionNoop( &aux );
388  }
389  // ========================================================================
390 
391  // ========================================================================
393  // ========================================================================
394  double NumericalDefiniteIntegral::operator()( const Argument& argument ) const
395  {
396  // reset the result and the error
397  m_result = -std::numeric_limits<double>::infinity();
398  m_error = std::numeric_limits<double>::infinity();
399 
400  // check the argument
401  if ( argument.dimension() != m_DIM ) {
402  Exception( "operator(): invalid argument size " );
403  };
404 
405  // copy the argument
406 
407  for ( size_t i = 0; i < m_DIM; ++i ) {
408  m_argument[i] = argument[i];
409  const size_t j = i < m_index ? i : i + 1;
410  m_argF[j] = argument[i];
411  };
412 
413  // create the helper object
414  GSL_Helper helper( *m_function, m_argF, m_index );
415 
416  // use GSL to evaluate the numerical derivative
417  gsl_function F;
418  F.function = &GSL_Adaptor;
419  F.params = &helper;
420  _Function F1;
421  F1.fn = &F;
422 
423  if ( GaudiMath::Integration::Infinite == category() ) {
424  return QAGI( &F1 );
425  } // RETURN
426 
427  if ( m_a == m_b ) {
428  m_result = 0;
429  m_error = 0; // EXACT
430  return m_result; // RETURN
431  }
432 
433  switch ( category() ) {
434  case GaudiMath::Integration::Singular:
435  return QAGP( &F1 );
436  case GaudiMath::Integration::Finite:
437  switch ( type() ) {
438  case GaudiMath::Integration::NonAdaptive:
439  return QNG( &F1 );
440  case GaudiMath::Integration::Adaptive:
441  return QAG( &F1 );
442  case GaudiMath::Integration::AdaptiveSingular:
443  return QAGS( &F1 );
444  default:
445  Exception( "::operator(): invalid type " );
446  }
447  break;
448  default:
449  Exception( "::operator(): invalid category " );
450  }
451 
452  return 0;
453  }
454  // ========================================================================
455 
456  // ========================================================================
458  // ========================================================================
459  NumericalDefiniteIntegral::_Workspace* NumericalDefiniteIntegral::allocate() const
460  {
461  if ( !m_ws ) {
462  m_ws.reset( new _Workspace() );
463  m_ws->ws = gsl_integration_workspace_alloc( size() );
464  if ( !m_ws->ws ) {
465  Exception( "allocate()::invalid workspace" );
466  };
467  }
468  return m_ws.get();
469  }
470  // ========================================================================
471 
472  // ========================================================================
473  // adaptive integration on infinite interval
474  // ========================================================================
475  double NumericalDefiniteIntegral::QAGI( _Function* F ) const
476  {
477  // check the argument
478  if ( !F ) {
479  Exception( "::QAGI: invalid function" );
480  }
481 
482  // allocate workspace
483  allocate();
484 
485  int ierror = 0;
486 
487  if ( m_a == -std::numeric_limits<double>::infinity() && m_b == std::numeric_limits<double>::infinity() ) {
488  ierror = gsl_integration_qagi( F->fn, m_epsabs, m_epsrel, size(), ws()->ws, &m_result, &m_error );
489  } else if ( m_a == -std::numeric_limits<double>::infinity() ) {
490  ierror = gsl_integration_qagil( F->fn, m_b, m_epsabs, m_epsrel, size(), ws()->ws, &m_result, &m_error );
491  } else if ( m_b == std::numeric_limits<double>::infinity() ) {
492  ierror = gsl_integration_qagiu( F->fn, m_a, m_epsabs, m_epsrel, size(), ws()->ws, &m_result, &m_error );
493  } else {
494  Exception( "::QAGI: invalid mode" );
495  };
496 
497  if ( ierror ) {
498  gsl_error( "NumericalDefiniteIntegral::QAGI", __FILE__, __LINE__, ierror );
499  }
500 
501  return m_result;
502  }
503  // ========================================================================
504 
505  // ========================================================================
506  // adaptive integration with known singular points
507  // ========================================================================
508  double NumericalDefiniteIntegral::QAGP( _Function* F ) const
509  {
510  if ( !F ) {
511  Exception( "QAGP::invalid function" );
512  }
513 
514  // no known singular points ?
515  if ( points().empty() ) {
516  return QAGS( F );
517  }
518 
519  // use GSL
520  int ierror = gsl_integration_qagp( F->fn, const_cast<double*>( points().data() ), points().size(), m_epsabs,
521  m_epsrel, size(), ws()->ws, &m_result, &m_error );
522 
523  if ( ierror ) {
524  gsl_error( "NumericalDefiniteIntegral::QAGI ", __FILE__, __LINE__, ierror );
525  }
526 
527  // the sign
528  if ( m_b < m_a ) m_result = -m_result;
529 
530  return m_result;
531  }
532  // ========================================================================
533 
534  // ========================================================================
535  // non-adaptive integration
536  // ========================================================================
537  double NumericalDefiniteIntegral::QNG( _Function* F ) const
538  {
539  if ( !F ) {
540  Exception( "QNG::invalid function" );
541  }
542 
543  // integration limits
544  const auto lim = std::minmax( m_a, m_b );
545 
546  size_t neval = 0;
547  int ierror = gsl_integration_qng( F->fn, lim.first, lim.second, m_epsabs, m_epsrel, &m_result, &m_error, &neval );
548 
549  if ( ierror ) {
550  gsl_error( "NumericalIndefiniteIntegral::QNG ", __FILE__, __LINE__, ierror );
551  }
552  // sign
553  if ( m_b < m_a ) m_result = -m_result;
554 
555  return m_result;
556  }
557  // ========================================================================
558 
559  // ========================================================================
560  // simple adaptive integration
561  // ========================================================================
562  double NumericalDefiniteIntegral::QAG( _Function* F ) const
563  {
564  if ( !F ) {
565  Exception( "QAG::invalid function" );
566  }
567 
568  // allocate workspace
569  allocate();
570 
571  // integration limits
572  const auto lim = std::minmax( m_a, m_b );
573 
574  int ierror = gsl_integration_qag( F->fn, lim.first, lim.second, m_epsabs, m_epsrel, size(), (int)rule(), ws()->ws,
575  &m_result, &m_error );
576 
577  if ( ierror ) {
578  gsl_error( "NumericalDefiniteIntegral::QAG ", __FILE__, __LINE__, ierror );
579  }
580  // sign
581  if ( m_b < m_a ) m_result = -m_result;
582 
583  return m_result;
584  }
585  // ========================================================================
586 
587  // ========================================================================
588  // adaptive integration with singularities
589  // ========================================================================
590  double NumericalDefiniteIntegral::QAGS( _Function* F ) const
591  {
592  if ( !F ) {
593  Exception( "QAG::invalid function" );
594  }
595 
596  if ( m_a == m_b ) {
597  m_result = 0;
598  m_error = 0; // EXACT !
599  return m_result;
600  }
601 
602  // allocate workspace
603  allocate();
604 
605  // integration limits
606  const auto lim = std::minmax( m_a, m_b );
607 
608  int ierror = gsl_integration_qags( F->fn, lim.first, lim.second, m_epsabs, m_epsrel, size(), ws()->ws, &m_result,
609  &m_error );
610 
611  if ( ierror ) {
612  gsl_error( "NumericalDefiniteIntegral::QAGS ", __FILE__, __LINE__, ierror );
613  }
614 
615  // sign
616  if ( m_b < m_a ) m_result = -m_result;
617 
618  return m_result;
619  }
620  // ========================================================================
621 
622  } // end of namespace GaudiMathImplementation
623 } // end of namespace Genfun
624 
625 // ============================================================================
626 // The END
627 // ============================================================================
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::category
GaudiMath::Integration::Category category() const
integration category
Definition: NumericalDefiniteIntegral.h:316
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::Exception
StatusCode Exception(const std::string &message, const StatusCode &sc=StatusCode::FAILURE) const
Definition: NumericalDefiniteIntegral.cpp:351
GaudiException
Define general base for Gaudi exception.
Definition: GaudiException.h:22
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::_Workspace
Definition: NumericalDefiniteIntegral.h:82
std::minmax
T minmax(T...args)
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::operator()
double operator()(double argument) const override
Function value.
Definition: NumericalDefiniteIntegral.cpp:361
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::partial
Genfun::Derivative partial(unsigned int index) const override
Derivatives.
Definition: NumericalDefiniteIntegral.cpp:380
std::vector::end
T end(T...args)
GaudiMath::Derivative
Genfun::GaudiMathImplementation::NumericalDerivative Derivative
Definition: GaudiMath.h:27
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:279
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::_Function::fn
gsl_function * fn
Definition: NumericalDefiniteIntegral.cpp:46
NumericalDefiniteIntegral.h
StatusCode
This class is used for returning status codes from appropriate routines.
Definition: StatusCode.h:26
GaudiMath::Integration::Type
Type
type of integration (for finite limits)
Definition: Integration.h:24
GaudiException.h
std::vector::erase
T erase(T...args)
end
auto end(reverse_wrapper< T > &w)
Definition: reverse.h:64
Genfun::GaudiMathImplementation::GSL_Adaptor
double GSL_Adaptor(double x, void *params)
Definition: Helpers.cpp:31
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::_Function
Definition: NumericalDefiniteIntegral.cpp:45
std::numeric_limits::infinity
T infinity(T...args)
Genfun::GaudiMathImplementation::NumericalDerivative
Numerical derivative (using GSL adaptive numerical differentiation)
Definition: NumericalDerivative.h:44
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::allocate
_Workspace * allocate() const
allocate the integration workspace
Definition: NumericalDefiniteIntegral.cpp:459
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::type
GaudiMath::Integration::Type type() const
integration type
Definition: NumericalDefiniteIntegral.h:314
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:28
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:294
Genfun::GaudiMathImplementation::GSL_Helper
the simple structure to be used for adaption interface Genfun::AbsFunction to gsl_function structure ...
Definition: Helpers.h:25
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::rule
GaudiMath::Integration::KronrodRule rule() const
integration rule
Definition: NumericalDefiniteIntegral.h:318
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:79