Splines.h

Go to the documentation of this file.

#ifndef GAUDIMATH_SPLINES_H
#define GAUDIMATH_SPLINES_H 1
// ============================================================================
// Include files
// ============================================================================
// STD/STL
// ============================================================================
#include <algorithm>
#include <memory>
#include <utility>
#include <vector>
// ============================================================================
// from CLHEP
// ============================================================================
#include "CLHEP/GenericFunctions/AbsFunction.hh"
#include "CLHEP/GenericFunctions/Argument.hh"
#include "CLHEP/GenericFunctions/GenericFunctions.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 );
        }
      }; 
    }    // namespace details

    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:
      // mandatory macro from CLHEP/GenericFunctions
      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;
    };
    // mandatory macro from CLHEP/GenericFunctions
    FUNCTION_OBJECT_IMP( GSLSpline )

    class GAUDI_API GSLSplineDeriv : public AbsFunction {
    public:
      typedef SplineBase::Data1D Data1D;
      typedef SplineBase::Data2D Data2D;

    public:
      // mandatory macro from CLHEP/GenericFunctions
      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;
    };
    // mandatory macro from CLHEP/GenericFunctions
    FUNCTION_OBJECT_IMP( GSLSplineDeriv )

    class GAUDI_API GSLSplineDeriv2 : public AbsFunction {
    public:
      typedef SplineBase::Data1D Data1D;
      typedef SplineBase::Data2D Data2D;

    public:
      // mandatory macro from CLHEP/GenericFunctions
      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;
    };
    // mandatory macro from CLHEP/GenericFunctions
    FUNCTION_OBJECT_IMP( GSLSplineDeriv2 )

    class GAUDI_API GSLSplineInteg : public AbsFunction {
    public:
      typedef SplineBase::Data1D Data1D;
      typedef SplineBase::Data2D Data2D;

    public:
      // mandatory macro from CLHEP/GenericFunctions
      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;
    };
    // mandatory macro from CLHEP/GenericFunctions
    FUNCTION_OBJECT_IMP( GSLSplineInteg )
  } // namespace GaudiMathImplementation
} // namespace Genfun

#ifdef __clang__
#  pragma clang diagnostic pop
#endif

#endif // GAUDIMATH_SPLINES_H