NumericalDefiniteIntegral.h

Go to the documentation of this file.

#ifndef GAUDIMATH_NUMERICALDEFINITEINTEGRAL_H
#define GAUDIMATH_NUMERICALDEFINITEINTEGRAL_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 NumericalDefiniteIntegral : public AbsFunction
    {
    public:
      struct _Workspace { gsl_integration_workspace* ws ; };
      struct _Function  ;
    public:

      typedef std::vector<double>   Points ;

    public:

      FUNCTION_OBJECT_DEF( NumericalDefiniteIntegral )

    public:

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


      NumericalDefiniteIntegral
      ( const AbsFunction&                  function           ,
        const size_t                        index              ,
        const double                        a                  ,
        const double                        b                  ,
        const Points&                       points             ,
        const double                        epsabs    = 1e-9   ,
        const double                        epsrel    = 1.e-6  ,
        const size_t                        size      = 1000   ) ;


      NumericalDefiniteIntegral
      ( const AbsFunction&                      function           ,
        const size_t                            index              ,
        const double                            a                  ,
        const GaudiMath::Integration::Inf       b         =
        GaudiMath::Integration::Infinity                           ,
        const double                            epsabs    = 1e-9   ,
        const double                            epsrel    = 1.e-6  ,
        const size_t                            size      = 1000   ) ;

      NumericalDefiniteIntegral
      ( const AbsFunction&                      function           ,
        const size_t                            index              ,
        const GaudiMath::Integration::Inf       a                  ,
        const double                            b                  ,
        const double                            epsabs    = 1e-9   ,
        const double                            epsrel    = 1.e-6  ,
        const size_t                            size      = 1000   ) ;

      NumericalDefiniteIntegral
      ( const AbsFunction&                      function           ,
        const size_t                            index              ,
        // FIXME: The next two arguments should be "double" but are "float" to resolve the
        // ambiguity with the first constructor when providing
        // '(const Genfun::AbsFunction, const unsigned int, const double, const double)'
        const float                             epsabs    = 1e-9   ,
        const float                             epsrel    = 1.e-6  ,
        const size_t                            size      = 1000   ) ;

      NumericalDefiniteIntegral ( const NumericalDefiniteIntegral& ) ;

      virtual ~NumericalDefiniteIntegral() = 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 ; }
      double             b        () const { return         m_b ; }
      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::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:

      // defautl constructor is disabled
      NumericalDefiniteIntegral();
      // assignement is disabled
      NumericalDefiniteIntegral& operator=( const NumericalDefiniteIntegral& );

    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        ;
      double                              m_b        ;
      bool                                m_ia       ;
      bool                                m_ib       ;

      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 ;
      mutable  Argument                   m_argF     ;

    };

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

// ============================================================================
// The END
// ============================================================================
#endif // GAUDIMATH_NUMERICALDEFINITEINTEGRAL_H
// ============================================================================