Splines.cpp

Go to the documentation of this file.

// $Id: Splines.cpp,v 1.2 2005/11/25 10:27:03 mato Exp $
// ============================================================================
// Include files
// ============================================================================
// GSL
// ============================================================================
#include "gsl/gsl_interp.h"
#include "gsl/gsl_spline.h"
// ============================================================================
// GaudiMath
// ============================================================================
#include "GaudiMath/Splines.h"
#include "GaudiMath/GaudiMath.h"
// ============================================================================
#include <cstring>


namespace Genfun
{
  namespace GaudiMathImplementation
  {
    SplineBase::SplineBase
    ( const SplineBase::Data1D&            x    ,
      const SplineBase::Data1D&            y    ,
      const GaudiMath::Interpolation::Type type )
      : m_init      ( false    )
      , m_dim       ( x.size() )
      , m_x         ( new double[ x.size() ] )
      , m_y         ( new double[ y.size() ] )
      , m_spline    ( 0 )
      , m_accel     ( 0 )
      , m_type      ( type )
    {
#ifdef WIN32
// Disable the warning
//    C4996: 'std::copy': Function call with parameters that may be unsafe
// The parameters are checked
#pragma warning(push)
#pragma warning(disable:4996)
#endif
      std::copy( x.begin() , x.end() , m_x ) ;
      std::copy( y.begin() , y.end() , m_y ) ;
#ifdef WIN32
#pragma warning(pop)
#endif

    }
    // ========================================================================

    // ========================================================================
    // ========================================================================
    SplineBase::SplineBase
    ( const SplineBase::Data2D&            data ,
      const GaudiMath::Interpolation::Type type )
      : m_init      ( false       )
      , m_dim       ( data.size() )
      , m_x         ( new double[ data.size() ] )
      , m_y         ( new double[ data.size() ] )
      , m_spline    ( 0    )
      , m_accel     ( 0    )
      , m_type      ( type )
    {
      double*  _x = m_x ;
      double*  _y = m_y ;
      for( Data2D::const_iterator it =
             data.begin() ; data.end() != it ; ++it )
      {
        *_x = it -> first  ; ++_x ;
        *_y = it -> second ; ++_y ;
      }
    }
    // ========================================================================

    // ========================================================================
    // ========================================================================
    SplineBase::SplineBase ( const SplineBase& right )
      : m_init      ( false       )
      , m_dim       ( right.m_dim )
      , m_x         ( new double[ right.m_dim ] )
      , m_y         ( new double[ right.m_dim ] )
      , m_spline    ( 0            )
      , m_accel     ( 0            )
      , m_type      ( right.m_type )
    {
      const std::size_t size = sizeof(double) * m_dim;
      std::memcpy(m_x, right.m_x, size);
      std::memcpy(m_y, right.m_y, size);
    }
    // ========================================================================

    // ========================================================================
    // ========================================================================
    SplineBase::~SplineBase()
    {
      if ( 0 != m_spline ) { gsl_spline_free       ( m_spline ) ; }
      if ( 0 != m_accel  ) { gsl_interp_accel_free ( m_accel  ) ; }

      if ( 0 != m_x      ) { delete[] m_x      ; }
      if ( 0 != m_y      ) { delete[] m_y      ; }
    }
    // ========================================================================

    // ========================================================================
    void SplineBase::initialize() const
    {
      if ( m_init ) { return ; }                                 // RETURN

      const gsl_interp_type* T = 0 ;

      switch ( m_type )
      {
      case GaudiMath::Interpolation::Linear           :
        T = gsl_interp_linear            ; break ;
      case GaudiMath::Interpolation::Polynomial       :
        T = gsl_interp_polynomial        ; break ;
      case GaudiMath::Interpolation::Cspline          :
        T = gsl_interp_cspline           ; break ;
      case GaudiMath::Interpolation::Cspline_Periodic :
        T = gsl_interp_cspline_periodic  ; break ;
      case GaudiMath::Interpolation::Akima            :
        T = gsl_interp_akima             ; break ;
      case GaudiMath::Interpolation::Akima_Periodic   :
        T = gsl_interp_akima_periodic    ; break ;
      default :
        T = gsl_interp_cspline           ; break ;
      };

      m_spline = gsl_spline_alloc( T , m_dim ) ;

      gsl_spline_init( m_spline , m_x , m_y , m_dim ) ;

      m_accel  = gsl_interp_accel_alloc() ;

      m_init   = true ;

    }
    // ========================================================================

    // ========================================================================
    double SplineBase::eval  ( const double x ) const
    {
      if ( !m_init ) { initialize() ; }
      return gsl_spline_eval ( m_spline , x , m_accel );
    }
    // ========================================================================

