All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
NumericalDefiniteIntegral.h
Go to the documentation of this file.
1 #ifndef GAUDIMATH_NUMERICALDEFINITEINTEGRAL_H
2 #define GAUDIMATH_NUMERICALDEFINITEINTEGRAL_H 1
3 // ============================================================================
4 // Include files
5 // ============================================================================
6 // STD & STL
7 // ============================================================================
8 #include <string>
9 // ============================================================================
10 // GaudiKernel
11 // ============================================================================
12 #include "GaudiKernel/StatusCode.h"
13 // ============================================================================
14 // GaudiMath
15 // ============================================================================
16 #include "GaudiMath/Integration.h"
17 // ============================================================================
18 // CLHEP
19 // ============================================================================
20 #include "CLHEP/GenericFunctions/AbsFunction.hh"
21 // ============================================================================
22 #include "gsl/gsl_integration.h"
23 
24 #if defined(__clang__) || defined(__CLING__)
25 #pragma clang diagnostic push
26 #pragma clang diagnostic ignored "-Winconsistent-missing-override"
27 #elif defined(__GNUC__) && __GNUC__ >= 5
28 #pragma GCC diagnostic push
29 #pragma GCC diagnostic ignored "-Wsuggest-override"
30 #endif
31 
32 namespace Genfun
33 {
34  namespace GaudiMathImplementation
35  {
79  class GAUDI_API NumericalDefiniteIntegral : public AbsFunction
80  {
81  public:
82  struct _Workspace { gsl_integration_workspace* ws ; };
83  struct _Function ;
84  public:
85 
87  typedef std::vector<double> Points ;
88 
89  public:
90 
92  FUNCTION_OBJECT_DEF( NumericalDefiniteIntegral )
93 
94  public:
95 
148  ( const AbsFunction& function ,
149  const size_t index ,
150  const double a ,
151  const double b ,
152  const GaudiMath::Integration::Type type =
154  const GaudiMath::Integration::KronrodRule rule =
156  const double epsabs = 1.e-10 ,
157  const double epsrel = 1.e-7 ,
158  const size_t size = 1000 );
159 
160 
190  ( const AbsFunction& function ,
191  const size_t index ,
192  const double a ,
193  const double b ,
194  const Points& points ,
195  const double epsabs = 1e-9 ,
196  const double epsrel = 1.e-6 ,
197  const size_t size = 1000 ) ;
198 
199 
226  ( const AbsFunction& function ,
227  const size_t index ,
228  const double a ,
229  const GaudiMath::Integration::Inf b =
231  const double epsabs = 1e-9 ,
232  const double epsrel = 1.e-6 ,
233  const size_t size = 1000 ) ;
234 
261  ( const AbsFunction& function ,
262  const size_t index ,
263  const GaudiMath::Integration::Inf a ,
264  const double b ,
265  const double epsabs = 1e-9 ,
266  const double epsrel = 1.e-6 ,
267  const size_t size = 1000 ) ;
268 
293  ( const AbsFunction& function ,
294  const size_t index ,
295  // FIXME: The next two arguments should be "double" but are "float" to resolve the
296  // ambiguity with the first constructor when providing
297  // '(const Genfun::AbsFunction, const unsigned int, const double, const double)'
298  const float epsabs = 1e-9 ,
299  const float epsrel = 1.e-6 ,
300  const size_t size = 1000 ) ;
301 
304 
306  ~NumericalDefiniteIntegral() override = default;
307 
308  public:
309 
311  unsigned int dimensionality() const override { return m_DIM ; }
312 
314  double operator() ( double argument ) const override ;
316  double operator() ( const Argument& argument ) const override ;
317 
319  bool hasAnalyticDerivative() const override { return true ;}
320 
322  Genfun::Derivative partial ( unsigned int index ) const override;
323 
324  public:
325 
327  const AbsFunction& function () const { return *m_function ; }
329  double a () const { return m_a ; }
330  double b () const { return m_b ; }
332  const Points& points () const { return m_points ; }
334  double epsabs () const { return m_epsabs ; }
336  double epsrel () const { return m_epsrel ; }
337 
339  double result () const { return m_result ; }
341  double error () const { return m_error ; }
342 
343  // maximal number of bisection integvals for adaptive algorithms
344  size_t size () const { return m_size ; }
345 
348  type () const { return m_type ; }
351  category () const { return m_category ; }
354  rule () const { return m_rule ; }
355 
356  protected:
357 
358  // adaptive integration on infinite intervals
359  double QAGI ( _Function* fun ) const ;
360  // adaptive integration with known singular points
361  double QAGP ( _Function* fun ) const ;
362  // non-adaptive integration
363  double QNG ( _Function* fun ) const ;
364  // adaptive integration
365  double QAG ( _Function* fun ) const ;
366  // adaptive integral with singularities
367  double QAGS ( _Function* fun ) const ;
368 
369  // allocate the integration workspace
370  _Workspace* allocate () const ;
371  // the integration workspace
372  _Workspace* ws () const
373  { return m_ws.get() ; };
374 
375  // throw the exception
376  StatusCode Exception
377  ( const std::string& message ,
378  const StatusCode& sc = StatusCode::FAILURE ) const ;
379 
380  private:
381 
382  // defautl constructor is disabled
384  // assignement is disabled
386 
387  public:
388  struct gsl_ws_deleter {
389  void operator()(_Workspace *p) const {
390  if (p) gsl_integration_workspace_free ( p->ws ) ;
391  delete p;
392  }
393  };
394  private:
395 
397  size_t m_DIM ;
398  size_t m_index ;
399 
400  double m_a ;
401  double m_b ;
402  bool m_ia ;
403  bool m_ib ;
404 
408 
409  Points m_points ;
411 
412  double m_epsabs ;
413  double m_epsrel ;
414 
415  mutable double m_result ;
416  mutable double m_error ;
417 
418  size_t m_size ;
420 
421  mutable Argument m_argument ;
422  mutable Argument m_argF ;
423 
424  };
426  FUNCTION_OBJECT_IMP( NumericalDefiniteIntegral )
427 
428  } // end of namespace GaudiMathImplementation
429 } // end of namespace Genfun
430 
431 #if defined(__clang__) || defined(__CLING__)
432 #pragma clang diagnostic pop
433 #elif defined(__GNUC__) && __GNUC__ >= 5
434 #pragma GCC diagnostic pop
435 #endif
436 
437 // ============================================================================
438 // The END
439 // ============================================================================
440 #endif // GAUDIMATH_NUMERICALDEFINITEINTEGRAL_H
441 // ============================================================================
GaudiMath::Integration::Category category() const
integration category
Category
integration category
Definition: Integration.h:30
Genfun::GaudiMathImplementation::NumericalDerivative Derivative
Definition: GaudiMath.h:29
PropertyMgr & operator=(const PropertyMgr &)=delete
STL class.
bool hasAnalyticDerivative() const override
Does this function have an analytic derivative?
This class is used for returning status codes from appropriate routines.
Definition: StatusCode.h:26
Type
type of integration (for finite limits)
Definition: Integration.h:25
GaudiMath::Integration::Type type() const
integration type
std::vector< double > Points
typedef for vector of singular points
virtual Out operator()(const vector_of_const_< In > &inputs) const =0
double fun(const std::vector< double > &x)
Definition: PFuncTest.cpp:26
GaudiMath.h GaudiMath/GaudiMath.h.
Definition: Adapters.h:13
KronrodRule
integration rule
Definition: Integration.h:34
CLHEP.
Definition: IEqSolver.h:13
#define GAUDI_API
Definition: Kernel.h:107
Type
the list of available types for ntuples
Definition: TupleObj.h:79
GaudiMath::Integration::KronrodRule rule() const
integration rule
collection of common types for classes NumericalIndefiniteIntegral and NumericalDefiniteIntegral ...
This class allows the numerical evaluation of the following functions: