•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
NumericalDefiniteIntegral.cpp
Go to the documentation of this file.
1 // $Id: NumericalDefiniteIntegral.cpp,v 1.4 2007/11/20 13:00:17 marcocle Exp $
2 // ============================================================================
3 // Include
4 // ============================================================================
5 // STD & STL
6 // ============================================================================
7 #include <vector>
8 #include <algorithm>
9 // ============================================================================
10 // GSL
11 // ============================================================================
12 #include "gsl/gsl_errno.h"
13 #include "gsl/gsl_integration.h"
14 // ============================================================================
15 // GaudiKernel
16 // ============================================================================
17 #include "GaudiKernel/GaudiException.h"
18 // ============================================================================
19 // GaudiMath
20 // ============================================================================
21 #include "GaudiMath/NumericalDefiniteIntegral.h"
22 #include "GaudiMath/NumericalDerivative.h"
23 // ============================================================================
24 // Local
25 // ============================================================================
26 #include "Helpers.h"
27 // ============================================================================
28 
29 
30 #ifdef __ICC
31 // disable icc remark #1572: floating-point equality and inequality comparisons are unreliable
32 // The comparison are meant
33 #pragma warning(disable:1572)
34 #endif
35 #ifdef WIN32
36 // Disable the warning
37 // C4996: 'std::copy': Function call with parameters that may be unsafe
38 // The parameters are checked
39 #pragma warning(disable:4996)
40 #endif
41 
42 namespace Genfun
43 {
44  namespace GaudiMathImplementation
45  {
46 
47  struct NumericalDefiniteIntegral::_Workspace
48  { gsl_integration_workspace* ws ; };
49  struct NumericalDefiniteIntegral::_Function
50  { gsl_function* fn ; };
51 
52  // ========================================================================
54  // ========================================================================
55  FUNCTION_OBJECT_IMP( NumericalDefiniteIntegral )
56  // ========================================================================
57 
58 
109  NumericalDefiniteIntegral::NumericalDefiniteIntegral
110  ( const AbsFunction& function ,
111  const size_t index ,
112  const double a ,
113  const double b ,
114  const GaudiMath::Integration::Type type ,
115  const GaudiMath::Integration::KronrodRule rule ,
116  const double epsabs ,
117  const double epsrel ,
118  const size_t size )
119  : AbsFunction ()
120  , m_function ( function.clone () )
121  , m_DIM ( 0 )
122  , m_index ( index )
123  , m_a ( a )
124  , m_b ( b )
125  , m_ia ( false )
126  , m_ib ( false )
127  , m_type ( type )
128  , m_category ( GaudiMath::Integration::Finite )
129  , m_rule ( rule )
130  , m_points ( )
131  , m_pdata ( 0 )
132  , m_epsabs ( epsabs )
133  , m_epsrel ( epsrel )
134  , m_result ( GSL_NEGINF )
135  , m_error ( GSL_POSINF )
136  , m_size ( size )
137  , m_ws ()
138  , m_argument ()
139  , m_argF ()
140  {
141  if ( GaudiMath::Integration::Fixed == m_rule )
142  { m_rule = GaudiMath::Integration::Default ; }
143  if ( function.dimensionality() < 2 )
144  { Exception("::constructor: invalid dimensionality ") ; }
145  if ( m_index >= function.dimensionality() )
146  { Exception("::constructor: invalid variable index") ; }
147 
148  m_DIM = function.dimensionality() - 1 ;
149  m_argument = Argument( m_DIM ) ;
150  m_argF = Argument( m_DIM + 1 ) ;
151 
152  }
153 
154 
164  NumericalDefiniteIntegral::NumericalDefiniteIntegral
165  ( const AbsFunction& function ,
166  const size_t index ,
167  const double a ,
168  const double b ,
169  const NumericalDefiniteIntegral::Points& points ,
170  const double epsabs ,
171  const double epsrel ,
172  const size_t size )
173  : AbsFunction ()
174  , m_function ( function.clone() )
175  , m_DIM ( 0 )
176  , m_index ( index )
177  , m_a ( a )
178  , m_b ( b )
179  , m_ia ( false )
180  , m_ib ( false )
181  , m_type ( GaudiMath::Integration:: Other )
182  , m_category ( GaudiMath::Integration:: Singular )
183  , m_rule ( GaudiMath::Integration:: Fixed )
184  , m_points ( points )
185  , m_pdata ( 0 )
186  , m_epsabs ( epsabs )
187  , m_epsrel ( epsrel )
188  //
189  , m_result ( GSL_NEGINF )
190  , m_error ( GSL_POSINF )
191  //
192  , m_size ( size )
193  , m_ws ( 0 )
194  , m_argument ()
195  , m_argF ()
196  {
197  if ( function.dimensionality() < 2 )
198  { Exception("::constructor: invalid dimensionality ") ; }
199  if ( m_index >= function.dimensionality() )
200  { Exception("::constructor: invalid variable index") ; }
201 
202  m_DIM = function.dimensionality() - 1 ;
203  m_argument = Argument( m_DIM ) ;
204  m_argF = Argument( m_DIM + 1 ) ;
205 
206  const double l1 = std::min ( a , b ) ;
207  const double l2 = std::max ( a , b ) ;
208  m_points.push_back ( l1 ) ;
209  m_points.push_back ( l2 ) ;
210  std::sort ( m_points.begin() , m_points.end() ) ;
211  m_points.erase( std::unique( m_points.begin () ,
212  m_points.end () ) ,
213  m_points.end() );
214 
215  Points::iterator lower =
216  std::lower_bound ( m_points.begin () , m_points.end () , l1 ) ;
217  m_points.erase ( m_points.begin () , lower ) ;
218  Points::iterator upper =
219  std::upper_bound ( m_points.begin () , m_points.end () , l2 ) ;
220  m_points.erase ( upper , m_points.end () ) ;
221 
222  m_pdata = new double[ m_points.size() ] ;
223  std::copy( m_points.begin() , m_points.end() , m_pdata );
224  }
225 
251  NumericalDefiniteIntegral::NumericalDefiniteIntegral
252  ( const AbsFunction& function ,
253  const size_t index ,
254  const double a ,
255  const GaudiMath::Integration::Inf /* b */ ,
256  const double epsabs ,
257  const double epsrel ,
258  const size_t size )
259  : AbsFunction()
260  , m_function ( function.clone() )
261  , m_DIM ( 0 )
262  , m_index ( index )
263  , m_a ( a )
264  , m_b ( GSL_POSINF )
265  , m_ia ( false )
266  , m_ib ( true )
267  , m_type ( GaudiMath::Integration:: Other )
268  , m_category ( GaudiMath::Integration:: Infinite )
269  , m_rule ( GaudiMath::Integration:: Fixed )
270  , m_points ( )
271  , m_pdata ( 0 )
272  , m_epsabs ( epsabs )
273  , m_epsrel ( epsrel )
274  //
275  , m_result ( GSL_NEGINF )
276  , m_error ( GSL_POSINF )
277  //
278  , m_size ( size )
279  , m_ws ( 0 )
280  , m_argument ()
281  , m_argF ()
282  {
283  if ( function.dimensionality() < 2 )
284  { Exception("::constructor: invalid dimensionality ") ; }
285  if ( m_index >= function.dimensionality() )
286  { Exception("::constructor: invalid variable index") ; }
287 
288  m_DIM = function.dimensionality() - 1 ;
289  m_argument = Argument( m_DIM ) ;
290  m_argF = Argument( m_DIM + 1 ) ;
291 
292  }
293 
294 
320  NumericalDefiniteIntegral::NumericalDefiniteIntegral
321  ( const AbsFunction& function ,
322  const size_t index ,
323  const GaudiMath::Integration::Inf /* a */ ,
324  const double b ,
325  const double epsabs ,
326  const double epsrel ,
327  const size_t size )
328  : AbsFunction()
329  , m_function ( function.clone() )
330  , m_DIM ( 0 )
331  , m_index ( index )
332  , m_a ( GSL_NEGINF )
333  , m_b ( b )
334  , m_ia ( true )
335  , m_ib ( false )
336  , m_type ( GaudiMath::Integration:: Other )
337  , m_category ( GaudiMath::Integration:: Infinite )
338  , m_rule ( GaudiMath::Integration:: Fixed )
339  , m_points ( )
340  , m_pdata ( 0 )
341  , m_epsabs ( epsabs )
342  , m_epsrel ( epsrel )
343  //
344  , m_result ( GSL_NEGINF )
345  , m_error ( GSL_POSINF )
346  //
347  , m_size ( size )
348  , m_ws ( 0 )
349  , m_argument ()
350  , m_argF ()
351  {
352  if ( function.dimensionality() < 2 )
353  { Exception("::constructor: invalid dimensionality ") ; }
354  if ( m_index >= function.dimensionality() )
355  { Exception("::constructor: invalid variable index") ; }
356 
357  m_DIM = function.dimensionality() - 1 ;
358  m_argument = Argument( m_DIM ) ;
359  m_argF = Argument( m_DIM + 1 ) ;
360  }
361 
385  NumericalDefiniteIntegral::NumericalDefiniteIntegral
386  ( const AbsFunction& function ,
387  const size_t index ,
388  const float epsabs ,
389  const float epsrel ,
390  const size_t size )
391  : AbsFunction()
392  , m_function ( function.clone() )
393  , m_DIM ( 0 )
394  , m_index ( index )
395  , m_a ( GSL_NEGINF )
396  , m_b ( GSL_POSINF )
397  , m_ia ( true )
398  , m_ib ( true )
399  , m_type ( GaudiMath::Integration:: Other )
400  , m_category ( GaudiMath::Integration:: Infinite )
401  , m_rule ( GaudiMath::Integration:: Fixed )
402  , m_points ( )
403  , m_pdata ( 0 )
404  , m_epsabs ( epsabs )
405  , m_epsrel ( epsrel )
406  //
407  , m_result ( GSL_NEGINF )
408  , m_error ( GSL_POSINF )
409  //
410  , m_size ( size )
411  , m_ws ( 0 )
412  , m_argument ()
413  , m_argF ()
414  {
415  if ( function.dimensionality() < 2 )
416  { Exception("::constructor: invalid dimensionality ") ; }
417  if ( m_index >= function.dimensionality() )
418  { Exception("::constructor: invalid variable index") ; }
419 
420  m_DIM = function.dimensionality() - 1 ;
421  m_argument = Argument( m_DIM ) ;
422  m_argF = Argument( m_DIM + 1 ) ;
423 
424  }
425 
427  NumericalDefiniteIntegral::NumericalDefiniteIntegral
428  ( const NumericalDefiniteIntegral& right )
429  : AbsFunction ()
430  , m_function ( right.m_function->clone() )
431  , m_DIM ( right.m_DIM )
432  , m_index ( right.m_index )
433  , m_a ( right.m_a )
434  , m_b ( right.m_b )
435  , m_ia ( right.m_ia )
436  , m_ib ( right.m_ib )
437  , m_type ( right.m_type )
438  , m_category ( right.m_category )
439  , m_rule ( right.m_rule )
440  , m_points ( right.m_points )
441  , m_pdata ( 0 )
442  , m_epsabs ( right.m_epsabs )
443  , m_epsrel ( right.m_epsrel )
444  , m_result ( GSL_NEGINF )
445  , m_error ( GSL_POSINF )
446  , m_size ( right.m_size )
447  , m_ws ( 0 )
448  , m_argument ( right.m_argument )
449  , m_argF ( right.m_argF )
450  {
451  if( 0 != right.m_pdata )
452  {
453  m_pdata = new double[m_points.size()] ;
454  std::copy( m_points.begin() , m_points.end() , m_pdata );
455  }
456  }
457 
458  NumericalDefiniteIntegral::~NumericalDefiniteIntegral()
459  {
460  if( 0 != m_function ) { delete m_function ; m_function = 0 ; }
461  if( 0 != m_pdata ) { delete[] m_pdata ; m_pdata = 0 ; }
462  }
463 
464  // ========================================================================
465  // throw the exception
466  // ========================================================================
467  StatusCode NumericalDefiniteIntegral::Exception
468  ( const std::string& message ,
469  const StatusCode& sc ) const
470  {
471  throw GaudiException( "NumericalDefiniteIntegral::" + message ,
472  "*GaudiMath*" , sc ) ;
473  return sc ;
474  }
475  // ========================================================================
476 
477  // ========================================================================
479  // ========================================================================
480  double NumericalDefiniteIntegral::operator()
481  ( const double argument ) const
482  {
483  // reset the result and the error
484  m_result = GSL_NEGINF ;
485  m_error = GSL_POSINF ;
486 
487  // check the argument
488  if( 1 != m_DIM ) { Exception ( "operator(): invalid argument size " ) ; };
489 
490  m_argument[0] = argument ;
491  return (*this) ( m_argument );
492  }
493  // ========================================================================
494 
495  // ========================================================================
497  // ========================================================================
498  Genfun::Derivative
499  NumericalDefiniteIntegral::partial ( unsigned int idx ) const
500  {
501  if ( idx >= m_DIM )
502  { Exception ( "::partial(i): invalid variable index " ) ; };
503  //
504  const AbsFunction& aux = NumericalDerivative( *this , idx ) ;
505  return Genfun::FunctionNoop( &aux ) ;
506  }
507  // ========================================================================
508 
509 
510  // ========================================================================
512  // ========================================================================
513  double NumericalDefiniteIntegral::operator()
514  ( const Argument& argument ) const
515  {
516  // reset the result and the error
517  m_result = GSL_NEGINF ;
518  m_error = GSL_POSINF ;
519 
520  // check the argument
521  if( argument.dimension() != m_DIM )
522  { Exception ( "operator(): invalid argument size " ) ; };
523 
524  // copy the argument
525 
526  for( size_t i = 0 ; i < m_DIM ; ++i )
527  {
528  m_argument [i] = argument [i] ;
529  const size_t j = i < m_index ? i : i + 1 ;
530  m_argF [j] = argument [i] ;
531  };
532 
533  // create the helper object
534  GSL_Helper helper( *m_function , m_argF , m_index );
535 
536  // use GSL to evaluate the numerical derivative
537  gsl_function F ;
538  F.function = &GSL_Adaptor ;
539  F.params = &helper ;
540  _Function F1 ;
541  F1.fn = &F ;
542 
543  if ( GaudiMath::Integration::Infinite == category () )
544  { return QAGI ( &F1 ) ; } // RETURN
545 
546  if ( m_a == m_b )
547  {
548  m_result = 0 ; m_error = 0 ; // EXACT
549  return m_result ; // RETURN
550  }
551 
552  if ( GaudiMath::Integration::Singular == category () )
553  { return QAGP ( &F1 ) ; } // RETURN
554  else if ( GaudiMath::Integration::Finite == category () )
555  if ( GaudiMath::Integration::NonAdaptive == type () )
556  { return QNG ( &F1 ) ; } // RETURN
557  else if ( GaudiMath::Integration::Adaptive == type () )
558  { return QAG ( &F1 ) ; } // RETURN
559  else if ( GaudiMath::Integration::AdaptiveSingular == type () )
560  { return QAGS ( &F1 ) ; } // RETURN
561  else
562  { Exception ( "::operator(): invalid type " ); }
563  else
564  { Exception ( "::operator(): invalid category " ); }
565 
566  return 0 ;
567  }
568  // ========================================================================
569 
570  // ========================================================================
572  // ========================================================================
573  NumericalDefiniteIntegral::_Workspace*
574  NumericalDefiniteIntegral::allocate() const
575  {
576  if ( 0 != m_ws ) { return m_ws; }
577  gsl_integration_workspace* aux =
578  gsl_integration_workspace_alloc( size () );
579  if ( 0 == aux ) { Exception ( "allocate()::invalid workspace" ) ; };
580  m_ws = new _Workspace() ;
581  m_ws->ws = aux ;
582  return m_ws ;
583  }
584  // ========================================================================
585 
586  // ========================================================================
587  // adaptive integration on infinite interval
588  // ========================================================================
589  double NumericalDefiniteIntegral::QAGI ( _Function* F ) const
590  {
591  // check the argument
592  if ( 0 == F ) { Exception("::QAGI: invalid function"); }
593 
594  // allocate workspace
595  if ( 0 == ws() ) { allocate() ; }
596 
597  int ierror = 0 ;
598 
599  if ( m_ia && m_ib )
600  {
601  ierror = gsl_integration_qagi ( F->fn ,
602  m_epsabs , m_epsrel ,
603  size () , ws()->ws ,
604  &m_result , &m_error ) ;
605  }
606  else if ( m_ia )
607  {
608  ierror = gsl_integration_qagil ( F->fn , m_b ,
609  m_epsabs , m_epsrel ,
610  size () , ws()->ws ,
611  &m_result , &m_error ) ;
612  }
613  else if ( m_ib )
614  {
615  ierror = gsl_integration_qagiu ( F->fn , m_a ,
616  m_epsabs , m_epsrel ,
617  size () , ws()->ws ,
618  &m_result , &m_error ) ;
619  }
620  else
621  { Exception ( "::QAGI: invalid mode" ) ; };
622 
623  if( ierror ) { gsl_error( "NumericalDefiniteIntegral::QAGI" ,
624  __FILE__ , __LINE__ , ierror ) ;}
625 
626  return m_result ;
627  }
628  // ========================================================================
629 
630  // ========================================================================
631  // adaptive integration with known singular points
632  // ========================================================================
633  double NumericalDefiniteIntegral::QAGP( _Function* F ) const
634  {
635  if( 0 == F ) { Exception("QAGP::invalid function") ; }
636 
637  // no known singular points ?
638  if( points().empty() || 0 == m_pdata ) { return QAGS( F ) ; }
639 
640  const size_t npts = points().size();
641 
642  // use GSL
643  int ierror =
644  gsl_integration_qagp ( F->fn ,
645  m_pdata , npts ,
646  m_epsabs , m_epsrel ,
647  size () , ws()->ws ,
648  &m_result , &m_error ) ;
649 
650  if( ierror ) { gsl_error( "NumericalDefiniteIntegral::QAGI " ,
651  __FILE__ , __LINE__ , ierror ) ; }
652 
653  // the sign
654  if ( m_a > m_b ) { m_result *= -1 ; }
655 
656  return m_result ;
657  }
658  // ========================================================================
659 
660  // ========================================================================
661  // non-adaptive integration
662  // ========================================================================
663  double NumericalDefiniteIntegral::QNG ( _Function* F ) const
664  {
665  if( 0 == F ) { Exception("QNG::invalid function") ; }
666 
667  // integration limits
668  const double low = std::min ( m_a , m_b ) ;
669  const double high = std::max ( m_a , m_b ) ;
670 
671  size_t neval = 0 ;
672  int ierror =
673  gsl_integration_qng ( F->fn ,
674  low , high ,
675  m_epsabs , m_epsrel ,
676  &m_result , &m_error , &neval ) ;
677 
678  if( ierror ) { gsl_error( "NumericalIndefiniteIntegral::QNG " ,
679  __FILE__ , __LINE__ , ierror ) ; }
680 
681  // sign
682  if ( m_a > m_b ) { m_result *= -1 ; }
683 
684  return m_result ;
685  }
686  // ========================================================================
687 
688 
689  // ========================================================================
690  // simple adaptive integration
691  // ========================================================================
692  double NumericalDefiniteIntegral::QAG ( _Function* F ) const
693  {
694  if( 0 == F ) { Exception("QAG::invalid function") ; }
695 
696  // allocate workspace
697  if( 0 == ws () ) { allocate () ; }
698 
699  // integration limits
700  const double low = std::min ( m_a , m_b ) ;
701  const double high = std::max ( m_a , m_b ) ;
702 
703  int ierror =
704  gsl_integration_qag ( F->fn ,
705  low , high ,
706  m_epsabs , m_epsrel ,
707  size () , (int) rule() , ws ()->ws ,
708  &m_result , &m_error );
709 
710  if( ierror ) { gsl_error( "NumericalDefiniteIntegral::QAG " ,
711  __FILE__ , __LINE__ , ierror ) ; }
712 
713  // sign
714  if ( m_a > m_b ) { m_result *= -1 ; }
715 
716  return m_result ;
717  }
718  // ========================================================================
719 
720  // ========================================================================
721  // adaptive integration with singularities
722  // ========================================================================
723  double NumericalDefiniteIntegral::QAGS ( _Function* F ) const
724  {
725  if( 0 == F ) { Exception("QAG::invalid function") ; }
726 
727  if ( m_a == m_b )
728  {
729  m_result = 0 ;
730  m_error = 0 ; // EXACT !
731  return m_result ;
732  }
733 
734  // allocate workspace
735  if( 0 == ws () ) { allocate () ; }
736 
737  // integration limits
738  const double low = std::min ( m_a , m_b ) ;
739  const double high = std::max ( m_a , m_b ) ;
740 
741  int ierror =
742  gsl_integration_qags ( F->fn ,
743  low , high ,
744  m_epsabs , m_epsrel ,
745  size () , ws()->ws ,
746  &m_result , &m_error );
747 
748  if( ierror ) { gsl_error( "NumericalDefiniteIntegral::QAGS " ,
749  __FILE__ , __LINE__ , ierror ) ; }
750 
751  // sign
752  if ( m_a > m_b ) { m_result *= -1 ; }
753 
754  return m_result ;
755  }
756  // ========================================================================
757 
758 
759  } // end of namespace GaudiMathImplementation
760 } // end of namespace Genfun
761 
762 
763 
764 // ============================================================================
765 // The END
766 // ============================================================================
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::Exception
StatusCode Exception(const std::string &message, const StatusCode &sc=StatusCode::FAILURE) const
Definition: NumericalDefiniteIntegral.cpp:468
GaudiException
Define general base for Gaudi exception.
Definition: GaudiException.h:24
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::_Workspace
Definition: NumericalDefiniteIntegral.cpp:47
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::_Workspace::ws
gsl_integration_workspace * ws
Definition: NumericalDefiniteIntegral.cpp:48
GaudiMath::Derivative
Genfun::GaudiMathImplementation::NumericalDerivative Derivative
Definition: GaudiMath.h:31
merge_files.sc
int sc
Definition: merge_files.py:153
NumericalDerivative.h
gaudirun.type
string type
Definition: gaudirun.py:126
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::_Function::fn
gsl_function * fn
Definition: NumericalDefiniteIntegral.cpp:50
NumericalDefiniteIntegral.h
StatusCode
This class is used for returning status codes from appropriate routines.
Definition: StatusCode.h:30
GaudiException.h
Genfun::GaudiMathImplementation::GSL_Adaptor
double GSL_Adaptor(double x, void *params)
Definition: Helpers.cpp:38
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::_Function
Definition: NumericalDefiniteIntegral.cpp:49
Genfun::GaudiMathImplementation::NumericalDerivative
Numerical derivative (using GSL adaptive numerical differentiation)
Definition: NumericalDerivative.h:36
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::allocate
_Workspace * allocate() const
allocate the integration workspace
Definition: NumericalDefiniteIntegral.cpp:574
min
#define min(a, b)
Definition: CallbackStreamBuf.cpp:12
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::Points
std::vector< double > Points
typedef for vector of singular points
Definition: NumericalDefiniteIntegral.h:76
Helpers.h
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::partial
virtual Genfun::Derivative partial(unsigned int index) const
Derivatives.
Definition: NumericalDefiniteIntegral.cpp:499
GaudiMath::Integration::KronrodRule
KronrodRule
integration rule
Definition: Integration.h:36
Genfun::GaudiMathImplementation::GSL_Helper
the simple structure to be used for adaption interface Genfun::AbsFunction to gsl_function structure ...
Definition: Helpers.h:28
ana.i
list i
Definition: ana.py:128
Tuples::Type
Type
the list of available types for ntuples
Definition: TupleObj.h:63
Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::rule
GaudiMath::Integration::KronrodRule rule() const
integration rule
Definition: NumericalDefiniteIntegral.h:347
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:72