Splines.h

Go to the documentation of this file.

#ifndef GAUDIMATH_SPLINES_H
#define GAUDIMATH_SPLINES_H 1
// ============================================================================
// Include files
// ============================================================================
// STD/STL
// ============================================================================
#include <vector>
#include <utility>
#include <algorithm>
#include <memory>
// ============================================================================
// from CLHEP
// ============================================================================
#include "CLHEP/GenericFunctions/GenericFunctions.hh"
#include "CLHEP/GenericFunctions/Argument.hh"
#include "CLHEP/GenericFunctions/AbsFunction.hh"
// ============================================================================
// GaudiGSL/GaudiMath
// ============================================================================
#include "GaudiMath/Interpolation.h"
// ============================================================================
// GSL
// ============================================================================
#include "gsl/gsl_interp.h"
#include "gsl/gsl_spline.h"
// ============================================================================
#include "GaudiKernel/Kernel.h"

namespace Genfun
{
  namespace GaudiMathImplementation
  {

    class  GAUDI_API SplineBase
    {
    public:
      typedef std::vector<double>                    Data1D ;
      typedef std::vector<std::pair<double,double> > Data2D ;
    public:
      SplineBase
      ( const Data1D&                        x    ,
        const Data1D&                        y    ,
        const GaudiMath::Interpolation::Type type ) ;
      SplineBase
      ( const Data2D&                        data ,
        const GaudiMath::Interpolation::Type type ) ;
      template <class DATAX, class DATAY>
      SplineBase
      ( const GaudiMath::Interpolation::Type type    ,
        DATAX                                begin_x ,
        DATAX                                end_x   ,
        DATAY                                begin_y )
        : m_dim       ( end_x - begin_x )
        , m_x         ( new double[ end_x - begin_x ] )
        , m_y         ( new double[ end_x - begin_x ] )
        , m_type      ( type )
      {
        std::copy ( begin_x , end_x                         , m_x.get() ) ;
        std::copy ( begin_y , begin_y + ( end_x - begin_x ) , m_y.get() ) ;
      }
      template <class DATA>
      SplineBase
      ( const GaudiMath::Interpolation::Type type  ,
        DATA                                 begin ,
        DATA                                 end   )
        : m_dim       ( end - begin  )
        , m_x         ( new double[ end - begin ] )
        , m_y         ( new double[ end - begin ] )
        , m_type      ( type )
      {
        double* _x = m_x.get() ;
        double* _y = m_y.get() ;
        for ( auto it = begin ; end != it ; ++ it )
        {
          *_x++ = it -> first  ;
          *_y++ = it -> second ;
        };
      }
      SplineBase( const SplineBase& ) ;
      virtual ~SplineBase();
    public:
      double eval   ( const double x ) const ;
      double deriv  ( const double x ) const ;
      double deriv2 ( const double x ) const ;
      double integ  ( const double a , const double b ) const ;
    protected:
      // initialize
      void initialize () const ;
    private:
      // default constructor is disabled
      SplineBase() = delete;
      // assigenement operator  is disabled
      SplineBase& operator=( const SplineBase& ) = delete;
    private:
      mutable bool                   m_init   = false;
      size_t                         m_dim    ;
      std::unique_ptr<double[]>      m_x      ;
      std::unique_ptr<double[]>      m_y      ;
      mutable gsl_spline*            m_spline = nullptr;
      mutable gsl_interp_accel*      m_accel  = nullptr;
      GaudiMath::Interpolation::Type m_type   ;
    };

