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