The Gaudi Framework  v28r3 (cc1cf868)
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_epsabs ( epsabs )
94  , m_epsrel ( epsrel )
95  //
96  , m_size ( size )
97  , m_argument ( function.dimensionality() )
98  {
99  if ( GaudiMath::Integration::Fixed == m_rule )
100  { m_rule = GaudiMath::Integration::Default ; }
101  if ( m_index >= m_DIM )
102  { Exception("::constructor: invalid variable index") ; }
103  }
104  // ========================================================================
105 
116  NumericalIndefiniteIntegral::NumericalIndefiniteIntegral
117  ( const AbsFunction& function ,
118  const size_t index ,
119  const double a ,
120  const Points& points ,
121  const GaudiMath::Integration::Limit limit ,
122  const double epsabs ,
123  const double epsrel ,
124  const size_t size )
125  : AbsFunction ()
126  , m_function ( function.clone() )
127  , m_DIM ( function.dimensionality() )
128  , m_index ( index )
129  , m_a ( a )
130  , m_limit ( limit )
131  , m_type ( GaudiMath::Integration:: Other )
132  , m_category ( GaudiMath::Integration:: Singular )
133  , m_rule ( GaudiMath::Integration:: Fixed )
134  , m_points ( points )
135  , m_epsabs ( epsabs )
136  , m_epsrel ( epsrel )
137  //
138  , m_size ( size )
139  , m_argument ( function.dimensionality() )
140  {
141  if ( m_index >= m_DIM )
142  { Exception("::constructor: invalid variable index") ; }
143  m_points.push_back( a ) ;
144  std::sort( m_points.begin() , m_points.end() ) ;
145  m_points.erase ( std::unique( m_points.begin () ,
146  m_points.end () ) , m_points.end() );
147  }
148 
149  // ========================================================================
157  // ========================================================================
158  NumericalIndefiniteIntegral::NumericalIndefiniteIntegral
159  ( const AbsFunction& function ,
160  const size_t index ,
161  const GaudiMath::Integration::Limit limit ,
162  const double epsabs ,
163  const double epsrel ,
164  const size_t size )
165  : AbsFunction ()
166  , m_function ( function.clone() )
167  , m_DIM ( function.dimensionality() )
168  , m_index ( index )
169  , m_a ( GSL_NEGINF ) // should not be used!
170  , m_limit ( limit )
171  , m_type ( GaudiMath::Integration:: Other )
172  , m_category ( GaudiMath::Integration:: Infinite )
173  , m_rule ( GaudiMath::Integration:: Fixed )
174  , m_epsabs ( epsabs )
175  , m_epsrel ( epsrel )
176  , m_size ( size )
177  , m_argument ( function.dimensionality() )
178  {
179  if ( m_index >= m_DIM )
180  { Exception("::constructor: invalid variable index") ; }
181  }
182  // ========================================================================
183 
184 
185  // ========================================================================
187  // ========================================================================
188  NumericalIndefiniteIntegral::
189  NumericalIndefiniteIntegral( const NumericalIndefiniteIntegral& right )
190  : AbsFunction()
191  , m_function ( right.m_function->clone() )
192  , m_DIM ( right.m_DIM )
193  , m_index ( right.m_index )
194  , m_a ( right.m_a )
195  , m_limit ( right.m_limit )
196  , m_type ( right.m_type )
197  , m_category ( right.m_category )
198  , m_rule ( right.m_rule )
199  , m_points ( right.m_points )
200  , m_epsabs ( right.m_epsabs )
201  , m_epsrel ( right.m_epsrel )
202  , m_size ( right.m_size )
203  , m_argument ( right.m_argument )
204  {
205  }
206  // ========================================================================
207 
208 
209  // ========================================================================
210  // throw the exception
211  // ========================================================================
212  StatusCode NumericalIndefiniteIntegral::Exception
213  ( const std::string& message ,
214  const StatusCode& sc ) const
215  {
216  throw GaudiException( "NumericalIndefiniteIntegral::" + message ,
217  "*GaudiMath*" , sc ) ;
218  return sc ;
219  }
220  // ========================================================================
221 
222  // ========================================================================
224  // ========================================================================
225  double NumericalIndefiniteIntegral::operator()
226  ( const double argument ) const
227  {
228  // reset the result and the error
229  m_result = GSL_NEGINF ;
230  m_error = GSL_POSINF ;
231 
232  // check the argument
233  if( 1 != m_DIM ) { Exception ( "operator(): invalid argument size " ) ; };
234 
235  m_argument[0] = argument ;
236  return (*this) ( m_argument );
237  }
238  // ========================================================================
239 
240  // ========================================================================
242  // ========================================================================
243  Genfun::Derivative
244  NumericalIndefiniteIntegral::partial ( unsigned int idx ) const
245  {
246  if ( idx >= m_DIM ) {
247  Exception ( "::partial(i): invalid variable index " ) ;
248  };
249  if ( idx != m_index ) {
250  const AbsFunction& aux = NumericalDerivative( *this , idx ) ;
251  return Genfun::FunctionNoop( &aux ) ;
252  } else if ( GaudiMath::Integration::VariableLowLimit == limit () ) {
253  const AbsFunction& aux = -1 * function() ;
254  return Genfun::FunctionNoop( &aux ) ;
255  }
256  const AbsFunction& aux = function() ;
257  return Genfun::FunctionNoop( &aux ) ;
258  }
259 
260  // ========================================================================
262  // ========================================================================
263  double NumericalIndefiniteIntegral::operator()( const Argument& argument ) const
264  {
265  // reset the result and the error
266  m_result = GSL_NEGINF ;
267  m_error = GSL_POSINF ;
268 
269  // check the argument
270  if( argument.dimension() != m_DIM ) {
271  Exception ( "operator(): invalid argument size " ) ;
272  };
273 
274  // copy the argument
275  for( size_t i = 0 ; i < m_DIM ; ++i ){ m_argument[i] = argument[i];}
276 
277  // create the helper object
278  GSL_Helper helper( *m_function , m_argument , m_index );
279 
280  // use GSL to evaluate the numerical derivative
281  gsl_function F ;
282  F.function = &GSL_Adaptor ;
283  F.params = &helper ;
284  _Function F1 ;
285  F1.fn = &F ;
286 
287  switch (category() ) {
288  case GaudiMath::Integration::Infinite : return QAGI( &F1 ) ;
289  case GaudiMath::Integration::Singular : return QAGP( &F1 ) ;
290  case GaudiMath::Integration::Finite :
291  switch (type()) {
292  case GaudiMath::Integration::NonAdaptive : return QNG( &F1 ) ;
293  case GaudiMath::Integration::Adaptive : return QAG( &F1 ) ;
294  case GaudiMath::Integration::AdaptiveSingular : return QAGS( &F1 ) ;
295  default: Exception ( "::operator(): invalid type " );
296  }
297  break;
298  default : Exception ( "::operator(): invalid category " );
299  }
300  return 0 ;
301  }
302  // ========================================================================
303 
304  // ========================================================================
306  // ========================================================================
307  NumericalIndefiniteIntegral::_Workspace*
308  NumericalIndefiniteIntegral::allocate() const
309  {
310  if ( !m_ws ) {
311  m_ws.reset( new _Workspace() );
312  m_ws->ws = gsl_integration_workspace_alloc( size () );
313  if ( !m_ws->ws ) { Exception ( "allocate()::invalid workspace" ) ; };
314  }
315  return m_ws.get() ;
316  }
317  // ========================================================================
318 
319  // ========================================================================
320  // adaptive integration on infinite interval
321  // ========================================================================
322  double NumericalIndefiniteIntegral::QAGI ( _Function* F ) const
323  {
324  // check the argument
325  if( !F ) { Exception("::QAGI: invalid function"); }
326 
327  const double x = m_argument[m_index] ;
328 
329  // allocate workspace
330  allocate() ;
331 
332  int ierror = 0 ;
333  switch ( limit() ) {
334  case GaudiMath::Integration::VariableLowLimit :
335  ierror = gsl_integration_qagiu ( F->fn , x ,
336  m_epsabs , m_epsrel ,
337  size () , ws()->ws ,
338  &m_result , &m_error ) ; break ;
339  case GaudiMath::Integration::VariableHighLimit :
340  ierror = gsl_integration_qagil ( F->fn , x ,
341  m_epsabs , m_epsrel ,
342  size () , ws()->ws ,
343  &m_result , &m_error ) ; break ;
344  default :
345  Exception ( "::QAGI: invalid mode" ) ;
346  };
347 
348  if( ierror ) { gsl_error( "NumericalIndefiniteIntegral::QAGI" ,
349  __FILE__ , __LINE__ , ierror ) ;}
350 
351  return m_result ;
352  }
353  // ========================================================================
354 
355  // ========================================================================
356  // adaptive integration with known singular points
357  // ========================================================================
358  double NumericalIndefiniteIntegral::QAGP( _Function* F ) const
359  {
360  if( !F ) { Exception("QAGP::invalid function") ; }
361 
362  const double x = m_argument[m_index] ;
363  if ( m_a == x ) {
364  m_result = 0 ;
365  m_error = 0 ; // EXACT !
366  return m_result ;
367  }
368 
369  // no known singular points ?
370  if( points().empty() ) { return QAGS( F ) ; }
371 
372  // integration limits
373  const auto limits = std::minmax( m_a, x );
374 
375  // "active" singular points
376  auto lower = std::lower_bound( points().begin(), points().end(), limits.first ) ;
377  auto upper = std::upper_bound( points().begin(), points().end(), limits.second ) ;
378 
379  Points pnts; pnts.reserve( std::distance(lower,upper) );
380  std::copy( lower , upper , std::back_inserter(pnts) );
381  if ( *lower != limits.first ) { pnts.insert( pnts.begin () , limits.first ) ; }
382  if ( *upper != limits.second ) { pnts.insert( pnts.end () , limits.second ) ; }
383 
384  // use GSL
385  int ierror = gsl_integration_qagp( F->fn ,
386  pnts.data(), pnts.size(),
387  m_epsabs , m_epsrel ,
388  size () , ws()->ws ,
389  &m_result , &m_error ) ;
390 
391  if( ierror ) { gsl_error( "NumericalIndefiniteIntegral::QAGI " ,
392  __FILE__ , __LINE__ , ierror ) ; }
393 
394  // sign
395  if ( GaudiMath::Integration::VariableHighLimit == limit()
396  && x < m_a ) { m_result *= -1 ; }
397  else if ( GaudiMath::Integration::VariableLowLimit == limit()
398  && x > m_a ) { m_result *= -1 ; }
399 
400  return m_result ;
401  }
402  // ========================================================================
403 
404  // ========================================================================
405  // non-adaptive integration
406  // ========================================================================
407  double NumericalIndefiniteIntegral::QNG ( _Function* F ) const
408  {
409  if( !F ) { Exception("QNG::invalid function") ; }
410 
411  const double x = m_argument[m_index] ;
412  if ( m_a == x ) {
413  m_result = 0 ;
414  m_error = 0 ; // EXACT !
415  return m_result ;
416  }
417 
418  // integration limits
419  const double a = std::min ( m_a , x ) ;
420  const double b = std::max ( m_a , x ) ;
421 
422  size_t neval = 0 ;
423  int ierror =
424  gsl_integration_qng ( F->fn ,
425  a , b ,
426  m_epsabs , m_epsrel ,
427  &m_result , &m_error , &neval ) ;
428 
429  if( ierror ) { gsl_error( "NumericalIndefiniteIntegral::QNG " ,
430  __FILE__ , __LINE__ , ierror ) ; }
431 
432  // sign
433  if ( GaudiMath::Integration::VariableHighLimit == limit()
434  && x < m_a ) { m_result *= -1 ; }
435  else if ( GaudiMath::Integration::VariableLowLimit == limit()
436  && x > m_a ) { m_result *= -1 ; }
437 
438  return m_result ;
439  }
440 
441  // ========================================================================
442  // simple adaptive integration
443  // ========================================================================
444  double NumericalIndefiniteIntegral::QAG ( _Function* F ) const
445  {
446  if( !F ) { Exception("QAG::invalid function") ; }
447 
448  const double x = m_argument[m_index] ;
449  if ( m_a == x )
450  {
451  m_result = 0 ;
452  m_error = 0 ; // EXACT !
453  return m_result ;
454  }
455 
456  // allocate workspace
457  allocate () ;
458 
459  // integration limits
460  const auto limits = std::minmax( m_a, x );
461 
462  int ierror =
463  gsl_integration_qag ( F->fn ,
464  limits.first, limits.second,
465  m_epsabs , m_epsrel ,
466  size () , (int) rule() , ws ()->ws ,
467  &m_result , &m_error );
468 
469  if( ierror ) { gsl_error( "NumericalIndefiniteIntegral::QAG " ,
470  __FILE__ , __LINE__ , ierror ) ; }
471 
472  // sign
473  if ( GaudiMath::Integration::VariableHighLimit == limit()
474  && x < m_a ) { m_result *= -1 ; }
475  else if ( GaudiMath::Integration::VariableLowLimit == limit()
476  && x > m_a ) { m_result *= -1 ; }
477 
478  return m_result ;
479  }
480  // ========================================================================
481 
482  // ========================================================================
483  // adaptive integration with singularities
484  // ========================================================================
485  double NumericalIndefiniteIntegral::QAGS ( _Function* F ) const
486  {
487  if( !F ) { Exception("QAG::invalid function") ; }
488 
489  const double x = m_argument[m_index] ;
490  if ( m_a == x ) {
491  m_result = 0 ;
492  m_error = 0 ; // EXACT !
493  return m_result ;
494  }
495 
496  // allocate workspace
497  allocate () ;
498 
499  // integration limits
500  const auto limits = std::minmax( m_a, x );
501 
502  int ierror =
503  gsl_integration_qags ( F->fn ,
504  limits.first, limits.second,
505  m_epsabs , m_epsrel ,
506  size () , ws()->ws ,
507  &m_result , &m_error );
508 
509  if( ierror ) { gsl_error( "NumericalIndefiniteIntegral::QAGS " ,
510  __FILE__ , __LINE__ , ierror ) ; }
511 
512  // sign
513  if ( GaudiMath::Integration::VariableHighLimit == limit()
514  && x < m_a ) { m_result *= -1 ; }
515  else if ( GaudiMath::Integration::VariableLowLimit == limit()
516  && x > m_a ) { m_result *= -1 ; }
517 
518  return m_result ;
519  }
520  // ========================================================================
521  } // end of namespace GaudiMathImplementation
522 } // end of namespace Genfun
523 
524 // ============================================================================
525 // The END
526 // ============================================================================
Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::partial
Genfun::Derivative partial(unsigned int index) const override
Derivatives.
Definition: NumericalIndefiniteIntegral.cpp:244
std::copy
T copy(T...args)
GaudiException
Define general base for Gaudi exception.
Definition: GaudiException.h:23
std::distance
T distance(T...args)
std::minmax
T minmax(T...args)
NumericalIndefiniteIntegral.h
Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::Exception
StatusCode Exception(const std::string &message, const StatusCode &sc=StatusCode::FAILURE) const
Definition: NumericalIndefiniteIntegral.cpp:213
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)
std::vector::data
T data(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:308
Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::operator()
double operator()(double argument) const override
Function value.
Definition: NumericalIndefiniteIntegral.cpp:226
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:33
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::vector::insert
T insert(T...args)
std::vector::size
T size(T...args)
std::vector::begin
T begin(T...args)
Helpers.h
std::back_inserter
T back_inserter(T...args)
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
std::vector::reserve
T reserve(T...args)