The Gaudi Framework  v28r3 (cc1cf868)
NumericalDefiniteIntegral.cpp
Go to the documentation of this file.
1 // ============================================================================
2 // Include
3 // ============================================================================
4 // STD & STL
5 // ============================================================================
6 #include <vector>
7 #include <algorithm>
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 
29 #ifdef __ICC
30 // disable icc remark #1572: floating-point equality and inequality comparisons are unreliable
31 // The comparison are meant
32 #pragma warning(disable:1572)
33 #endif
34 #ifdef WIN32
35 // Disable the warning
36 // C4996: 'std::copy': Function call with parameters that may be unsafe
37 // The parameters are checked
38 #pragma warning(disable:4996)
39 #endif
40 
41 namespace Genfun
42 {
43  namespace GaudiMathImplementation
44  {
45 
46  struct NumericalDefiniteIntegral::_Function
47  { gsl_function* fn ; };
48 
49  // ========================================================================
50 
102  NumericalDefiniteIntegral::NumericalDefiniteIntegral
103  ( const AbsFunction& function ,
104  const size_t index ,
105  const double a ,
106  const double b ,
107  const GaudiMath::Integration::Type type ,
108  const GaudiMath::Integration::KronrodRule rule ,
109  const double epsabs ,
110  const double epsrel ,
111  const size_t size )
112  : m_function ( function.clone () )
113  , m_index ( index )
114  , m_a ( a )
115  , m_b ( b )
116  , m_type ( type )
117  , m_category ( GaudiMath::Integration::Finite )
118  , m_rule ( rule )
119  , m_epsabs ( epsabs )
120  , m_epsrel ( epsrel )
121  , m_size ( size )
122  {
123  if ( GaudiMath::Integration::Fixed == m_rule )
124  { m_rule = GaudiMath::Integration::Default ; }
125  if ( function.dimensionality() < 2 )
126  { Exception("::constructor: invalid dimensionality ") ; }
127  if ( m_index >= function.dimensionality() )
128  { Exception("::constructor: invalid variable index") ; }
129 
130  m_DIM = function.dimensionality() - 1 ;
131  m_argument = Argument( m_DIM ) ;
132  m_argF = Argument( m_DIM + 1 ) ;
133 
134  }
135 
136 
146  NumericalDefiniteIntegral::NumericalDefiniteIntegral
147  ( const AbsFunction& function ,
148  const size_t index ,
149  const double a ,
150  const double b ,
151  const NumericalDefiniteIntegral::Points& points ,
152  const double epsabs ,
153  const double epsrel ,
154  const size_t size )
155  : m_function ( function.clone() )
156  , m_index ( index )
157  , m_a ( a )
158  , m_b ( b )
159  , m_type ( GaudiMath::Integration:: Other )
160  , m_category ( GaudiMath::Integration:: Singular )
161  , m_rule ( GaudiMath::Integration:: Fixed )
162  , m_points ( points )
163  , m_epsabs ( epsabs )
164  , m_epsrel ( epsrel )
165  , m_size ( size )
166  {
167  if ( function.dimensionality() < 2 )
168  { Exception("::constructor: invalid dimensionality ") ; }
169  if ( m_index >= function.dimensionality() )
170  { Exception("::constructor: invalid variable index") ; }
171 
172  m_DIM = function.dimensionality() - 1 ;
173  m_argument = Argument( m_DIM ) ;
174  m_argF = Argument( m_DIM + 1 ) ;
175 
176  const auto lim = std::minmax( a, b );
177  m_points.push_back ( lim.first ) ;
178  m_points.push_back ( lim.second ) ;
179  auto end = std::remove_if( m_points.begin(), m_points.end(),
180  [&](double x) { return x<lim.first || lim.second<x; } );
181  std::sort( m_points.begin() , end );
182  m_points.erase( std::unique( m_points.begin () , end ) ,
183  m_points.end() );
184 
185  }
186 
212  NumericalDefiniteIntegral::NumericalDefiniteIntegral
213  ( const AbsFunction& function ,
214  const size_t index ,
215  const double a ,
216  const GaudiMath::Integration::Inf /* b */ ,
217  const double epsabs ,
218  const double epsrel ,
219  const size_t size )
220  : m_function ( function.clone() )
221  , m_DIM ( 0 )
222  , m_index ( index )
223  , m_a ( a )
224  , m_type ( GaudiMath::Integration:: Other )
225  , m_category ( GaudiMath::Integration:: Infinite )
226  , m_rule ( GaudiMath::Integration:: Fixed )
227  , m_epsabs ( epsabs )
228  , m_epsrel ( epsrel )
229  , m_size ( size )
230  {
231  if ( function.dimensionality() < 2 )
232  { Exception("::constructor: invalid dimensionality ") ; }
233  if ( m_index >= function.dimensionality() )
234  { Exception("::constructor: invalid variable index") ; }
235 
236  m_DIM = function.dimensionality() - 1 ;
237  m_argument = Argument( m_DIM ) ;
238  m_argF = Argument( m_DIM + 1 ) ;
239 
240  }
241 
242 
268  NumericalDefiniteIntegral::NumericalDefiniteIntegral
269  ( const AbsFunction& function ,
270  const size_t index ,
271  const GaudiMath::Integration::Inf /* a */ ,
272  const double b ,
273  const double epsabs ,
274  const double epsrel ,
275  const size_t size )
276  : m_function ( function.clone() )
277  , m_DIM ( 0 )
278  , m_index ( index )
279  , m_b ( b )
280  , m_type ( GaudiMath::Integration:: Other )
281  , m_category ( GaudiMath::Integration:: Infinite )
282  , m_rule ( GaudiMath::Integration:: Fixed )
283  , m_epsabs ( epsabs )
284  , m_epsrel ( epsrel )
285  //
286  , m_size ( size )
287  {
288  if ( function.dimensionality() < 2 )
289  { Exception("::constructor: invalid dimensionality ") ; }
290  if ( m_index >= function.dimensionality() )
291  { Exception("::constructor: invalid variable index") ; }
292 
293  m_DIM = function.dimensionality() - 1 ;
294  m_argument = Argument( m_DIM ) ;
295  m_argF = Argument( m_DIM + 1 ) ;
296  }
297 
321  NumericalDefiniteIntegral::NumericalDefiniteIntegral
322  ( const AbsFunction& function ,
323  const size_t index ,
324  const float epsabs ,
325  const float epsrel ,
326  const size_t size )
327  : m_function ( function.clone() )
328  , m_DIM ( 0 )
329  , m_index ( index )
330  , m_type ( GaudiMath::Integration:: Other )
331  , m_category ( GaudiMath::Integration:: Infinite )
332  , m_rule ( GaudiMath::Integration:: Fixed )
333  , m_epsabs ( epsabs )
334  , m_epsrel ( epsrel )
335  , m_size ( size )
336  {
337  if ( function.dimensionality() < 2 )
338  { Exception("::constructor: invalid dimensionality ") ; }
339  if ( m_index >= function.dimensionality() )
340  { Exception("::constructor: invalid variable index") ; }
341 
342  m_DIM = function.dimensionality() - 1 ;
343  m_argument = Argument( m_DIM ) ;
344  m_argF = Argument( m_DIM + 1 ) ;
345 
346  }
347 
349  NumericalDefiniteIntegral::NumericalDefiniteIntegral
350  ( const NumericalDefiniteIntegral& right )
351  : AbsFunction()
352  , m_function ( right.m_function->clone() )
353  , m_DIM ( right.m_DIM )
354  , m_index ( right.m_index )
355  , m_a ( right.m_a )
356  , m_b ( right.m_b )
357  , m_type ( right.m_type )
358  , m_category ( right.m_category )
359  , m_rule ( right.m_rule )
360  , m_points ( right.m_points )
361  , m_epsabs ( right.m_epsabs )
362  , m_epsrel ( right.m_epsrel )
363  , m_size ( right.m_size )
364  , m_argument ( right.m_argument )
365  , m_argF ( right.m_argF )
366  {
367  }
368 
369  // ========================================================================
370  // throw the exception
371  // ========================================================================
372  StatusCode NumericalDefiniteIntegral::Exception
373  ( const std::string& message ,
374  const StatusCode& sc ) const
375  {
376  throw GaudiException( "NumericalDefiniteIntegral::" + message ,
377  "*GaudiMath*" , sc ) ;
378  return sc ;
379  }
380  // ========================================================================
381 
382  // ========================================================================
384  // ========================================================================
385  double NumericalDefiniteIntegral::operator()( const double argument ) const
386  {
387  // reset the result and the error
388  m_result = -std::numeric_limits<double>::infinity();
389  m_error = std::numeric_limits<double>::infinity();
390 
391  // check the argument
392  if( 1 != m_DIM ) { Exception ( "operator(): invalid argument size " ) ; };
393 
394  m_argument[0] = argument ;
395  return (*this) ( m_argument );
396  }
397  // ========================================================================
398 
399  // ========================================================================
401  // ========================================================================
402  Genfun::Derivative
403  NumericalDefiniteIntegral::partial ( unsigned int idx ) const
404  {
405  if ( idx >= m_DIM )
406  { Exception ( "::partial(i): invalid variable index " ) ; };
407  //
408  const AbsFunction& aux = NumericalDerivative( *this , idx ) ;
409  return Genfun::FunctionNoop( &aux ) ;
410  }
411  // ========================================================================
412 
413 
414  // ========================================================================
416  // ========================================================================
417  double NumericalDefiniteIntegral::operator()( const Argument& argument ) const
418  {
419  // reset the result and the error
420  m_result = -std::numeric_limits<double>::infinity();
421  m_error = std::numeric_limits<double>::infinity() ;
422 
423  // check the argument
424  if( argument.dimension() != m_DIM )
425  { Exception ( "operator(): invalid argument size " ) ; };
426 
427  // copy the argument
428 
429  for( size_t i = 0 ; i < m_DIM ; ++i ) {
430  m_argument [i] = argument [i] ;
431  const size_t j = i < m_index ? i : i + 1 ;
432  m_argF [j] = argument [i] ;
433  };
434 
435  // create the helper object
436  GSL_Helper helper( *m_function , m_argF , m_index );
437 
438  // use GSL to evaluate the numerical derivative
439  gsl_function F ;
440  F.function = &GSL_Adaptor ;
441  F.params = &helper ;
442  _Function F1 ;
443  F1.fn = &F ;
444 
445  if ( GaudiMath::Integration::Infinite == category () )
446  { return QAGI ( &F1 ) ; } // RETURN
447 
448  if ( m_a == m_b ) {
449  m_result = 0 ; m_error = 0 ; // EXACT
450  return m_result ; // RETURN
451  }
452 
453  switch( category() ) {
454  case GaudiMath::Integration::Singular : return QAGP ( &F1 ) ;
455  case GaudiMath::Integration::Finite :
456  switch (type()) {
457  case GaudiMath::Integration::NonAdaptive :return QNG( &F1 ) ;
458  case GaudiMath::Integration::Adaptive :return QAG( &F1 ) ;
459  case GaudiMath::Integration::AdaptiveSingular: return QAGS ( &F1 ) ;
460  default: Exception ( "::operator(): invalid type " );
461  }
462  break;
463  default: Exception ( "::operator(): invalid category " );
464  }
465 
466  return 0 ;
467  }
468  // ========================================================================
469 
470  // ========================================================================
472  // ========================================================================
473  NumericalDefiniteIntegral::_Workspace*
474  NumericalDefiniteIntegral::allocate() const
475  {
476  if ( !m_ws ) {
477  m_ws.reset( new _Workspace() );
478  m_ws->ws = gsl_integration_workspace_alloc( size () );
479  if ( !m_ws->ws ) { Exception ( "allocate()::invalid workspace" ) ; };
480  }
481  return m_ws.get() ;
482  }
483  // ========================================================================
484 
485  // ========================================================================
486  // adaptive integration on infinite interval
487  // ========================================================================
488  double NumericalDefiniteIntegral::QAGI ( _Function* F ) const
489  {
490  // check the argument
491  if ( !F ) { Exception("::QAGI: invalid function"); }
492 
493  // allocate workspace
494  allocate() ;
495 
496  int ierror = 0 ;
497 
498  if ( m_a == -std::numeric_limits<double>::infinity() && m_b == std::numeric_limits<double>::infinity() ) {
499  ierror = gsl_integration_qagi ( F->fn ,
500  m_epsabs , m_epsrel ,
501  size () , ws()->ws ,
502  &m_result , &m_error ) ;
503  } else if ( m_a == -std::numeric_limits<double>::infinity() ) {
504  ierror = gsl_integration_qagil ( F->fn , m_b ,
505  m_epsabs , m_epsrel ,
506  size () , ws()->ws ,
507  &m_result , &m_error ) ;
508  } else if ( m_b == std::numeric_limits<double>::infinity() ) {
509  ierror = gsl_integration_qagiu ( F->fn , m_a ,
510  m_epsabs , m_epsrel ,
511  size () , ws()->ws ,
512  &m_result , &m_error ) ;
513  } else { Exception ( "::QAGI: invalid mode" ) ; };
514 
515  if( ierror ) { gsl_error( "NumericalDefiniteIntegral::QAGI" ,
516  __FILE__ , __LINE__ , ierror ) ;}
517 
518  return m_result ;
519  }
520  // ========================================================================
521 
522  // ========================================================================
523  // adaptive integration with known singular points
524  // ========================================================================
525  double NumericalDefiniteIntegral::QAGP( _Function* F ) const
526  {
527  if( !F ) { Exception("QAGP::invalid function") ; }
528 
529  // no known singular points ?
530  if( points().empty() ) { return QAGS( F ) ; }
531 
532  // use GSL
533  int ierror =
534  gsl_integration_qagp ( F->fn ,
535  const_cast<double*>(points().data()) , points().size() ,
536  m_epsabs , m_epsrel ,
537  size () , ws()->ws ,
538  &m_result , &m_error ) ;
539 
540  if( ierror ) { gsl_error( "NumericalDefiniteIntegral::QAGI " ,
541  __FILE__ , __LINE__ , ierror ) ; }
542 
543  // the sign
544  if ( m_b < m_a ) m_result = -m_result;
545 
546  return m_result;
547  }
548  // ========================================================================
549 
550  // ========================================================================
551  // non-adaptive integration
552  // ========================================================================
553  double NumericalDefiniteIntegral::QNG ( _Function* F ) const
554  {
555  if( !F ) { Exception("QNG::invalid function") ; }
556 
557  // integration limits
558  const auto lim = std::minmax( m_a, m_b );
559 
560  size_t neval = 0 ;
561  int ierror =
562  gsl_integration_qng ( F->fn ,
563  lim.first , lim.second,
564  m_epsabs , m_epsrel ,
565  &m_result , &m_error , &neval ) ;
566 
567  if( ierror ) { gsl_error( "NumericalIndefiniteIntegral::QNG " ,
568  __FILE__ , __LINE__ , ierror ) ; }
569  // sign
570  if ( m_b < m_a ) m_result = -m_result;
571 
572  return m_result;
573  }
574  // ========================================================================
575 
576 
577  // ========================================================================
578  // simple adaptive integration
579  // ========================================================================
580  double NumericalDefiniteIntegral::QAG ( _Function* F ) const
581  {
582  if( !F ) { Exception("QAG::invalid function") ; }
583 
584  // allocate workspace
585  allocate () ;
586 
587  // integration limits
588  const auto lim = std::minmax( m_a, m_b );
589 
590  int ierror =
591  gsl_integration_qag ( F->fn ,
592  lim.first , lim.second,
593  m_epsabs , m_epsrel ,
594  size () , (int) rule() , ws ()->ws ,
595  &m_result , &m_error );
596 
597  if( ierror ) { gsl_error( "NumericalDefiniteIntegral::QAG " ,
598  __FILE__ , __LINE__ , ierror ) ; }
599  // sign
600  if ( m_b < m_a ) m_result = -m_result;
601 
602  return m_result;
603  }
604  // ========================================================================
605 
606  // ========================================================================
607  // adaptive integration with singularities
608  // ========================================================================
609  double NumericalDefiniteIntegral::QAGS ( _Function* F ) const
610  {
611  if( !F ) { Exception("QAG::invalid function") ; }
612 
613  if ( m_a == m_b ) {
614  m_result = 0 ;
615  m_error = 0 ; // EXACT !
616  return m_result ;
617  }
618 
619  // allocate workspace
620  allocate () ;
621 
622  // integration limits
623  const auto lim = std::minmax( m_a, m_b );
624 
625  int ierror =
626  gsl_integration_qags ( F->fn ,
627  lim.first , lim.second,
628  m_epsabs , m_epsrel ,
629  size () , ws()->ws ,
630  &m_result , &m_error );
631 
632  if( ierror ) { gsl_error( "NumericalDefiniteIntegral::QAGS " ,
633  __FILE__ , __LINE__ , ierror ) ; }
634 
635  // sign
636  if ( m_b < m_a ) m_result = -m_result;
637 
638  return m_result;
639  }
640  // ========================================================================
641 
642  } // end of namespace GaudiMathImplementation
643 } // end of namespace Genfun
644 
645 // ============================================================================
646 // The END
647 // ============================================================================
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::category
GaudiMath::Integration::Category category() const
integration category
Definition: NumericalDefiniteIntegral.h:351
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::Exception
StatusCode Exception(const std::string &message, const StatusCode &sc=StatusCode::FAILURE) const
Definition: NumericalDefiniteIntegral.cpp:373
GaudiException
Define general base for Gaudi exception.
Definition: GaudiException.h:23
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:385
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::partial
Genfun::Derivative partial(unsigned int index) const override
Derivatives.
Definition: NumericalDefiniteIntegral.cpp:403
std::vector::end
T end(T...args)
GaudiMath::Derivative
Genfun::GaudiMathImplementation::NumericalDerivative Derivative
Definition: GaudiMath.h:29
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:311
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::_Function::fn
gsl_function * fn
Definition: NumericalDefiniteIntegral.cpp:47
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:25
GaudiException.h
std::vector::erase
T erase(T...args)
end
auto end(reverse_wrapper< T > &w)
Definition: reverse.h:50
Genfun::GaudiMathImplementation::GSL_Adaptor
double GSL_Adaptor(double x, void *params)
Definition: Helpers.cpp:33
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::_Function
Definition: NumericalDefiniteIntegral.cpp:46
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:474
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::type
GaudiMath::Integration::Type type() const
integration type
Definition: NumericalDefiniteIntegral.h:348
std::vector::begin
T begin(T...args)
Helpers.h
GaudiMath::Integration::KronrodRule
KronrodRule
integration rule
Definition: Integration.h:34
std::sort
T sort(T...args)
Genfun
CLHEP.
Definition: IEqSolver.h:13
Genfun::GaudiMathImplementation::GSL_Helper
the simple structure to be used for adaption interface Genfun::AbsFunction to gsl_function structure ...
Definition: Helpers.h:26
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::rule
GaudiMath::Integration::KronrodRule rule() const
integration rule
Definition: NumericalDefiniteIntegral.h:354
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral
This class allows the numerical evaluation of the following functions:
Definition: NumericalDefiniteIntegral.h:79