NumericalDefiniteIntegral.h

Go to the documentation of this file.

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

#if defined( __clang__ ) || defined( __CLING__ )
#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Winconsistent-missing-override"
#elif defined( __GNUC__ ) && __GNUC__ >= 5
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wsuggest-override"
#endif

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& );

    public:
      unsigned int dimensionality() const override { return m_DIM; }

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

      bool hasAnalyticDerivative() const override { return true; }

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

    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 = 0;
      size_t                             m_index;

      double m_a = -std::numeric_limits<double>::infinity();
      double m_b = std::numeric_limits<double>::infinity();

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

      Points m_points;

      double m_epsabs;
      double m_epsrel;

      mutable double m_result = -std::numeric_limits<double>::infinity();
      mutable double m_error  = std::numeric_limits<double>::infinity();

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

      mutable Argument m_argument;
      mutable Argument m_argF;
    };
    FUNCTION_OBJECT_IMP( NumericalDefiniteIntegral )

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

#if defined( __clang__ ) || defined( __CLING__ )
#pragma clang diagnostic pop
#elif defined( __GNUC__ ) && __GNUC__ >= 5
#pragma GCC diagnostic pop
#endif

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