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"

#ifdef __clang__
#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Winconsistent-missing-override"
#endif

namespace Genfun
{
  namespace GaudiMathImplementation
  {
      namespace details {
      struct gsl_deleter {
          void operator()(gsl_spline* p) { if (p) gsl_spline_free(p); }
          void operator()(gsl_interp_accel* p) { if (p) gsl_interp_accel_free(p); }
      }; 
      }

    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_type      ( type )
      {
        auto size = std::distance(begin_x,end_x);
        m_x.reserve(size);
        std::copy_n ( begin_x , size , std::back_inserter(m_x) ) ;
        m_y.reserve(size);
        std::copy_n ( begin_y , size , std::back_inserter(m_y) ) ;
      }
      template <class DATA>
      SplineBase
      ( const GaudiMath::Interpolation::Type type  ,
        DATA                                 begin ,
        DATA                                 end   )
        : m_type      ( type )
      {
        m_x.reserve(end-begin);
        m_y.reserve(end-begin);
        for ( ; begin != end ; ++begin ) {
          m_x.push_back( begin -> first ) ;
          m_y.push_back( begin -> second) ;
        };
      }

    SplineBase( SplineBase&& ) = default;
    SplineBase& operator=( SplineBase&& rhs ) = default;
    SplineBase( const SplineBase& rhs)
      : m_x         ( rhs.m_x )
      , m_y         ( rhs.m_y )
      , m_type      ( rhs.m_type )
      // note that we do NOT copy m_spline, m_accel!
    {}
    SplineBase& operator=( const SplineBase& rhs ) { *this = SplineBase(rhs); return *this; }

    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:
      std::vector<double>            m_x      ;
      std::vector<double>            m_y      ;
      mutable std::unique_ptr<gsl_spline,details::gsl_deleter> m_spline;
      mutable std::unique_ptr<gsl_interp_accel,details::gsl_deleter> m_accel;
      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 )
        : m_spline( type , begin_x , end_x , begin_y )
      {}
      template <class DATA>
      GSLSpline( const GaudiMath::Interpolation::Type type  ,
                DATA                                 begin ,
                DATA                                 end   )
      : m_spline( type , begin , end )
      {}
      GSLSpline ( const SplineBase& ) ;
    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() ; }
      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& ) ;
    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& ) = default;
    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& 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&      ) = default;
    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& operator=( const GSLSplineInteg& ) = delete;
    private:
      // the actual spline function
      SplineBase m_spline ;
      double     m_low    ;
    };

  }
}

#ifdef __clang__
#pragma clang diagnostic pop
#endif

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