    class GAUDI_API GSLSpline : public AbsFunction
    {
    public:
      typedef SplineBase::Data1D Data1D ;
      typedef SplineBase::Data2D Data2D ;
    public:
      FUNCTION_OBJECT_DEF( GSLSpline )
    public:
      GSLSpline
      ( const Data1D&                        x    ,
        const Data1D&                        y    ,
        const GaudiMath::Interpolation::Type type ) ;
      GSLSpline
      ( const Data2D&                        data ,
        const GaudiMath::Interpolation::Type type ) ;
      template <class DATAX, class DATAY>
      GSLSpline
      ( const GaudiMath::Interpolation::Type type    ,
        DATAX                                begin_x ,
        DATAX                                end_x   ,
        DATAY                                begin_y )
        : AbsFunction ( )
        , m_spline( type , begin_x , end_x , begin_y )
      {}
      template <class DATA>
      GSLSpline
      ( const GaudiMath::Interpolation::Type type  ,
        DATA                                 begin ,
        DATA                                 end   )
        : AbsFunction ( )
        , m_spline( type , begin , end )
      {}
      GSLSpline ( const SplineBase& ) ;
      GSLSpline ( const GSLSpline& ) ;
      virtual ~GSLSpline() ;
    public:
      virtual double operator() ( double a          ) const ;
      virtual double operator() ( const Argument& x ) const ;
      virtual unsigned int dimensionality () const { return 1   ; }
      virtual bool  hasAnalyticDerivative() const { return true ; }
      virtual Genfun::Derivative partial( unsigned int i  ) const ;
    public:
      inline   const SplineBase& spline() const { return m_spline ; }
      operator const SplineBase&       () const { return spline() ; }
    private:
      GSLSpline() = delete;
      GSLSpline& operator=( const GSLSpline& ) = delete;
    private:
      // the actual spline function
      SplineBase m_spline ;
    };

    class GAUDI_API GSLSplineDeriv : public AbsFunction
    {
    public:
      typedef SplineBase::Data1D Data1D ;
      typedef SplineBase::Data2D Data2D ;
    public:
      FUNCTION_OBJECT_DEF( GSLSplineDeriv )
    public:
      GSLSplineDeriv
      ( const Data1D&                        x    ,
        const Data1D&                        y    ,
        const GaudiMath::Interpolation::Type type ) ;
      GSLSplineDeriv
      ( const Data2D&                        data ,
        const GaudiMath::Interpolation::Type type ) ;
      template <class DATAX, class DATAY>
      GSLSplineDeriv
      ( const GaudiMath::Interpolation::Type type    ,
        DATAX                                begin_x ,
        DATAX                                end_x   ,
        DATAY                                begin_y )
        : AbsFunction ( )
        , m_spline( type , begin_x , end_x , begin_y )
      {}
      template <class DATA>
      GSLSplineDeriv
      ( const GaudiMath::Interpolation::Type type  ,
        DATA                                 begin ,
        DATA                                 end   )
        : AbsFunction ( )
        , m_spline( type , begin , end )
      {}
      GSLSplineDeriv ( const SplineBase&     ) ;
      GSLSplineDeriv ( const GSLSplineDeriv& ) ;
      ~GSLSplineDeriv() override;
    public:
      double operator() ( double a          ) const override;
      double operator() ( const Argument& x ) const override;
      unsigned int dimensionality () const override { return 1   ; }
      bool  hasAnalyticDerivative() const override { return true ; }
      Genfun::Derivative partial( unsigned int i  ) const override ;
    public:
      inline   const SplineBase& spline() const { return m_spline ; }
      operator const SplineBase&       () const { return spline() ; }
    private:
      GSLSplineDeriv() ;
      GSLSplineDeriv& operator=( const GSLSplineDeriv& ) ;
    private:
      // the actual spline function
      SplineBase m_spline ;
    };

