Loading [MathJax]/extensions/tex2jax.js
The Gaudi Framework  v31r0 (aeb156f0)
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
Splines.cpp
Go to the documentation of this file.
1 // ============================================================================
2 // Include files
3 // ============================================================================
4 // GSL
5 // ============================================================================
6 #include "gsl/gsl_interp.h"
7 #include "gsl/gsl_spline.h"
8 // ============================================================================
9 // GaudiMath
10 // ============================================================================
11 #include "GaudiMath/GaudiMath.h"
12 #include "GaudiMath/Splines.h"
13 // ============================================================================
14 #include <cstring>
15 
22 namespace Genfun {
23  namespace GaudiMathImplementation {
29  SplineBase::SplineBase( const SplineBase::Data1D& x, const SplineBase::Data1D& y,
30  const GaudiMath::Interpolation::Type type )
31  : m_x( x ), m_y( y ), m_type( type ) {}
32  // ========================================================================
33 
34  // ========================================================================
39  // ========================================================================
40  SplineBase::SplineBase( const SplineBase::Data2D& data, const GaudiMath::Interpolation::Type type )
41  : m_type( type ) {
42  m_x.reserve( data.size() );
43  m_y.reserve( data.size() );
44  for ( const auto& i : data ) {
45  m_x.push_back( i.first );
46  m_y.push_back( i.second );
47  }
48  }
49  // ========================================================================
50 
51  // ========================================================================
52  void SplineBase::initialize() const {
53  if ( m_spline ) { return; } // RETURN
54 
55  const gsl_interp_type* T = nullptr;
56 
57  switch ( m_type ) {
58  case GaudiMath::Interpolation::Linear:
59  T = gsl_interp_linear;
60  break;
61  case GaudiMath::Interpolation::Polynomial:
62  T = gsl_interp_polynomial;
63  break;
64  case GaudiMath::Interpolation::Cspline:
65  T = gsl_interp_cspline;
66  break;
67  case GaudiMath::Interpolation::Cspline_Periodic:
68  T = gsl_interp_cspline_periodic;
69  break;
70  case GaudiMath::Interpolation::Akima:
71  T = gsl_interp_akima;
72  break;
73  case GaudiMath::Interpolation::Akima_Periodic:
74  T = gsl_interp_akima_periodic;
75  break;
76  default:
77  T = gsl_interp_cspline;
78  break;
79  };
80 
81  m_spline.reset( gsl_spline_alloc( T, m_x.size() ) );
82  gsl_spline_init( m_spline.get(), m_x.data(), m_y.data(), m_x.size() );
83  m_accel.reset( gsl_interp_accel_alloc() );
84 
85  if ( !m_accel ) m_spline.reset();
86  }
87  // ========================================================================
88 
89  // ========================================================================
90  double SplineBase::eval( const double x ) const {
91  if ( !m_spline ) { initialize(); }
92  return gsl_spline_eval( m_spline.get(), x, m_accel.get() );
93  }
94  // ========================================================================
95 
96  // ========================================================================
97  double SplineBase::deriv( const double x ) const {
98  if ( !m_spline ) { initialize(); }
99  return gsl_spline_eval_deriv( m_spline.get(), x, m_accel.get() );
100  }
101  // ========================================================================
102 
103  // ========================================================================
104  double SplineBase::deriv2( const double x ) const {
105  if ( !m_spline ) { initialize(); }
106  return gsl_spline_eval_deriv2( m_spline.get(), x, m_accel.get() );
107  }
108  // ========================================================================
109 
110  // ========================================================================
111  double SplineBase::integ( const double a, const double b ) const {
112  if ( !m_spline ) { initialize(); }
113  return gsl_spline_eval_integ( m_spline.get(), a, b, m_accel.get() );
114  }
115  // ========================================================================
116 
117  // ========================================================================
138  // ========================================================================
139  GSLSpline::GSLSpline( const GSLSpline::Data1D& x, const GSLSpline::Data1D& y,
140  const GaudiMath::Interpolation::Type type )
141  : AbsFunction(), m_spline( x, y, type ) {}
142  // ========================================================================
143 
144  // ========================================================================
163  // ========================================================================
164  GSLSpline::GSLSpline( const GSLSpline::Data2D& data, const GaudiMath::Interpolation::Type type )
165  : AbsFunction(), m_spline( data, type ) {}
166  // ========================================================================
167 
168  // ========================================================================
170  // ========================================================================
171  GSLSpline::GSLSpline( const SplineBase& right ) : AbsFunction(), m_spline( right ) {}
172  // ========================================================================
173 
174  // ========================================================================
175  double GSLSpline::operator()( double x ) const { return m_spline.eval( x ); }
176  // ========================================================================
177  double GSLSpline::operator()( const Argument& x ) const { return m_spline.eval( x[0] ); }
178  // ========================================================================
179 
180  // ========================================================================
182  // ========================================================================
183  Genfun::Derivative GSLSpline::partial( unsigned int i ) const {
184  if ( i >= 1 ) {
185  const AbsFunction& aux = GaudiMath::Constant( 0.0, 1 );
186  return Genfun::FunctionNoop( &aux );
187  }
188  const AbsFunction& aux = GSLSplineDeriv( *this );
189  return Genfun::FunctionNoop( &aux );
190  }
191  // ========================================================================
192 
193  // ========================================================================
214  // ========================================================================
215  GSLSplineDeriv::GSLSplineDeriv( const GSLSplineDeriv::Data1D& x, const GSLSplineDeriv::Data1D& y,
216  const GaudiMath::Interpolation::Type type )
217  : AbsFunction(), m_spline( x, y, type ) {}
218  // ========================================================================
219 
220  // ========================================================================
239  // ========================================================================
240  GSLSplineDeriv::GSLSplineDeriv( const GSLSplineDeriv::Data2D& data, const GaudiMath::Interpolation::Type type )
241  : AbsFunction(), m_spline( data, type ) {}
242  // ========================================================================
243 
244  // ========================================================================
246  // ========================================================================
247  GSLSplineDeriv::GSLSplineDeriv( const SplineBase& right ) : AbsFunction(), m_spline( right ) {}
248  // ========================================================================
249 
250  // ========================================================================
252  // ========================================================================
253  GSLSplineDeriv::GSLSplineDeriv( const GSLSplineDeriv& right ) : AbsFunction(), m_spline( right ) {}
254  // ========================================================================
255 
256  // ========================================================================
257  double GSLSplineDeriv::operator()( double x ) const { return m_spline.deriv( x ); }
258  // ========================================================================
259  double GSLSplineDeriv::operator()( const Argument& x ) const { return m_spline.deriv( x[0] ); }
260  // ========================================================================
261 
262  // ========================================================================
264  // ========================================================================
265  Genfun::Derivative GSLSplineDeriv::partial( unsigned int i ) const {
266  if ( i >= 1 ) {
267  const AbsFunction& aux = GaudiMath::Constant( 0.0, 1 );
268  return Genfun::FunctionNoop( &aux );
269  }
270  const AbsFunction& aux = GSLSplineDeriv2( *this );
271  return Genfun::FunctionNoop( &aux );
272  }
273  // ========================================================================
274 
275  // ========================================================================
296  // ========================================================================
297  GSLSplineDeriv2::GSLSplineDeriv2( const GSLSplineDeriv2::Data1D& x, const GSLSplineDeriv2::Data1D& y,
298  const GaudiMath::Interpolation::Type type )
299  : AbsFunction(), m_spline( x, y, type ) {}
300  // ========================================================================
301 
302  // ========================================================================
321  // ========================================================================
322  GSLSplineDeriv2::GSLSplineDeriv2( const GSLSplineDeriv2::Data2D& data, const GaudiMath::Interpolation::Type type )
323  : AbsFunction(), m_spline( data, type ) {}
324  // ========================================================================
325 
326  // ========================================================================
328  // ========================================================================
329  GSLSplineDeriv2::GSLSplineDeriv2( const SplineBase& right ) : AbsFunction(), m_spline( right ) {}
330  // ========================================================================
331 
332  // ========================================================================
333  double GSLSplineDeriv2::operator()( double x ) const { return m_spline.deriv2( x ); }
334  // ========================================================================
335  double GSLSplineDeriv2::operator()( const Argument& x ) const { return m_spline.deriv2( x[0] ); }
336  // ========================================================================
337 
338  // ========================================================================
340  // ========================================================================
341  Genfun::Derivative GSLSplineDeriv2::partial( unsigned int i ) const {
342  if ( i >= 1 ) {
343  const AbsFunction& aux = GaudiMath::Constant( 0.0, 1 );
344  return Genfun::FunctionNoop( &aux );
345  }
346  const AbsFunction& aux = GaudiMath::Derivative( *this, i );
347  return Genfun::FunctionNoop( &aux );
348  }
349  // ========================================================================
350 
351  // ========================================================================
372  // ========================================================================
373  GSLSplineInteg::GSLSplineInteg( const GSLSplineInteg::Data1D& x, const GSLSplineInteg::Data1D& y,
374  const GaudiMath::Interpolation::Type type, const double low )
375  : AbsFunction(), m_spline( x, y, type ), m_low( low ) {}
376  // ========================================================================
377 
378  // ========================================================================
397  // ========================================================================
398  GSLSplineInteg::GSLSplineInteg( const GSLSplineInteg::Data2D& data, const GaudiMath::Interpolation::Type type,
399  const double low )
400  : AbsFunction(), m_spline( data, type ), m_low( low ) {}
401  // ========================================================================
402 
403  // ========================================================================
405  // ========================================================================
406  GSLSplineInteg::GSLSplineInteg( const SplineBase& right, const double low )
407  : AbsFunction(), m_spline( right ), m_low( low ) {}
408  // ========================================================================
409 
410  // ========================================================================
411  double GSLSplineInteg::operator()( double x ) const { return m_spline.integ( m_low, x ); }
412  // ========================================================================
413  double GSLSplineInteg::operator()( const Argument& x ) const { return m_spline.integ( m_low, x[0] ); }
414  // ========================================================================
415 
416  // ========================================================================
418  // ========================================================================
419  Genfun::Derivative GSLSplineInteg::partial( unsigned int i ) const {
420  if ( i >= 1 ) {
421  const AbsFunction& aux = GaudiMath::Constant( 0.0, 1 );
422  return Genfun::FunctionNoop( &aux );
423  }
424  const AbsFunction& aux = GSLSpline( *this );
425  return Genfun::FunctionNoop( &aux );
426  }
427  // ========================================================================
428  } // namespace GaudiMathImplementation
429 } // namespace Genfun
Genfun::GaudiMathImplementation::SplineBase::integ
double integ(const double a, const double b) const
evaluate the integral on [a,b]
Definition: Splines.cpp:111
Genfun::GaudiMathImplementation::GSLSplineDeriv::operator()
double operator()(double a) const override
main method: evaluate the function
Definition: Splines.cpp:257
Genfun::GaudiMathImplementation::SplineBase::SplineBase
SplineBase(const Data1D &x, const Data1D &y, const GaudiMath::Interpolation::Type type)
constructor from vectors and type
Definition: Splines.cpp:29
Genfun::GaudiMathImplementation::SplineBase::eval
double eval(const double x) const
evaluate the function
Definition: Splines.cpp:90
Genfun::GaudiMathImplementation::GSLSplineDeriv2::operator()
double operator()(double a) const override
main method: evaluate the function
Definition: Splines.cpp:333
Genfun::GaudiMathImplementation::GSLSplineInteg::operator()
double operator()(double a) const override
main method: evaluate the function
Definition: Splines.cpp:411
Splines.h
Genfun::GaudiMathImplementation::GSLSplineInteg::partial
Genfun::Derivative partial(unsigned int i) const override
Derivatives.
Definition: Splines.cpp:419
Genfun::GaudiMathImplementation::SplineBase::m_spline
std::unique_ptr< gsl_spline, details::gsl_deleter > m_spline
Definition: Splines.h:131
GaudiMath::Derivative
Genfun::GaudiMathImplementation::NumericalDerivative Derivative
Definition: GaudiMath.h:26
std::vector::push_back
T push_back(T...args)
std::vector::data
T data(T...args)
Genfun::GaudiMathImplementation::SplineBase::m_type
GaudiMath::Interpolation::Type m_type
transient
Definition: Splines.h:133
Genfun::GaudiMathImplementation::GSLSplineDeriv::partial
Genfun::Derivative partial(unsigned int i) const override
Derivatives.
Definition: Splines.cpp:265
Genfun::GaudiMathImplementation::SplineBase::deriv2
double deriv2(const double x) const
evaluate the second derivative
Definition: Splines.cpp:104
Genfun::GaudiMathImplementation::SplineBase::deriv
double deriv(const double x) const
evaluate the first derivative
Definition: Splines.cpp:97
Genfun::GaudiMathImplementation::GSLSplineDeriv2::GSLSplineDeriv2
GSLSplineDeriv2(const Data1D &x, const Data1D &y, const GaudiMath::Interpolation::Type type)
constructor from vectors and type
Definition: Splines.cpp:297
Genfun::GaudiMathImplementation::GSLSplineDeriv::GSLSplineDeriv
GSLSplineDeriv()
default construtor is desabled ;
std::unique_ptr::reset
T reset(T...args)
Genfun::GaudiMathImplementation::GSLSpline::operator()
double operator()(double a) const override
main methgod: evaluate teh function
Definition: Splines.cpp:175
std::unique_ptr::get
T get(T...args)
std::vector::size
T size(T...args)
Genfun::GaudiMathImplementation::GSLSpline::GSLSpline
GSLSpline(const Data1D &x, const Data1D &y, const GaudiMath::Interpolation::Type type)
constructor from vectors and type
Definition: Splines.cpp:139
Genfun::GaudiMathImplementation::GSLSpline::partial
Genfun::Derivative partial(unsigned int i) const override
Derivatives.
Definition: Splines.cpp:183
Genfun::GaudiMathImplementation::SplineBase::m_accel
std::unique_ptr< gsl_interp_accel, details::gsl_deleter > m_accel
transient
Definition: Splines.h:132
Genfun
CLHEP.
Definition: IEqSolver.h:13
Genfun::GaudiMathImplementation::GSLSplineInteg::GSLSplineInteg
GSLSplineInteg(const Data1D &x, const Data1D &y, const GaudiMath::Interpolation::Type type, const double low=0)
constructor from vectors and type
Definition: Splines.cpp:373
GaudiMath::Constant
Genfun::GaudiMathImplementation::Constant Constant
Definition: GaudiMath.h:25
GaudiMath.h
Genfun::GaudiMathImplementation::GSLSplineDeriv2::partial
Genfun::Derivative partial(unsigned int i) const override
Derivatives.
Definition: Splines.cpp:341
std::vector::reserve
T reserve(T...args)