NumericalIndefiniteIntegral.h

Go to the documentation of this file.

#ifndef GAUDIMATH_NUMERICALINDEFINITEINTEGRAL_H
#define GAUDIMATH_NUMERICALINDEFINITEINTEGRAL_H 1
// ============================================================================
// Include files
// ============================================================================
// STD & STL
// ============================================================================
#include <string>
// ============================================================================
// GaudiKernel
// ============================================================================
#include "GaudiKernel/StatusCode.h"
// ============================================================================
// GaudiMath
// ============================================================================
#include "GaudiMath/Integration.h"
// ============================================================================
// CLHEP
// ============================================================================
#include "CLHEP/GenericFunctions/AbsFunction.hh"
// ============================================================================
#include "gsl/gsl_integration.h"



namespace Genfun
{
  namespace GaudiMathImplementation
  {
    class GAUDI_API NumericalIndefiniteIntegral : public AbsFunction
    {
    public:
      struct _Workspace { gsl_integration_workspace* ws ; };
      struct _Function  ;
    public:

      typedef std::vector<double>   Points ;

    public:

      FUNCTION_OBJECT_DEF( NumericalIndefiniteIntegral )

    public:

      NumericalIndefiniteIntegral
      ( const AbsFunction&                         function                   ,
        const size_t                               index                      ,
        const double                               a                          ,
        const GaudiMath::Integration::Limit        limit  =
        GaudiMath::Integration::VariableHighLimit ,
        const GaudiMath::Integration::Type         type   =
        GaudiMath::Integration::Adaptive          ,
        const GaudiMath::Integration::KronrodRule  rule   =
        GaudiMath::Integration::Default           ,
        const double                               epsabs = 1e-10             ,
        const double                               epsrel = 1.e-7             ,
        const size_t                               size   = 1000              );

      NumericalIndefiniteIntegral
      ( const AbsFunction&                  function                      ,
        const size_t                        index                         ,
        const double                        a                             ,
        const Points&                       points                        ,
        const GaudiMath::Integration::Limit limit     =
        GaudiMath::Integration::VariableHighLimit ,
        const double                        epsabs    = 1e-9              ,
        const double                        epsrel    = 1.e-6             ,
        const size_t                        size      = 1000              ) ;

      NumericalIndefiniteIntegral
      ( const AbsFunction&                  function                      ,
        const size_t                        index                         ,
        const GaudiMath::Integration::Limit limit     =
        GaudiMath::Integration::VariableHighLimit ,
        const double                        epsabs    = 1e-9              ,
        const double                        epsrel    = 1.e-6             ,
        const size_t                        size      = 1000              ) ;

      NumericalIndefiniteIntegral
      ( const NumericalIndefiniteIntegral& ) ;

      virtual ~NumericalIndefiniteIntegral() = default;

    public:

      virtual unsigned int dimensionality() const { return m_DIM ; }

      virtual double operator() ( double          argument ) const ;
      virtual double operator() ( const Argument& argument ) const ;

      virtual bool hasAnalyticDerivative() const { return true ;}

      virtual Genfun::Derivative partial ( unsigned int index ) const;

    public:

      const AbsFunction& function () const { return *m_function ; }
      double             a        () const { return         m_a ; }
      const Points&      points   () const { return    m_points ; }
      double             epsabs   () const { return    m_epsabs ; }
      double             epsrel   () const { return    m_epsrel ; }

      double             result   () const { return    m_result ; }
      double             error    () const { return     m_error ; }

      // maximal number of bisection integvals for adaptive algorithms
      size_t             size     () const { return      m_size ; }

      GaudiMath::Integration::Limit
      limit    () const { return     m_limit ; }
      GaudiMath::Integration::Type
      type     () const { return      m_type ; }
      GaudiMath::Integration::Category
      category () const { return  m_category ; }
      GaudiMath::Integration::KronrodRule
      rule     () const { return      m_rule ; }

    protected:

      // adaptive integration on infinite intervals
      double                     QAGI ( _Function* fun ) const ;
      // adaptive integration  with known singular points
      double                     QAGP ( _Function* fun ) const ;
      // non-adaptive integration
      double                     QNG  ( _Function* fun ) const ;
      // adaptive integration
      double                     QAG  ( _Function* fun ) const ;
      // adaptive integral with singularities
      double                     QAGS ( _Function* fun ) const ;

      // allocate the integration workspace
      _Workspace*                allocate                () const ;
      // the integration workspace
      _Workspace*                ws                      () const
      { return m_ws.get(); }

      // throw the exception
      StatusCode Exception
      ( const std::string& message ,
        const StatusCode&  sc = StatusCode::FAILURE ) const ;

   private:

      // default constructor is disabled
      NumericalIndefiniteIntegral() ;
      // assignement operator is disabled
      NumericalIndefiniteIntegral& operator=
      ( const NumericalIndefiniteIntegral& ) ;

    public:

      struct gsl_ws_deleter {
          void operator()(_Workspace* p) const {
            if (p) gsl_integration_workspace_free ( p->ws ) ;
            delete p;
          }
      };
    private:

      std::unique_ptr<const AbsFunction>  m_function ;
      size_t                              m_DIM      ;
      size_t                              m_index    ;

      double                              m_a        ;

      GaudiMath::Integration::Limit       m_limit    ;
      GaudiMath::Integration::Type        m_type     ;
      GaudiMath::Integration::Category    m_category ;
      GaudiMath::Integration::KronrodRule m_rule     ;

      Points                              m_points   ;
      std::unique_ptr<double[]>           m_pdata    ;

      double                              m_epsabs   ;
      double                              m_epsrel   ;

      mutable double                      m_result   ;
      mutable double                      m_error    ;

      size_t                              m_size     ;
      mutable std::unique_ptr<_Workspace,gsl_ws_deleter> m_ws       ;

      mutable  Argument                   m_argument ;

    };

  } // end of namespace GaudiMathImplementation
} // end of namespace Genfun

// ============================================================================
// The END
// ============================================================================
#endif // GAUDIMATH_NUMERICALINDEFINITEINTEGRAL_H
// ============================================================================