•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 <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_DIM ( 0 )
114  , m_index ( index )
115  , m_a ( a )
116  , m_b ( b )
117  , m_ia ( false )
118  , m_ib ( false )
119  , m_type ( type )
120  , m_category ( GaudiMath::Integration::Finite )
121  , m_rule ( rule )
122  , m_epsabs ( epsabs )
123  , m_epsrel ( epsrel )
124  , m_result ( GSL_NEGINF )
125  , m_error ( GSL_POSINF )
126  , m_size ( size )
127  {
128  if ( GaudiMath::Integration::Fixed == m_rule )
129  { m_rule = GaudiMath::Integration::Default ; }
130  if ( function.dimensionality() < 2 )
131  { Exception("::constructor: invalid dimensionality ") ; }
132  if ( m_index >= function.dimensionality() )
133  { Exception("::constructor: invalid variable index") ; }
134 
135  m_DIM = function.dimensionality() - 1 ;
136  m_argument = Argument( m_DIM ) ;
137  m_argF = Argument( m_DIM + 1 ) ;
138 
139  }
140 
141 
151  NumericalDefiniteIntegral::NumericalDefiniteIntegral
152  ( const AbsFunction& function ,
153  const size_t index ,
154  const double a ,
155  const double b ,
156  const NumericalDefiniteIntegral::Points& points ,
157  const double epsabs ,
158  const double epsrel ,
159  const size_t size )
160  : m_function ( function.clone() )
161  , m_DIM ( 0 )
162  , m_index ( index )
163  , m_a ( a )
164  , m_b ( b )
165  , m_ia ( false )
166  , m_ib ( false )
167  , m_type ( GaudiMath::Integration:: Other )
168  , m_category ( GaudiMath::Integration:: Singular )
169  , m_rule ( GaudiMath::Integration:: Fixed )
170  , m_points ( points )
171  , m_epsabs ( epsabs )
172  , m_epsrel ( epsrel )
173  //
174  , m_result ( GSL_NEGINF )
175  , m_error ( GSL_POSINF )
176  //
177  , m_size ( size )
178  {
179  if ( function.dimensionality() < 2 )
180  { Exception("::constructor: invalid dimensionality ") ; }
181  if ( m_index >= function.dimensionality() )
182  { Exception("::constructor: invalid variable index") ; }
183 
184  m_DIM = function.dimensionality() - 1 ;
185  m_argument = Argument( m_DIM ) ;
186  m_argF = Argument( m_DIM + 1 ) ;
187 
188  const double l1 = std::min ( a , b ) ;
189  const double l2 = std::max ( a , b ) ;
190  m_points.push_back ( l1 ) ;
191  m_points.push_back ( l2 ) ;
192  // TODO: replace by remove_if (<l1,>l2), unique, erase and sort...
193  std::sort ( m_points.begin() , m_points.end() ) ;
194  m_points.erase( std::unique( m_points.begin () , m_points.end () ) ,
195  m_points.end() );
196 
197  auto lower = std::lower_bound ( m_points.begin () , m_points.end () , l1 ) ;
198  m_points.erase( m_points.begin () , lower ) ;
199  auto upper = std::upper_bound ( m_points.begin () , m_points.end () , l2 ) ;
200  m_points.erase( upper , m_points.end () ) ;
201 
202  m_pdata.reset( new double[ m_points.size() ] );
203  std::copy( m_points.begin() , m_points.end() , m_pdata.get() );
204  }
205 
231  NumericalDefiniteIntegral::NumericalDefiniteIntegral
232  ( const AbsFunction& function ,
233  const size_t index ,
234  const double a ,
235  const GaudiMath::Integration::Inf /* b */ ,
236  const double epsabs ,
237  const double epsrel ,
238  const size_t size )
239  : m_function ( function.clone() )
240  , m_DIM ( 0 )
241  , m_index ( index )
242  , m_a ( a )
243  , m_b ( GSL_POSINF )
244  , m_ia ( false )
245  , m_ib ( true )
246  , m_type ( GaudiMath::Integration:: Other )
247  , m_category ( GaudiMath::Integration:: Infinite )
248  , m_rule ( GaudiMath::Integration:: Fixed )
249  , m_epsabs ( epsabs )
250  , m_epsrel ( epsrel )
251  //
252  , m_result ( GSL_NEGINF )
253  , m_error ( GSL_POSINF )
254  //
255  , m_size ( size )
256  {
257  if ( function.dimensionality() < 2 )
258  { Exception("::constructor: invalid dimensionality ") ; }
259  if ( m_index >= function.dimensionality() )
260  { Exception("::constructor: invalid variable index") ; }
261 
262  m_DIM = function.dimensionality() - 1 ;
263  m_argument = Argument( m_DIM ) ;
264  m_argF = Argument( m_DIM + 1 ) ;
265 
266  }
267 
268 
294  NumericalDefiniteIntegral::NumericalDefiniteIntegral
295  ( const AbsFunction& function ,
296  const size_t index ,
297  const GaudiMath::Integration::Inf /* a */ ,
298  const double b ,
299  const double epsabs ,
300  const double epsrel ,
301  const size_t size )
302  : m_function ( function.clone() )
303  , m_DIM ( 0 )
304  , m_index ( index )
305  , m_a ( GSL_NEGINF )
306  , m_b ( b )
307  , m_ia ( true )
308  , m_ib ( false )
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  //
315  , m_result ( GSL_NEGINF )
316  , m_error ( GSL_POSINF )
317  //
318  , m_size ( size )
319  {
320  if ( function.dimensionality() < 2 )
321  { Exception("::constructor: invalid dimensionality ") ; }
322  if ( m_index >= function.dimensionality() )
323  { Exception("::constructor: invalid variable index") ; }
324 
325  m_DIM = function.dimensionality() - 1 ;
326  m_argument = Argument( m_DIM ) ;
327  m_argF = Argument( m_DIM + 1 ) ;
328  }
329 
353  NumericalDefiniteIntegral::NumericalDefiniteIntegral
354  ( const AbsFunction& function ,
355  const size_t index ,
356  const float epsabs ,
357  const float epsrel ,
358  const size_t size )
359  : AbsFunction()
360  , m_function ( function.clone() )
361  , m_DIM ( 0 )
362  , m_index ( index )
363  , m_a ( GSL_NEGINF )
364  , m_b ( GSL_POSINF )
365  , m_ia ( true )
366  , m_ib ( true )
367  , m_type ( GaudiMath::Integration:: Other )
368  , m_category ( GaudiMath::Integration:: Infinite )
369  , m_rule ( GaudiMath::Integration:: Fixed )
370  , m_epsabs ( epsabs )
371  , m_epsrel ( epsrel )
372  //
373  , m_result ( GSL_NEGINF )
374  , m_error ( GSL_POSINF )
375  //
376  , m_size ( size )
377  {
378  if ( function.dimensionality() < 2 )
379  { Exception("::constructor: invalid dimensionality ") ; }
380  if ( m_index >= function.dimensionality() )
381  { Exception("::constructor: invalid variable index") ; }
382 
383  m_DIM = function.dimensionality() - 1 ;
384  m_argument = Argument( m_DIM ) ;
385  m_argF = Argument( m_DIM + 1 ) ;
386 
387  }
388 
390  NumericalDefiniteIntegral::NumericalDefiniteIntegral
391  ( const NumericalDefiniteIntegral& right )
392  : AbsFunction ()
393  , m_function ( right.m_function->clone() )
394  , m_DIM ( right.m_DIM )
395  , m_index ( right.m_index )
396  , m_a ( right.m_a )
397  , m_b ( right.m_b )
398  , m_ia ( right.m_ia )
399  , m_ib ( right.m_ib )
400  , m_type ( right.m_type )
401  , m_category ( right.m_category )
402  , m_rule ( right.m_rule )
403  , m_points ( right.m_points )
404  , m_epsabs ( right.m_epsabs )
405  , m_epsrel ( right.m_epsrel )
406  , m_result ( GSL_NEGINF )
407  , m_error ( GSL_POSINF )
408  , m_size ( right.m_size )
409  , m_argument ( right.m_argument )
410  , m_argF ( right.m_argF )
411  {
412  if( right.m_pdata )
413  {
414  m_pdata.reset( new double[m_points.size()] );
415  std::copy( m_points.begin() , m_points.end() , m_pdata.get() );
416  }
417  }
418 
419  // ========================================================================
420  // throw the exception
421  // ========================================================================
422  StatusCode NumericalDefiniteIntegral::Exception
423  ( const std::string& message ,
424  const StatusCode& sc ) const
425  {
426  throw GaudiException( "NumericalDefiniteIntegral::" + message ,
427  "*GaudiMath*" , sc ) ;
428  return sc ;
429  }
430  // ========================================================================
431 
432  // ========================================================================
434  // ========================================================================
435  double NumericalDefiniteIntegral::operator()
436  ( const double argument ) const
437  {
438  // reset the result and the error
439  m_result = GSL_NEGINF ;
440  m_error = GSL_POSINF ;
441 
442  // check the argument
443  if( 1 != m_DIM ) { Exception ( "operator(): invalid argument size " ) ; };
444 
445  m_argument[0] = argument ;
446  return (*this) ( m_argument );
447  }
448  // ========================================================================
449 
450  // ========================================================================
452  // ========================================================================
453  Genfun::Derivative
454  NumericalDefiniteIntegral::partial ( unsigned int idx ) const
455  {
456  if ( idx >= m_DIM )
457  { Exception ( "::partial(i): invalid variable index " ) ; };
458  //
459  const AbsFunction& aux = NumericalDerivative( *this , idx ) ;
460  return Genfun::FunctionNoop( &aux ) ;
461  }
462  // ========================================================================
463 
464 
465  // ========================================================================
467  // ========================================================================
468  double NumericalDefiniteIntegral::operator()
469  ( const Argument& argument ) const
470  {
471  // reset the result and the error
472  m_result = GSL_NEGINF ;
473  m_error = GSL_POSINF ;
474 
475  // check the argument
476  if( argument.dimension() != m_DIM )
477  { Exception ( "operator(): invalid argument size " ) ; };
478 
479  // copy the argument
480 
481  for( size_t i = 0 ; i < m_DIM ; ++i )
482  {
483  m_argument [i] = argument [i] ;
484  const size_t j = i < m_index ? i : i + 1 ;
485  m_argF [j] = argument [i] ;
486  };
487 
488  // create the helper object
489  GSL_Helper helper( *m_function , m_argF , m_index );
490 
491  // use GSL to evaluate the numerical derivative
492  gsl_function F ;
493  F.function = &GSL_Adaptor ;
494  F.params = &helper ;
495  _Function F1 ;
496  F1.fn = &F ;
497 
498  if ( GaudiMath::Integration::Infinite == category () )
499  { return QAGI ( &F1 ) ; } // RETURN
500 
501  if ( m_a == m_b )
502  {
503  m_result = 0 ; m_error = 0 ; // EXACT
504  return m_result ; // RETURN
505  }
506 
507  if ( GaudiMath::Integration::Singular == category () )
508  { return QAGP ( &F1 ) ; } // RETURN
509  else if ( GaudiMath::Integration::Finite == category () )
510  if ( GaudiMath::Integration::NonAdaptive == type () )
511  { return QNG ( &F1 ) ; } // RETURN
512  else if ( GaudiMath::Integration::Adaptive == type () )
513  { return QAG ( &F1 ) ; } // RETURN
514  else if ( GaudiMath::Integration::AdaptiveSingular == type () )
515  { return QAGS ( &F1 ) ; } // RETURN
516  else
517  { Exception ( "::operator(): invalid type " ); }
518  else
519  { Exception ( "::operator(): invalid category " ); }
520 
521  return 0 ;
522  }
523  // ========================================================================
524 
525  // ========================================================================
527  // ========================================================================
528  NumericalDefiniteIntegral::_Workspace*
529  NumericalDefiniteIntegral::allocate() const
530  {
531  if ( !m_ws ) {
532  m_ws.reset( new _Workspace() );
533  m_ws->ws = gsl_integration_workspace_alloc( size () );
534  if ( !m_ws->ws ) { Exception ( "allocate()::invalid workspace" ) ; };
535  }
536  return m_ws.get() ;
537  }
538  // ========================================================================
539 
540  // ========================================================================
541  // adaptive integration on infinite interval
542  // ========================================================================
543  double NumericalDefiniteIntegral::QAGI ( _Function* F ) const
544  {
545  // check the argument
546  if ( !F ) { Exception("::QAGI: invalid function"); }
547 
548  // allocate workspace
549  allocate() ;
550 
551  int ierror = 0 ;
552 
553  if ( m_ia && m_ib )
554  {
555  ierror = gsl_integration_qagi ( F->fn ,
556  m_epsabs , m_epsrel ,
557  size () , ws()->ws ,
558  &m_result , &m_error ) ;
559  }
560  else if ( m_ia )
561  {
562  ierror = gsl_integration_qagil ( F->fn , m_b ,
563  m_epsabs , m_epsrel ,
564  size () , ws()->ws ,
565  &m_result , &m_error ) ;
566  }
567  else if ( m_ib )
568  {
569  ierror = gsl_integration_qagiu ( F->fn , m_a ,
570  m_epsabs , m_epsrel ,
571  size () , ws()->ws ,
572  &m_result , &m_error ) ;
573  }
574  else
575  { Exception ( "::QAGI: invalid mode" ) ; };
576 
577  if( ierror ) { gsl_error( "NumericalDefiniteIntegral::QAGI" ,
578  __FILE__ , __LINE__ , ierror ) ;}
579 
580  return m_result ;
581  }
582  // ========================================================================
583 
584  // ========================================================================
585  // adaptive integration with known singular points
586  // ========================================================================
587  double NumericalDefiniteIntegral::QAGP( _Function* F ) const
588  {
589  if( !F ) { Exception("QAGP::invalid function") ; }
590 
591  // no known singular points ?
592  if( points().empty() || !m_pdata ) { return QAGS( F ) ; }
593 
594  const size_t npts = points().size();
595 
596  // use GSL
597  int ierror =
598  gsl_integration_qagp ( F->fn ,
599  m_pdata.get() , npts ,
600  m_epsabs , m_epsrel ,
601  size () , ws()->ws ,
602  &m_result , &m_error ) ;
603 
604  if( ierror ) { gsl_error( "NumericalDefiniteIntegral::QAGI " ,
605  __FILE__ , __LINE__ , ierror ) ; }
606 
607  // the sign
608  if ( m_a > m_b ) { m_result *= -1 ; }
609 
610  return m_result ;
611  }
612  // ========================================================================
613 
614  // ========================================================================
615  // non-adaptive integration
616  // ========================================================================
617  double NumericalDefiniteIntegral::QNG ( _Function* F ) const
618  {
619  if( !F ) { Exception("QNG::invalid function") ; }
620 
621  // integration limits
622  const double low = std::min ( m_a , m_b ) ;
623  const double high = std::max ( m_a , m_b ) ;
624 
625  size_t neval = 0 ;
626  int ierror =
627  gsl_integration_qng ( F->fn ,
628  low , high ,
629  m_epsabs , m_epsrel ,
630  &m_result , &m_error , &neval ) ;
631 
632  if( ierror ) { gsl_error( "NumericalIndefiniteIntegral::QNG " ,
633  __FILE__ , __LINE__ , ierror ) ; }
634 
635  // sign
636  if ( m_a > m_b ) { m_result *= -1 ; }
637 
638  return m_result ;
639  }
640  // ========================================================================
641 
642 
643  // ========================================================================
644  // simple adaptive integration
645  // ========================================================================
646  double NumericalDefiniteIntegral::QAG ( _Function* F ) const
647  {
648  if( !F ) { Exception("QAG::invalid function") ; }
649 
650  // allocate workspace
651  allocate () ;
652 
653  // integration limits
654  const double low = std::min ( m_a , m_b ) ;
655  const double high = std::max ( m_a , m_b ) ;
656 
657  int ierror =
658  gsl_integration_qag ( F->fn ,
659  low , high ,
660  m_epsabs , m_epsrel ,
661  size () , (int) rule() , ws ()->ws ,
662  &m_result , &m_error );
663 
664  if( ierror ) { gsl_error( "NumericalDefiniteIntegral::QAG " ,
665  __FILE__ , __LINE__ , ierror ) ; }
666 
667  // sign
668  if ( m_a > m_b ) { m_result *= -1 ; }
669 
670  return m_result ;
671  }
672  // ========================================================================
673 
674  // ========================================================================
675  // adaptive integration with singularities
676  // ========================================================================
677  double NumericalDefiniteIntegral::QAGS ( _Function* F ) const
678  {
679  if( !F ) { Exception("QAG::invalid function") ; }
680 
681  if ( m_a == m_b )
682  {
683  m_result = 0 ;
684  m_error = 0 ; // EXACT !
685  return m_result ;
686  }
687 
688  // allocate workspace
689  allocate () ;
690 
691  // integration limits
692  const double low = std::min ( m_a , m_b ) ;
693  const double high = std::max ( m_a , m_b ) ;
694 
695  int ierror =
696  gsl_integration_qags ( F->fn ,
697  low , high ,
698  m_epsabs , m_epsrel ,
699  size () , ws()->ws ,
700  &m_result , &m_error );
701 
702  if( ierror ) { gsl_error( "NumericalDefiniteIntegral::QAGS " ,
703  __FILE__ , __LINE__ , ierror ) ; }
704 
705  // sign
706  if ( m_a > m_b ) { m_result *= -1 ; }
707 
708  return m_result ;
709  }
710  // ========================================================================
711 
712 
713  } // end of namespace GaudiMathImplementation
714 } // end of namespace Genfun
715 
716 
717 
718 // ============================================================================
719 // The END
720 // ============================================================================
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:423
std::copy
T copy(T...args)
GaudiException
Define general base for Gaudi exception.
Definition: GaudiException.h:23
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::_Workspace
Definition: NumericalDefiniteIntegral.h:82
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::_Workspace::ws
gsl_integration_workspace * ws
Definition: NumericalDefiniteIntegral.h:82
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::partial
Genfun::Derivative partial(unsigned int index) const override
Derivatives.
Definition: NumericalDefiniteIntegral.cpp:454
std::upper_bound
T upper_bound(T...args)
std::vector::end
T end(T...args)
GaudiMath::Derivative
Genfun::GaudiMathImplementation::NumericalDerivative Derivative
Definition: GaudiMath.h:29
std::lower_bound
T lower_bound(T...args)
std::unique
T unique(T...args)
std::string
STL class.
std::min
T min(T...args)
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)
Genfun::GaudiMathImplementation::GSL_Adaptor
double GSL_Adaptor(double x, void *params)
Definition: Helpers.cpp:37
std::unique_ptr::reset
T reset(T...args)
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::_Function
Definition: NumericalDefiniteIntegral.cpp:46
std::max
T max(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:529
std::unique_ptr::get
T get(T...args)
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::type
GaudiMath::Integration::Type type() const
integration type
Definition: NumericalDefiniteIntegral.h:348
std::vector::size
T size(T...args)
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