•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
NumericalIndefiniteIntegral.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 // GaudiMath
15 // ============================================================================
16 #include "GaudiMath/NumericalIndefiniteIntegral.h"
17 #include "GaudiMath/NumericalDerivative.h"
18 // ============================================================================
19 // GaudiKernel
20 // ============================================================================
21 #include "GaudiKernel/GaudiException.h"
22 // ============================================================================
23 // local
24 // ============================================================================
25 #include "Helpers.h"
26 // ============================================================================
27 
28 // ============================================================================
33 // ============================================================================
34 
35 #ifdef __ICC
36 // disable icc remark #1572: floating-point equality and inequality comparisons are unreliable
37 // The comparison are meant
38 #pragma warning(disable:1572)
39 #endif
40 #ifdef WIN32
41 // Disable the warning
42 // C4996: 'std::copy': Function call with parameters that may be unsafe
43 // The parameters are checked
44 #pragma warning(disable:4996)
45 #endif
46 
47 namespace Genfun
48 {
49  namespace GaudiMathImplementation
50  {
51 
52  struct NumericalIndefiniteIntegral::_Function
53  { gsl_function* fn ; };
54 
55  // ========================================================================
56 
57  // ========================================================================
72  // ========================================================================
73  NumericalIndefiniteIntegral::NumericalIndefiniteIntegral
74  ( const AbsFunction& function ,
75  const size_t index ,
76  const double a ,
77  const GaudiMath::Integration::Limit limit ,
78  const GaudiMath::Integration::Type type ,
79  const GaudiMath::Integration::KronrodRule rule ,
80  const double epsabs ,
81  const double epsrel ,
82  const size_t size )
83  : AbsFunction ()
84  , m_function ( function.clone() )
85  , m_DIM ( function.dimensionality() )
86  , m_index ( index )
87  , m_a ( a )
88  , m_limit ( limit )
89  , m_type ( type )
90  , m_category ( GaudiMath::Integration::Finite )
91  , m_rule ( rule )
92  //
93  , m_points ( )
94  //
95  , m_epsabs ( epsabs )
96  , m_epsrel ( epsrel )
97  //
98  , m_result ( GSL_NEGINF )
99  , m_error ( GSL_POSINF )
100  //
101  , m_size ( size )
102  , m_argument ( function.dimensionality() )
103  {
104  if ( GaudiMath::Integration::Fixed == m_rule )
105  { m_rule = GaudiMath::Integration::Default ; }
106  if ( m_index >= m_DIM )
107  { Exception("::constructor: invalid variable index") ; }
108  }
109  // ========================================================================
110 
121  NumericalIndefiniteIntegral::NumericalIndefiniteIntegral
122  ( const AbsFunction& function ,
123  const size_t index ,
124  const double a ,
125  const Points& points ,
126  const GaudiMath::Integration::Limit limit ,
127  const double epsabs ,
128  const double epsrel ,
129  const size_t size )
130  : AbsFunction ()
131  , m_function ( function.clone() )
132  , m_DIM ( function.dimensionality() )
133  , m_index ( index )
134  , m_a ( a )
135  , m_limit ( limit )
136  , m_type ( GaudiMath::Integration:: Other )
137  , m_category ( GaudiMath::Integration:: Singular )
138  , m_rule ( GaudiMath::Integration:: Fixed )
139  , m_points ( points )
140  , m_epsabs ( epsabs )
141  , m_epsrel ( epsrel )
142  //
143  , m_result ( GSL_NEGINF )
144  , m_error ( GSL_POSINF )
145  //
146  , m_size ( size )
147  , m_argument ( function.dimensionality() )
148  {
149  if ( m_index >= m_DIM )
150  { Exception("::constructor: invalid variable index") ; }
151  m_pdata.reset( new double[ 2 + m_points.size() ] );
152  m_points.push_back( a ) ;
153  std::sort( m_points.begin() , m_points.end() ) ;
154  m_points.erase ( std::unique( m_points.begin () ,
155  m_points.end () ) , m_points.end() );
156  }
157 
158  // ========================================================================
166  // ========================================================================
167  NumericalIndefiniteIntegral::NumericalIndefiniteIntegral
168  ( const AbsFunction& function ,
169  const size_t index ,
170  const GaudiMath::Integration::Limit limit ,
171  const double epsabs ,
172  const double epsrel ,
173  const size_t size )
174  : AbsFunction ()
175  , m_function ( function.clone() )
176  , m_DIM ( function.dimensionality() )
177  , m_index ( index )
178  , m_a ( GSL_NEGINF ) // should not be used!
179  , m_limit ( limit )
180  , m_type ( GaudiMath::Integration:: Other )
181  , m_category ( GaudiMath::Integration:: Infinite )
182  , m_rule ( GaudiMath::Integration:: Fixed )
183  , m_points ( )
184  , m_epsabs ( epsabs )
185  , m_epsrel ( epsrel )
186  , m_result ( GSL_NEGINF )
187  , m_error ( GSL_POSINF )
188  , m_size ( size )
189  , m_argument ( function.dimensionality() )
190  {
191  if ( m_index >= m_DIM )
192  { Exception("::constructor: invalid variable index") ; }
193  }
194  // ========================================================================
195 
196 
197  // ========================================================================
199  // ========================================================================
200  NumericalIndefiniteIntegral::
201  NumericalIndefiniteIntegral( const NumericalIndefiniteIntegral& right )
202  : AbsFunction ()
203  , m_function ( right.m_function->clone() )
204  , m_DIM ( right.m_DIM )
205  , m_index ( right.m_index )
206  , m_a ( right.m_a )
207  , m_limit ( right.m_limit )
208  , m_type ( right.m_type )
209  , m_category ( right.m_category )
210  , m_rule ( right.m_rule )
211  , m_points ( right.m_points )
212  , m_epsabs ( right.m_epsabs )
213  , m_epsrel ( right.m_epsrel )
214  , m_result ( GSL_NEGINF )
215  , m_error ( GSL_POSINF )
216  , m_size ( right.m_size )
217  , m_argument ( right.m_argument )
218  {
219  m_pdata.reset( new double[ 2 + m_points.size() ] ) ;
220  }
221  // ========================================================================
222 
223 
224  // ========================================================================
225  // throw the exception
226  // ========================================================================
227  StatusCode NumericalIndefiniteIntegral::Exception
228  ( const std::string& message ,
229  const StatusCode& sc ) const
230  {
231  throw GaudiException( "NumericalIndefiniteIntegral::" + message ,
232  "*GaudiMath*" , sc ) ;
233  return sc ;
234  }
235  // ========================================================================
236 
237  // ========================================================================
239  // ========================================================================
240  double NumericalIndefiniteIntegral::operator()
241  ( const double argument ) const
242  {
243  // reset the result and the error
244  m_result = GSL_NEGINF ;
245  m_error = GSL_POSINF ;
246 
247  // check the argument
248  if( 1 != m_DIM ) { Exception ( "operator(): invalid argument size " ) ; };
249 
250  m_argument[0] = argument ;
251  return (*this) ( m_argument );
252  }
253  // ========================================================================
254 
255  // ========================================================================
257  // ========================================================================
258  Genfun::Derivative
259  NumericalIndefiniteIntegral::partial ( unsigned int idx ) const
260  {
261  if ( idx >= m_DIM )
262  { Exception ( "::partial(i): invalid variable index " ) ; };
263  if ( idx != m_index )
264  {
265  const AbsFunction& aux = NumericalDerivative( *this , idx ) ;
266  return Genfun::FunctionNoop( &aux ) ;
267  }
268  else if ( GaudiMath::Integration::VariableLowLimit == limit () )
269  {
270  const AbsFunction& aux = -1 * function() ;
271  return Genfun::FunctionNoop( &aux ) ;
272  }
273  const AbsFunction& aux = function() ;
274  return Genfun::FunctionNoop( &aux ) ;
275  }
276 
277  // ========================================================================
279  // ========================================================================
280  double NumericalIndefiniteIntegral::operator()
281  ( const Argument& argument ) const
282  {
283  // reset the result and the error
284  m_result = GSL_NEGINF ;
285  m_error = GSL_POSINF ;
286 
287  // check the argument
288  if( argument.dimension() != m_DIM )
289  { Exception ( "operator(): invalid argument size " ) ; };
290 
291  // copy the argument
292  {for( size_t i = 0 ; i < m_DIM ; ++i ){ m_argument[i] = argument[i];}}
293 
294  // create the helper object
295  GSL_Helper helper( *m_function , m_argument , m_index );
296 
297  // use GSL to evaluate the numerical derivative
298  gsl_function F ;
299  F.function = &GSL_Adaptor ;
300  F.params = &helper ;
301  _Function F1 ;
302  F1.fn = &F ;
303 
304  if ( GaudiMath::Integration::Infinite == category () )
305  { return QAGI ( &F1 ) ; } // RETURN
306  else if ( GaudiMath::Integration::Singular == category () )
307  { return QAGP ( &F1 ) ; } // RETURN
308  else if ( GaudiMath::Integration::Finite == category () )
309  if ( GaudiMath::Integration::NonAdaptive == type () )
310  { return QNG ( &F1 ) ; } // RETURN
311  else if ( GaudiMath::Integration::Adaptive == type () )
312  { return QAG ( &F1 ) ; } // RETURN
313  else if ( GaudiMath::Integration::AdaptiveSingular == type () )
314  { return QAGS ( &F1 ) ; } // RETURN
315  else
316  { Exception ( "::operator(): invalid type " ); }
317  else
318  { Exception ( "::operator(): invalid category " ); }
319 
320  return 0 ;
321  }
322  // ========================================================================
323 
324  // ========================================================================
326  // ========================================================================
327  NumericalIndefiniteIntegral::_Workspace*
328  NumericalIndefiniteIntegral::allocate() const
329  {
330  if ( !m_ws ) {
331  m_ws.reset( new _Workspace() );
332  m_ws->ws = gsl_integration_workspace_alloc( size () );
333  if ( !m_ws->ws ) { Exception ( "allocate()::invalid workspace" ) ; };
334  }
335  return m_ws.get() ;
336  }
337  // ========================================================================
338 
339  // ========================================================================
340  // adaptive integration on infinite interval
341  // ========================================================================
342  double NumericalIndefiniteIntegral::QAGI ( _Function* F ) const
343  {
344  // check the argument
345  if( !F ) { Exception("::QAGI: invalid function"); }
346 
347  const double x = m_argument[m_index] ;
348 
349  // allocate workspace
350  allocate() ;
351 
352  int ierror = 0 ;
353  switch ( limit() )
354  {
355  case GaudiMath::Integration::VariableLowLimit :
356  ierror = gsl_integration_qagiu ( F->fn , x ,
357  m_epsabs , m_epsrel ,
358  size () , ws()->ws ,
359  &m_result , &m_error ) ; break ;
360  case GaudiMath::Integration::VariableHighLimit :
361  ierror = gsl_integration_qagil ( F->fn , x ,
362  m_epsabs , m_epsrel ,
363  size () , ws()->ws ,
364  &m_result , &m_error ) ; break ;
365  default :
366  Exception ( "::QAGI: invalid mode" ) ;
367  };
368 
369  if( ierror ) { gsl_error( "NumericalIndefiniteIntegral::QAGI" ,
370  __FILE__ , __LINE__ , ierror ) ;}
371 
372  return m_result ;
373  }
374  // ========================================================================
375 
376  // ========================================================================
377  // adaptive integration with known singular points
378  // ========================================================================
379  double NumericalIndefiniteIntegral::QAGP( _Function* F ) const
380  {
381  if( !F ) { Exception("QAGP::invalid function") ; }
382 
383  const double x = m_argument[m_index] ;
384  if ( m_a == x )
385  {
386  m_result = 0 ;
387  m_error = 0 ; // EXACT !
388  return m_result ;
389  }
390 
391  // no known singular points ?
392  if( points().empty() ) { return QAGS( F ) ; }
393 
394  // integration limits
395  const double a = std::min ( m_a , x ) ;
396  const double b = std::max ( m_a , x ) ;
397 
398  // "active" singular points
399  auto lower = std::lower_bound( points().begin(), points().end(), a ) ;
400  auto upper = std::upper_bound( points().begin(), points().end(), b ) ;
401 
402  Points pnts ( upper - lower ) ;
403  std::copy( lower , upper , pnts.begin() );
404  if ( *lower != a ) { pnts.insert( pnts.begin () , a ) ; }
405  if ( *upper != b ) { pnts.insert( pnts.end () , b ) ; }
406  std::copy( pnts.begin() , pnts.end() , m_pdata.get() );
407  const size_t npts = pnts.size() ;
408 
409  // use GSL
410  int ierror =
411  gsl_integration_qagp ( F->fn ,
412  m_pdata.get() , npts ,
413  m_epsabs , m_epsrel ,
414  size () , ws()->ws ,
415  &m_result , &m_error ) ;
416 
417  if( ierror ) { gsl_error( "NumericalIndefiniteIntegral::QAGI " ,
418  __FILE__ , __LINE__ , ierror ) ; }
419 
420  // sign
421  if ( GaudiMath::Integration::VariableHighLimit == limit()
422  && x < m_a ) { m_result *= -1 ; }
423  else if ( GaudiMath::Integration::VariableLowLimit == limit()
424  && x > m_a ) { m_result *= -1 ; }
425 
426  return m_result ;
427  }
428  // ========================================================================
429 
430  // ========================================================================
431  // non-adaptive integration
432  // ========================================================================
433  double NumericalIndefiniteIntegral::QNG ( _Function* F ) const
434  {
435  if( !F ) { Exception("QNG::invalid function") ; }
436 
437  const double x = m_argument[m_index] ;
438  if ( m_a == x )
439  {
440  m_result = 0 ;
441  m_error = 0 ; // EXACT !
442  return m_result ;
443  }
444 
445  // integration limits
446  const double a = std::min ( m_a , x ) ;
447  const double b = std::max ( m_a , x ) ;
448 
449  size_t neval = 0 ;
450  int ierror =
451  gsl_integration_qng ( F->fn ,
452  a , b ,
453  m_epsabs , m_epsrel ,
454  &m_result , &m_error , &neval ) ;
455 
456  if( ierror ) { gsl_error( "NumericalIndefiniteIntegral::QNG " ,
457  __FILE__ , __LINE__ , ierror ) ; }
458 
459  // sign
460  if ( GaudiMath::Integration::VariableHighLimit == limit()
461  && x < m_a ) { m_result *= -1 ; }
462  else if ( GaudiMath::Integration::VariableLowLimit == limit()
463  && x > m_a ) { m_result *= -1 ; }
464 
465  return m_result ;
466  }
467 
468  // ========================================================================
469  // simple adaptive integration
470  // ========================================================================
471  double NumericalIndefiniteIntegral::QAG ( _Function* F ) const
472  {
473  if( !F ) { Exception("QAG::invalid function") ; }
474 
475  const double x = m_argument[m_index] ;
476  if ( m_a == x )
477  {
478  m_result = 0 ;
479  m_error = 0 ; // EXACT !
480  return m_result ;
481  }
482 
483  // allocate workspace
484  allocate () ;
485 
486  // integration limits
487  const double a = std::min ( m_a , x ) ;
488  const double b = std::max ( m_a , x ) ;
489 
490  int ierror =
491  gsl_integration_qag ( F->fn ,
492  a , b ,
493  m_epsabs , m_epsrel ,
494  size () , (int) rule() , ws ()->ws ,
495  &m_result , &m_error );
496 
497  if( ierror ) { gsl_error( "NumericalIndefiniteIntegral::QAG " ,
498  __FILE__ , __LINE__ , ierror ) ; }
499 
500  // sign
501  if ( GaudiMath::Integration::VariableHighLimit == limit()
502  && x < m_a ) { m_result *= -1 ; }
503  else if ( GaudiMath::Integration::VariableLowLimit == limit()
504  && x > m_a ) { m_result *= -1 ; }
505 
506  return m_result ;
507  }
508  // ========================================================================
509 
510  // ========================================================================
511  // adaptive integration with singularities
512  // ========================================================================
513  double NumericalIndefiniteIntegral::QAGS ( _Function* F ) const
514  {
515  if( !F ) { Exception("QAG::invalid function") ; }
516 
517  const double x = m_argument[m_index] ;
518  if ( m_a == x )
519  {
520  m_result = 0 ;
521  m_error = 0 ; // EXACT !
522  return m_result ;
523  }
524 
525  // allocate workspace
526  allocate () ;
527 
528  // integration limits
529  const double a = std::min ( m_a , x ) ;
530  const double b = std::max ( m_a , x ) ;
531 
532  int ierror =
533  gsl_integration_qags ( F->fn ,
534  a , b ,
535  m_epsabs , m_epsrel ,
536  size () , ws()->ws ,
537  &m_result , &m_error );
538 
539  if( ierror ) { gsl_error( "NumericalIndefiniteIntegral::QAGS " ,
540  __FILE__ , __LINE__ , ierror ) ; }
541 
542  // sign
543  if ( GaudiMath::Integration::VariableHighLimit == limit()
544  && x < m_a ) { m_result *= -1 ; }
545  else if ( GaudiMath::Integration::VariableLowLimit == limit()
546  && x > m_a ) { m_result *= -1 ; }
547 
548  return m_result ;
549  }
550  // ========================================================================
551 
552  } // end of namespace GaudiMathImplementation
553 } // end of namespace Genfun
554 
555 
556 // ============================================================================
557 // The END
558 // ============================================================================
Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::partial
Genfun::Derivative partial(unsigned int index) const override
Derivatives.
Definition: NumericalIndefiniteIntegral.cpp:259
std::copy
T copy(T...args)
GaudiException
Define general base for Gaudi exception.
Definition: GaudiException.h:23
NumericalIndefiniteIntegral.h
Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::Exception
StatusCode Exception(const std::string &message, const StatusCode &sc=StatusCode::FAILURE) const
Definition: NumericalIndefiniteIntegral.cpp:228
std::upper_bound
T upper_bound(T...args)
std::vector::end
T end(T...args)
Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::limit
GaudiMath::Integration::Limit limit() const
integration limit
Definition: NumericalIndefiniteIntegral.h:321
GaudiMath::Derivative
Genfun::GaudiMathImplementation::NumericalDerivative Derivative
Definition: GaudiMath.h:29
std::lower_bound
T lower_bound(T...args)
Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::_Workspace::ws
gsl_integration_workspace * ws
Definition: NumericalIndefiniteIntegral.h:86
std::unique
T unique(T...args)
Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::category
GaudiMath::Integration::Category category() const
integration category
Definition: NumericalIndefiniteIntegral.h:327
begin
auto begin(reverse_wrapper< T > &w)
Definition: reverse.h:48
Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::_Function::fn
gsl_function * fn
Definition: NumericalIndefiniteIntegral.cpp:53
std::string
STL class.
std::min
T min(T...args)
NumericalDerivative.h
std::vector::push_back
T push_back(T...args)
Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::rule
GaudiMath::Integration::KronrodRule rule() const
integration rule
Definition: NumericalIndefiniteIntegral.h:330
Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::_Function
Definition: NumericalIndefiniteIntegral.cpp:52
StatusCode
This class is used for returning status codes from appropriate routines.
Definition: StatusCode.h:26
Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::allocate
_Workspace * allocate() const
allocate the integration workspace
Definition: NumericalIndefiniteIntegral.cpp:328
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::NumericalIndefiniteIntegral::_Workspace
Definition: NumericalIndefiniteIntegral.h:86
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:37
std::unique_ptr::reset
T reset(T...args)
std::max
T max(T...args)
Genfun::GaudiMathImplementation::NumericalDerivative
Numerical derivative (using GSL adaptive numerical differentiation)
Definition: NumericalDerivative.h:42
Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::dimensionality
unsigned int dimensionality() const override
dimensionality of the problem
Definition: NumericalIndefiniteIntegral.h:285
Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::type
GaudiMath::Integration::Type type() const
integration type
Definition: NumericalIndefiniteIntegral.h:324
std::unique_ptr::get
T get(T...args)
std::vector::insert
T insert(T...args)
std::vector::size
T size(T...args)
std::vector::begin
T begin(T...args)
Helpers.h
GaudiMath::Integration::Limit
Limit
how to distinguish variable low and variable high limits
Definition: Integration.h:22
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