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