    // ========================================================================
    double SplineBase::deriv  ( const double x ) const
    {
      if ( !m_init ) { initialize() ; }
      return gsl_spline_eval_deriv  ( m_spline , x , m_accel );
    }
    // ========================================================================

    // ========================================================================
    double SplineBase::deriv2 ( const double x ) const
    {
      if ( !m_init ) { initialize() ; }
      return gsl_spline_eval_deriv2 ( m_spline , x , m_accel );
    }
    // ========================================================================

    // ========================================================================
    double SplineBase::integ  ( const double a ,
                                const double b ) const
    {
      if ( !m_init ) { initialize() ; }
      return gsl_spline_eval_integ ( m_spline , a , b , m_accel ) ;
    }
    // ========================================================================


    // ========================================================================
    FUNCTION_OBJECT_IMP( GSLSpline )
    // ========================================================================

    // ========================================================================
    // ========================================================================
    GSLSpline::GSLSpline
    ( const GSLSpline::Data1D&             x    ,
      const GSLSpline::Data1D&             y    ,
      const GaudiMath::Interpolation::Type type )
      : AbsFunction ()
      , m_spline ( x , y , type )
    {}
    // ========================================================================

    // ========================================================================
    // ========================================================================
    GSLSpline::GSLSpline
    ( const GSLSpline::Data2D&             data ,
      const GaudiMath::Interpolation::Type type )
      : AbsFunction ()
      , m_spline    ( data , type )
    {}
    // ========================================================================

    // ========================================================================
    // ========================================================================
    GSLSpline::GSLSpline
    ( const SplineBase& right )
      : AbsFunction ()
      , m_spline    ( right )
    {}
    // ========================================================================

    // ========================================================================
    // ========================================================================
    GSLSpline::GSLSpline
    ( const GSLSpline& right )
      : AbsFunction ()
      , m_spline    ( right  )
    {}
    // ========================================================================

    // ========================================================================
    // ========================================================================
    GSLSpline::~GSLSpline(){}
    // ========================================================================

    // ========================================================================
    double GSLSpline::operator() (       double    x ) const
    { return m_spline.eval( x    ) ; }
    // ========================================================================
    double GSLSpline::operator() ( const Argument& x ) const
    { return m_spline.eval( x[0] ) ; }
    // ========================================================================

    // ========================================================================
    // ========================================================================
    Genfun::Derivative GSLSpline::partial( unsigned int i  ) const
    {
      if ( i >= 1 )
      {
        const AbsFunction& aux = GaudiMath::Constant( 0.0 , 1 ) ;
        return Genfun::FunctionNoop( &aux ) ;
      }
      const AbsFunction& aux = GSLSplineDeriv( *this ) ;
      return Genfun::FunctionNoop( &aux ) ;
    }
    // ========================================================================

    // ========================================================================
    FUNCTION_OBJECT_IMP( GSLSplineDeriv )
    // ========================================================================

    // ========================================================================
    // ========================================================================
    GSLSplineDeriv::GSLSplineDeriv
    ( const GSLSplineDeriv::Data1D&        x    ,
      const GSLSplineDeriv::Data1D&        y    ,
      const GaudiMath::Interpolation::Type type )
      : AbsFunction ()
      , m_spline ( x , y , type )
    {}
    // ========================================================================

    // ========================================================================
    // ========================================================================
    GSLSplineDeriv::GSLSplineDeriv
    ( const GSLSplineDeriv::Data2D&        data ,
      const GaudiMath::Interpolation::Type type )
      : AbsFunction ()
      , m_spline    ( data , type )
    {}
    // ========================================================================

    // ========================================================================
    // ========================================================================
    GSLSplineDeriv::GSLSplineDeriv
    ( const SplineBase& right )
      : AbsFunction ()
      , m_spline    ( right )
    {}
    // ========================================================================

    // ========================================================================
    // ========================================================================
    GSLSplineDeriv::GSLSplineDeriv
    ( const GSLSplineDeriv& right )
      : AbsFunction ()
      , m_spline    ( right  )
    {}
    // ========================================================================

    // ========================================================================
    // ========================================================================
    GSLSplineDeriv::~GSLSplineDeriv(){}
    // ========================================================================

    // ========================================================================
    double GSLSplineDeriv::operator() (       double    x ) const
    { return m_spline.deriv ( x    ) ; }
    // ========================================================================
    double GSLSplineDeriv::operator() ( const Argument& x ) const
    { return m_spline.deriv ( x[0] ) ; }
    // ========================================================================

    // ========================================================================
    // ========================================================================
    Genfun::Derivative GSLSplineDeriv::partial( unsigned int i  ) const
    {
      if ( i >= 1 )
      {
        const AbsFunction& aux = GaudiMath::Constant( 0.0 , 1 ) ;
        return Genfun::FunctionNoop( &aux ) ;
      }
      const AbsFunction& aux = GSLSplineDeriv2( *this ) ;
      return Genfun::FunctionNoop( &aux ) ;
    }
    // ========================================================================

