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  // ========================================================================
51  // ========================================================================
52  FUNCTION_OBJECT_IMP( NumericalDefiniteIntegral )
53  // ========================================================================
54 
55 
106  NumericalDefiniteIntegral::NumericalDefiniteIntegral
107  ( const AbsFunction& function ,
108  const size_t index ,
109  const double a ,
110  const double b ,
111  const GaudiMath::Integration::Type type ,
112  const GaudiMath::Integration::KronrodRule rule ,
113  const double epsabs ,
114  const double epsrel ,
115  const size_t size )
116  : AbsFunction ()
117  , m_function ( function.clone () )
118  , m_DIM ( 0 )
119  , m_index ( index )
120  , m_a ( a )
121  , m_b ( b )
122  , m_ia ( false )
123  , m_ib ( false )
124  , m_type ( type )
125  , m_category ( GaudiMath::Integration::Finite )
126  , m_rule ( rule )
127  , m_points ( )
128  , m_epsabs ( epsabs )
129  , m_epsrel ( epsrel )
130  , m_result ( GSL_NEGINF )
131  , m_error ( GSL_POSINF )
132  , m_size ( size )
133  , m_argument ()
134  , m_argF ()
135  {
136  if ( GaudiMath::Integration::Fixed == m_rule )
137  { m_rule = GaudiMath::Integration::Default ; }
138  if ( function.dimensionality() < 2 )
139  { Exception("::constructor: invalid dimensionality ") ; }
140  if ( m_index >= function.dimensionality() )
141  { Exception("::constructor: invalid variable index") ; }
142 
143  m_DIM = function.dimensionality() - 1 ;
144  m_argument = Argument( m_DIM ) ;
145  m_argF = Argument( m_DIM + 1 ) ;
146 
147  }
148 
149 
159  NumericalDefiniteIntegral::NumericalDefiniteIntegral
160  ( const AbsFunction& function ,
161  const size_t index ,
162  const double a ,
163  const double b ,
164  const NumericalDefiniteIntegral::Points& points ,
165  const double epsabs ,
166  const double epsrel ,
167  const size_t size )
168  : AbsFunction ()
169  , m_function ( function.clone() )
170  , m_DIM ( 0 )
171  , m_index ( index )
172  , m_a ( a )
173  , m_b ( b )
174  , m_ia ( false )
175  , m_ib ( false )
176  , m_type ( GaudiMath::Integration:: Other )
177  , m_category ( GaudiMath::Integration:: Singular )
178  , m_rule ( GaudiMath::Integration:: Fixed )
179  , m_points ( points )
180  , m_epsabs ( epsabs )
181  , m_epsrel ( epsrel )
182  //
183  , m_result ( GSL_NEGINF )
184  , m_error ( GSL_POSINF )
185  //
186  , m_size ( size )
187  , m_argument ()
188  , m_argF ()
189  {
190  if ( function.dimensionality() < 2 )
191  { Exception("::constructor: invalid dimensionality ") ; }
192  if ( m_index >= function.dimensionality() )
193  { Exception("::constructor: invalid variable index") ; }
194 
195  m_DIM = function.dimensionality() - 1 ;
196  m_argument = Argument( m_DIM ) ;
197  m_argF = Argument( m_DIM + 1 ) ;
198 
199  const double l1 = std::min ( a , b ) ;
200  const double l2 = std::max ( a , b ) ;
201  m_points.push_back ( l1 ) ;
202  m_points.push_back ( l2 ) ;
203  std::sort ( m_points.begin() , m_points.end() ) ;
204  m_points.erase( std::unique( m_points.begin () ,
205  m_points.end () ) ,
206  m_points.end() );
207 
208  Points::iterator lower =
209  std::lower_bound ( m_points.begin () , m_points.end () , l1 ) ;
210  m_points.erase ( m_points.begin () , lower ) ;
211  Points::iterator upper =
212  std::upper_bound ( m_points.begin () , m_points.end () , l2 ) ;
213  m_points.erase ( upper , m_points.end () ) ;
214 
215  m_pdata.reset( new double[ m_points.size() ] );
216  std::copy( m_points.begin() , m_points.end() , m_pdata.get() );
217  }
218 
244  NumericalDefiniteIntegral::NumericalDefiniteIntegral
245  ( const AbsFunction& function ,
246  const size_t index ,
247  const double a ,
248  const GaudiMath::Integration::Inf /* b */ ,
249  const double epsabs ,
250  const double epsrel ,
251  const size_t size )
252  : AbsFunction()
253  , m_function ( function.clone() )
254  , m_DIM ( 0 )
255  , m_index ( index )
256  , m_a ( a )
257  , m_b ( GSL_POSINF )
258  , m_ia ( false )
259  , m_ib ( true )
260  , m_type ( GaudiMath::Integration:: Other )
261  , m_category ( GaudiMath::Integration:: Infinite )
262  , m_rule ( GaudiMath::Integration:: Fixed )
263  , m_points ( )
264  , m_epsabs ( epsabs )
265  , m_epsrel ( epsrel )
266  //
267  , m_result ( GSL_NEGINF )
268  , m_error ( GSL_POSINF )
269  //
270  , m_size ( size )
271  , m_argument ()
272  , m_argF ()
273  {
274  if ( function.dimensionality() < 2 )
275  { Exception("::constructor: invalid dimensionality ") ; }
276  if ( m_index >= function.dimensionality() )
277  { Exception("::constructor: invalid variable index") ; }
278 
279  m_DIM = function.dimensionality() - 1 ;
280  m_argument = Argument( m_DIM ) ;
281  m_argF = Argument( m_DIM + 1 ) ;
282 
283  }
284 
285 
311  NumericalDefiniteIntegral::NumericalDefiniteIntegral
312  ( const AbsFunction& function ,
313  const size_t index ,
314  const GaudiMath::Integration::Inf /* a */ ,
315  const double b ,
316  const double epsabs ,
317  const double epsrel ,
318  const size_t size )
319  : AbsFunction()
320  , m_function ( function.clone() )
321  , m_DIM ( 0 )
322  , m_index ( index )
323  , m_a ( GSL_NEGINF )
324  , m_b ( b )
325  , m_ia ( true )
326  , m_ib ( false )
327  , m_type ( GaudiMath::Integration:: Other )
328  , m_category ( GaudiMath::Integration:: Infinite )
329  , m_rule ( GaudiMath::Integration:: Fixed )
330  , m_points ( )
331  , m_epsabs ( epsabs )
332  , m_epsrel ( epsrel )
333  //
334  , m_result ( GSL_NEGINF )
335  , m_error ( GSL_POSINF )
336  //
337  , m_size ( size )
338  , m_argument ()
339  , m_argF ()
340  {
341  if ( function.dimensionality() < 2 )
342  { Exception("::constructor: invalid dimensionality ") ; }
343  if ( m_index >= function.dimensionality() )
344  { Exception("::constructor: invalid variable index") ; }
345 
346  m_DIM = function.dimensionality() - 1 ;
347  m_argument = Argument( m_DIM ) ;
348  m_argF = Argument( m_DIM + 1 ) ;
349  }
350 
374  NumericalDefiniteIntegral::NumericalDefiniteIntegral
375  ( const AbsFunction& function ,
376  const size_t index ,
377  const float epsabs ,
378  const float epsrel ,
379  const size_t size )
380  : AbsFunction()
381  , m_function ( function.clone() )
382  , m_DIM ( 0 )
383  , m_index ( index )
384  , m_a ( GSL_NEGINF )
385  , m_b ( GSL_POSINF )
386  , m_ia ( true )
387  , m_ib ( true )
388  , m_type ( GaudiMath::Integration:: Other )
389  , m_category ( GaudiMath::Integration:: Infinite )
390  , m_rule ( GaudiMath::Integration:: Fixed )
391  , m_points ( )
392  , m_epsabs ( epsabs )
393  , m_epsrel ( epsrel )
394  //
395  , m_result ( GSL_NEGINF )
396  , m_error ( GSL_POSINF )
397  //
398  , m_size ( size )
399  , m_argument ()
400  , m_argF ()
401  {
402  if ( function.dimensionality() < 2 )
403  { Exception("::constructor: invalid dimensionality ") ; }
404  if ( m_index >= function.dimensionality() )
405  { Exception("::constructor: invalid variable index") ; }
406 
407  m_DIM = function.dimensionality() - 1 ;
408  m_argument = Argument( m_DIM ) ;
409  m_argF = Argument( m_DIM + 1 ) ;
410 
411  }
412 
414  NumericalDefiniteIntegral::NumericalDefiniteIntegral
415  ( const NumericalDefiniteIntegral& right )
416  : AbsFunction ()
417  , m_function ( right.m_function->clone() )
418  , m_DIM ( right.m_DIM )
419  , m_index ( right.m_index )
420  , m_a ( right.m_a )
421  , m_b ( right.m_b )
422  , m_ia ( right.m_ia )
423  , m_ib ( right.m_ib )
424  , m_type ( right.m_type )
425  , m_category ( right.m_category )
426  , m_rule ( right.m_rule )
427  , m_points ( right.m_points )
428  , m_epsabs ( right.m_epsabs )
429  , m_epsrel ( right.m_epsrel )
430  , m_result ( GSL_NEGINF )
431  , m_error ( GSL_POSINF )
432  , m_size ( right.m_size )
433  , m_argument ( right.m_argument )
434  , m_argF ( right.m_argF )
435  {
436  if( right.m_pdata )
437  {
438  m_pdata.reset( new double[m_points.size()] );
439  std::copy( m_points.begin() , m_points.end() , m_pdata.get() );
440  }
441  }
442 
443  // ========================================================================
444  // throw the exception
445  // ========================================================================
446  StatusCode NumericalDefiniteIntegral::Exception
447  ( const std::string& message ,
448  const StatusCode& sc ) const
449  {
450  throw GaudiException( "NumericalDefiniteIntegral::" + message ,
451  "*GaudiMath*" , sc ) ;
452  return sc ;
453  }
454  // ========================================================================
455 
456  // ========================================================================
458  // ========================================================================
459  double NumericalDefiniteIntegral::operator()
460  ( const double argument ) const
461  {
462  // reset the result and the error
463  m_result = GSL_NEGINF ;
464  m_error = GSL_POSINF ;
465 
466  // check the argument
467  if( 1 != m_DIM ) { Exception ( "operator(): invalid argument size " ) ; };
468 
469  m_argument[0] = argument ;
470  return (*this) ( m_argument );
471  }
472  // ========================================================================
473 
474  // ========================================================================
476  // ========================================================================
477  Genfun::Derivative
478  NumericalDefiniteIntegral::partial ( unsigned int idx ) const
479  {
480  if ( idx >= m_DIM )
481  { Exception ( "::partial(i): invalid variable index " ) ; };
482  //
483  const AbsFunction& aux = NumericalDerivative( *this , idx ) ;
484  return Genfun::FunctionNoop( &aux ) ;
485  }
486  // ========================================================================
487 
488 
489  // ========================================================================
491  // ========================================================================
492  double NumericalDefiniteIntegral::operator()
493  ( const Argument& argument ) const
494  {
495  // reset the result and the error
496  m_result = GSL_NEGINF ;
497  m_error = GSL_POSINF ;
498 
499  // check the argument
500  if( argument.dimension() != m_DIM )
501  { Exception ( "operator(): invalid argument size " ) ; };
502 
503  // copy the argument
504 
505  for( size_t i = 0 ; i < m_DIM ; ++i )
506  {
507  m_argument [i] = argument [i] ;
508  const size_t j = i < m_index ? i : i + 1 ;
509  m_argF [j] = argument [i] ;
510  };
511 
512  // create the helper object
513  GSL_Helper helper( *m_function , m_argF , m_index );
514 
515  // use GSL to evaluate the numerical derivative
516  gsl_function F ;
517  F.function = &GSL_Adaptor ;
518  F.params = &helper ;
519  _Function F1 ;
520  F1.fn = &F ;
521 
522  if ( GaudiMath::Integration::Infinite == category () )
523  { return QAGI ( &F1 ) ; } // RETURN
524 
525  if ( m_a == m_b )
526  {
527  m_result = 0 ; m_error = 0 ; // EXACT
528  return m_result ; // RETURN
529  }
530 
531  if ( GaudiMath::Integration::Singular == category () )
532  { return QAGP ( &F1 ) ; } // RETURN
533  else if ( GaudiMath::Integration::Finite == category () )
534  if ( GaudiMath::Integration::NonAdaptive == type () )
535  { return QNG ( &F1 ) ; } // RETURN
536  else if ( GaudiMath::Integration::Adaptive == type () )
537  { return QAG ( &F1 ) ; } // RETURN
538  else if ( GaudiMath::Integration::AdaptiveSingular == type () )
539  { return QAGS ( &F1 ) ; } // RETURN
540  else
541  { Exception ( "::operator(): invalid type " ); }
542  else
543  { Exception ( "::operator(): invalid category " ); }
544 
545  return 0 ;
546  }
547  // ========================================================================
548 
549  // ========================================================================
551  // ========================================================================
552  NumericalDefiniteIntegral::_Workspace*
553  NumericalDefiniteIntegral::allocate() const
554  {
555  if ( !m_ws ) {
556  m_ws.reset( new _Workspace() );
557  m_ws->ws = gsl_integration_workspace_alloc( size () );
558  if ( !m_ws->ws ) { Exception ( "allocate()::invalid workspace" ) ; };
559  }
560  return m_ws.get() ;
561  }
562  // ========================================================================
563 
564  // ========================================================================
565  // adaptive integration on infinite interval
566  // ========================================================================
567  double NumericalDefiniteIntegral::QAGI ( _Function* F ) const
568  {
569  // check the argument
570  if ( 0 == F ) { Exception("::QAGI: invalid function"); }
571 
572  // allocate workspace
573  allocate() ;
574 
575  int ierror = 0 ;
576 
577  if ( m_ia && m_ib )
578  {
579  ierror = gsl_integration_qagi ( F->fn ,
580  m_epsabs , m_epsrel ,
581  size () , ws()->ws ,
582  &m_result , &m_error ) ;
583  }
584  else if ( m_ia )
585  {
586  ierror = gsl_integration_qagil ( F->fn , m_b ,
587  m_epsabs , m_epsrel ,
588  size () , ws()->ws ,
589  &m_result , &m_error ) ;
590  }
591  else if ( m_ib )
592  {
593  ierror = gsl_integration_qagiu ( F->fn , m_a ,
594  m_epsabs , m_epsrel ,
595  size () , ws()->ws ,
596  &m_result , &m_error ) ;
597  }
598  else
599  { Exception ( "::QAGI: invalid mode" ) ; };
600 
601  if( ierror ) { gsl_error( "NumericalDefiniteIntegral::QAGI" ,
602  __FILE__ , __LINE__ , ierror ) ;}
603 
604  return m_result ;
605  }
606  // ========================================================================
607 
608  // ========================================================================
609  // adaptive integration with known singular points
610  // ========================================================================
611  double NumericalDefiniteIntegral::QAGP( _Function* F ) const
612  {
613  if( 0 == F ) { Exception("QAGP::invalid function") ; }
614 
615  // no known singular points ?
616  if( points().empty() || !m_pdata ) { return QAGS( F ) ; }
617 
618  const size_t npts = points().size();
619 
620  // use GSL
621  int ierror =
622  gsl_integration_qagp ( F->fn ,
623  m_pdata.get() , npts ,
624  m_epsabs , m_epsrel ,
625  size () , ws()->ws ,
626  &m_result , &m_error ) ;
627 
628  if( ierror ) { gsl_error( "NumericalDefiniteIntegral::QAGI " ,
629  __FILE__ , __LINE__ , ierror ) ; }
630 
631  // the sign
632  if ( m_a > m_b ) { m_result *= -1 ; }
633 
634  return m_result ;
635  }
636  // ========================================================================
637 
638  // ========================================================================
639  // non-adaptive integration
640  // ========================================================================
641  double NumericalDefiniteIntegral::QNG ( _Function* F ) const
642  {
643  if( 0 == F ) { Exception("QNG::invalid function") ; }
644 
645  // integration limits
646  const double low = std::min ( m_a , m_b ) ;
647  const double high = std::max ( m_a , m_b ) ;
648 
649  size_t neval = 0 ;
650  int ierror =
651  gsl_integration_qng ( F->fn ,
652  low , high ,
653  m_epsabs , m_epsrel ,
654  &m_result , &m_error , &neval ) ;
655 
656  if( ierror ) { gsl_error( "NumericalIndefiniteIntegral::QNG " ,
657  __FILE__ , __LINE__ , ierror ) ; }
658 
659  // sign
660  if ( m_a > m_b ) { m_result *= -1 ; }
661 
662  return m_result ;
663  }
664  // ========================================================================
665 
666 
667  // ========================================================================
668  // simple adaptive integration
669  // ========================================================================
670  double NumericalDefiniteIntegral::QAG ( _Function* F ) const
671  {
672  if( 0 == F ) { Exception("QAG::invalid function") ; }
673 
674  // allocate workspace
675  allocate () ;
676 
677  // integration limits
678  const double low = std::min ( m_a , m_b ) ;
679  const double high = std::max ( m_a , m_b ) ;
680 
681  int ierror =
682  gsl_integration_qag ( F->fn ,
683  low , high ,
684  m_epsabs , m_epsrel ,
685  size () , (int) rule() , ws ()->ws ,
686  &m_result , &m_error );
687 
688  if( ierror ) { gsl_error( "NumericalDefiniteIntegral::QAG " ,
689  __FILE__ , __LINE__ , ierror ) ; }
690 
691  // sign
692  if ( m_a > m_b ) { m_result *= -1 ; }
693 
694  return m_result ;
695  }
696  // ========================================================================
697 
698  // ========================================================================
699  // adaptive integration with singularities
700  // ========================================================================
701  double NumericalDefiniteIntegral::QAGS ( _Function* F ) const
702  {
703  if( 0 == F ) { Exception("QAG::invalid function") ; }
704 
705  if ( m_a == m_b )
706  {
707  m_result = 0 ;
708  m_error = 0 ; // EXACT !
709  return m_result ;
710  }
711 
712  // allocate workspace
713  allocate () ;
714 
715  // integration limits
716  const double low = std::min ( m_a , m_b ) ;
717  const double high = std::max ( m_a , m_b ) ;
718 
719  int ierror =
720  gsl_integration_qags ( F->fn ,
721  low , high ,
722  m_epsabs , m_epsrel ,
723  size () , ws()->ws ,
724  &m_result , &m_error );
725 
726  if( ierror ) { gsl_error( "NumericalDefiniteIntegral::QAGS " ,
727  __FILE__ , __LINE__ , ierror ) ; }
728 
729  // sign
730  if ( m_a > m_b ) { m_result *= -1 ; }
731 
732  return m_result ;
733  }
734  // ========================================================================
735 
736 
737  } // end of namespace GaudiMathImplementation
738 } // end of namespace Genfun
739 
740 
741 
742 // ============================================================================
743 // The END
744 // ============================================================================
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::Exception
StatusCode Exception(const std::string &message, const StatusCode &sc=StatusCode::FAILURE) const
Definition: NumericalDefiniteIntegral.cpp:447
GaudiException
Define general base for Gaudi exception.
Definition: GaudiException.h:23
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::_Workspace
Definition: NumericalDefiniteIntegral.h:74
GaudiMath::Derivative
Genfun::GaudiMathImplementation::NumericalDerivative Derivative
Definition: GaudiMath.h:29
merge_files.sc
int sc
Definition: merge_files.py:167
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::_Function::fn
gsl_function * fn
Definition: NumericalDefiniteIntegral.cpp:47
StatusCode
This class is used for returning status codes from appropriate routines.
Definition: StatusCode.h:26
Genfun::GaudiMathImplementation::GSL_Adaptor
double GSL_Adaptor(double x, void *params)
Definition: Helpers.cpp:37
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::_Function
Definition: NumericalDefiniteIntegral.cpp:46
Genfun::GaudiMathImplementation::NumericalDerivative
Numerical derivative (using GSL adaptive numerical differentiation)
Definition: NumericalDerivative.h:34
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::allocate
_Workspace * allocate() const
allocate the integration workspace
Definition: NumericalDefiniteIntegral.cpp:553
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::Points
std::vector< double > Points
typedef for vector of singular points
Definition: NumericalDefiniteIntegral.h:75
Helpers.h
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::partial
virtual Genfun::Derivative partial(unsigned int index) const
Derivatives.
Definition: NumericalDefiniteIntegral.cpp:478
GaudiMath
GaudiMath.h GaudiMath/GaudiMath.h.
Definition: Adapters.h:13
GaudiMath::Integration::KronrodRule
KronrodRule
integration rule
Definition: Integration.h:34
Genfun
CLHEP.
Definition: IEqSolver.h:13
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::_Workspace::ws
gsl_integration_workspace * ws
Definition: NumericalDefiniteIntegral.h:74
Genfun::GaudiMathImplementation::GSL_Helper
the simple structure to be used for adaption interface Genfun::AbsFunction to gsl_function structure ...
Definition: Helpers.h:26
ana.i
list i
Definition: ana.py:128
Tuples::Type
Type
the list of available types for ntuples
Definition: TupleObj.h:73
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::rule
GaudiMath::Integration::KronrodRule rule() const
integration rule
Definition: NumericalDefiniteIntegral.h:346
gaudirun.type
string type
Definition: gaudirun.py:151
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:71