    class GAUDI_API GSLSplineDeriv2 : public AbsFunction
    {
    public:
      typedef SplineBase::Data1D Data1D ;
      typedef SplineBase::Data2D Data2D ;
    public:
      FUNCTION_OBJECT_DEF( GSLSplineDeriv2 )
    public:
      GSLSplineDeriv2
      ( const Data1D&                        x    ,
        const Data1D&                        y    ,
        const GaudiMath::Interpolation::Type type ) ;
      GSLSplineDeriv2
      ( const Data2D&                        data ,
        const GaudiMath::Interpolation::Type type ) ;
      template <class DATAX, class DATAY>
      GSLSplineDeriv2
      ( const GaudiMath::Interpolation::Type type    ,
        DATAX                                begin_x ,
        DATAX                                end_x   ,
        DATAY                                begin_y )
        : AbsFunction ( )
        , m_spline( type , begin_x , end_x , begin_y )
      {}
      template <class DATA>
      GSLSplineDeriv2
      ( const GaudiMath::Interpolation::Type type  ,
        DATA                                 begin ,
        DATA                                 end   )
        : AbsFunction ( )
        , m_spline( type , begin , end )
      {}
      GSLSplineDeriv2 ( const SplineBase&      ) ;
      GSLSplineDeriv2 ( const GSLSplineDeriv2& ) ;
      ~GSLSplineDeriv2() override;
    public:
      double operator() ( double a          ) const override;
      double operator() ( const Argument& x ) const override;
      unsigned int dimensionality () const override { return 1   ; }
      bool  hasAnalyticDerivative() const override { return true ; }
      Genfun::Derivative partial( unsigned int i  ) const override;
    public:
      inline   const SplineBase& spline() const { return m_spline ; }
      operator const SplineBase&       () const { return spline() ; }
    private:
      GSLSplineDeriv2() = delete;
      GSLSplineDeriv2& operator=( const GSLSplineDeriv2& ) = delete;
    private:
      // the actual spline function
      SplineBase m_spline ;
    };


    class GAUDI_API GSLSplineInteg : public AbsFunction
    {
    public:
      typedef SplineBase::Data1D Data1D ;
      typedef SplineBase::Data2D Data2D ;
    public:
      FUNCTION_OBJECT_DEF( GSLSplineInteg )
    public:
      GSLSplineInteg
      ( const Data1D&                        x        ,
        const Data1D&                        y        ,
        const GaudiMath::Interpolation::Type type     ,
        const double                         low  = 0 ) ;
      GSLSplineInteg
      ( const Data2D&                        data     ,
        const GaudiMath::Interpolation::Type type     ,
        const double                         low  = 0 ) ;
      template <class DATAX, class DATAY>
      GSLSplineInteg
      ( const GaudiMath::Interpolation::Type type    ,
        DATAX                                begin_x ,
        DATAX                                end_x   ,
        DATAY                                begin_y ,
        const double                         low     )
        : AbsFunction (      )
        , m_spline    ( type , begin_x , end_x , begin_y )
        , m_low       ( low  )
      {}
      template <class DATA>
      GSLSplineInteg
      ( const GaudiMath::Interpolation::Type type  ,
        DATA&&                               begin ,
        DATA&&                               end   ,
        const double                         low   )
        : m_spline    ( type , std::forward<DATA>(begin) , std::forward<DATA>(end) )
        , m_low       ( low  )
      {}
      GSLSplineInteg ( const SplineBase&              ,
                       const double           low = 0 ) ;
      GSLSplineInteg ( const GSLSplineInteg&      ) ;
      ~GSLSplineInteg () override;
    public:
      double operator() ( double a          ) const override;
      double operator() ( const Argument& x ) const override;
      unsigned int dimensionality () const override { return 1   ; }
      bool  hasAnalyticDerivative() const override { return true ; }
      Genfun::Derivative partial( unsigned int i  ) const override;
    public:
      inline   const SplineBase& spline() const { return m_spline ; }
      operator const SplineBase&       () const { return spline() ; }
    private:
      GSLSplineInteg () = delete;
      GSLSplineInteg& operator=( const GSLSplineInteg& ) = delete;
    private:
      // the actual spline function
      SplineBase m_spline ;
      double     m_low    ;
    };

  }
}

#endif // GAUDIMATH_SPLINES_H
// ============================================================================