    // ========================================================================
    FUNCTION_OBJECT_IMP( GSLSplineDeriv2 )
    // ========================================================================

    // ========================================================================
    // ========================================================================
    GSLSplineDeriv2::GSLSplineDeriv2
    ( const GSLSplineDeriv2::Data1D&       x    ,
      const GSLSplineDeriv2::Data1D&       y    ,
      const GaudiMath::Interpolation::Type type )
      : AbsFunction ()
      , m_spline ( x , y , type )
    {}
    // ========================================================================

    // ========================================================================
    // ========================================================================
    GSLSplineDeriv2::GSLSplineDeriv2
    ( const GSLSplineDeriv2::Data2D&       data ,
      const GaudiMath::Interpolation::Type type )
      : AbsFunction ()
      , m_spline    ( data , type )
    {}
    // ========================================================================

    // ========================================================================
    // ========================================================================
    GSLSplineDeriv2::GSLSplineDeriv2
    ( const SplineBase& right )
      : AbsFunction ()
      , m_spline    ( right )
    {}
    // ========================================================================

    // ========================================================================
    // ========================================================================
    GSLSplineDeriv2::GSLSplineDeriv2
    ( const GSLSplineDeriv2& right )
      : AbsFunction ()
      , m_spline    ( right  )
    {}
    // ========================================================================

    // ========================================================================
    // ========================================================================
    GSLSplineDeriv2::~GSLSplineDeriv2(){}
    // ========================================================================

    // ========================================================================
    double GSLSplineDeriv2::operator() (       double    x ) const
    { return m_spline.deriv2 ( x    ) ; }
    // ========================================================================
    double GSLSplineDeriv2::operator() ( const Argument& x ) const
    { return m_spline.deriv2 ( x[0] ) ; }
    // ========================================================================

    // ========================================================================
    // ========================================================================
    Genfun::Derivative GSLSplineDeriv2::partial( unsigned int i  ) const
    {
      if ( i >= 1 )
      {
        const AbsFunction& aux = GaudiMath::Constant( 0.0 , 1 ) ;
        return Genfun::FunctionNoop( &aux ) ;
      }
      const AbsFunction& aux = GaudiMath::Derivative( *this , i ) ;
      return Genfun::FunctionNoop( &aux ) ;
    }
    // ========================================================================


    // ========================================================================
    FUNCTION_OBJECT_IMP( GSLSplineInteg )
    // ========================================================================

    // ========================================================================
    // ========================================================================
    GSLSplineInteg::GSLSplineInteg
    ( const GSLSplineInteg::Data1D&        x     ,
      const GSLSplineInteg::Data1D&        y     ,
      const GaudiMath::Interpolation::Type type  ,
      const double                         low   )
      : AbsFunction (      )
      , m_spline    ( x    , y , type )
      , m_low       ( low  )
    {}
    // ========================================================================

    // ========================================================================
    // ========================================================================
    GSLSplineInteg::GSLSplineInteg
    ( const GSLSplineInteg::Data2D&        data  ,
      const GaudiMath::Interpolation::Type type  ,
      const double                         low   )
      : AbsFunction ()
      , m_spline    ( data , type )
      , m_low       ( low  )
    {}
    // ========================================================================

    // ========================================================================
    // ========================================================================
    GSLSplineInteg::GSLSplineInteg
    ( const SplineBase& right ,
      const double      low   )
      : AbsFunction ()
      , m_spline    ( right )
      , m_low       ( low   )
    {}
    // ========================================================================

    // ========================================================================
    // ========================================================================
    GSLSplineInteg::GSLSplineInteg
    ( const GSLSplineInteg& right )
      : AbsFunction ()
      , m_spline    ( right       )
      , m_low       ( right.m_low )
    {}
    // ========================================================================

    // ========================================================================
    // ========================================================================
    GSLSplineInteg::~GSLSplineInteg(){}
    // ========================================================================

    // ========================================================================
    double GSLSplineInteg::operator() (       double    x ) const
    { return m_spline.integ ( m_low , x    ) ; }
    // ========================================================================
    double GSLSplineInteg::operator() ( const Argument& x ) const
    { return m_spline.integ ( m_low , x[0] ) ; }
    // ========================================================================

    // ========================================================================
    // ========================================================================
    Genfun::Derivative GSLSplineInteg::partial( unsigned int i  ) const
    {
      if ( i >= 1 )
      {
        const AbsFunction& aux = GaudiMath::Constant( 0.0 , 1 ) ;
        return Genfun::FunctionNoop( &aux ) ;
      }
      const AbsFunction& aux = GSLSpline( *this ) ;
      return Genfun::FunctionNoop( &aux ) ;
    }
    // ========================================================================

  }
}