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"

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

      ~NumericalIndefiniteIntegral() override = default;

    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; }
      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;

      double m_epsabs;
      double m_epsrel;

      mutable double m_result = GSL_NEGINF;
      mutable double m_error  = GSL_POSINF;

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

      mutable Argument m_argument;
    };
    FUNCTION_OBJECT_IMP( NumericalIndefiniteIntegral )

  } // 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_NUMERICALINDEFINITEINTEGRAL_H
// ============================================================================