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