The Gaudi Framework  v30r3 (a5ef0a68)
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 {
24  namespace GaudiMathImplementation
25  {
31  SplineBase::SplineBase( const SplineBase::Data1D& x, const SplineBase::Data1D& y,
32  const GaudiMath::Interpolation::Type type )
33  : m_x( x ), m_y( y ), m_type( type )
34  {
35  }
36  // ========================================================================
37 
38  // ========================================================================
43  // ========================================================================
44  SplineBase::SplineBase( const SplineBase::Data2D& data, const GaudiMath::Interpolation::Type type ) : m_type( type )
45  {
46  m_x.reserve( data.size() );
47  m_y.reserve( data.size() );
48  for ( const auto& i : data ) {
49  m_x.push_back( i.first );
50  m_y.push_back( i.second );
51  }
52  }
53  // ========================================================================
54 
55  // ========================================================================
56  void SplineBase::initialize() const
57  {
58  if ( m_spline ) {
59  return;
60  } // RETURN
61 
62  const gsl_interp_type* T = nullptr;
63 
64  switch ( m_type ) {
65  case GaudiMath::Interpolation::Linear:
66  T = gsl_interp_linear;
67  break;
68  case GaudiMath::Interpolation::Polynomial:
69  T = gsl_interp_polynomial;
70  break;
71  case GaudiMath::Interpolation::Cspline:
72  T = gsl_interp_cspline;
73  break;
74  case GaudiMath::Interpolation::Cspline_Periodic:
75  T = gsl_interp_cspline_periodic;
76  break;
77  case GaudiMath::Interpolation::Akima:
78  T = gsl_interp_akima;
79  break;
80  case GaudiMath::Interpolation::Akima_Periodic:
81  T = gsl_interp_akima_periodic;
82  break;
83  default:
84  T = gsl_interp_cspline;
85  break;
86  };
87 
88  m_spline.reset( gsl_spline_alloc( T, m_x.size() ) );
89  gsl_spline_init( m_spline.get(), m_x.data(), m_y.data(), m_x.size() );
90  m_accel.reset( gsl_interp_accel_alloc() );
91 
92  if ( !m_accel ) m_spline.reset();
93  }
94  // ========================================================================
95 
96  // ========================================================================
97  double SplineBase::eval( const double x ) const
98  {
99  if ( !m_spline ) {
100  initialize();
101  }
102  return gsl_spline_eval( m_spline.get(), x, m_accel.get() );
103  }
104  // ========================================================================
105 
106  // ========================================================================
107  double SplineBase::deriv( const double x ) const
108  {
109  if ( !m_spline ) {
110  initialize();
111  }
112  return gsl_spline_eval_deriv( m_spline.get(), x, m_accel.get() );
113  }
114  // ========================================================================
115 
116  // ========================================================================
117  double SplineBase::deriv2( const double x ) const
118  {
119  if ( !m_spline ) {
120  initialize();
121  }
122  return gsl_spline_eval_deriv2( m_spline.get(), x, m_accel.get() );
123  }
124  // ========================================================================
125 
126  // ========================================================================
127  double SplineBase::integ( const double a, const double b ) const
128  {
129  if ( !m_spline ) {
130  initialize();
131  }
132  return gsl_spline_eval_integ( m_spline.get(), a, b, m_accel.get() );
133  }
134  // ========================================================================
135 
136  // ========================================================================
157  // ========================================================================
158  GSLSpline::GSLSpline( const GSLSpline::Data1D& x, const GSLSpline::Data1D& y,
159  const GaudiMath::Interpolation::Type type )
160  : AbsFunction(), m_spline( x, y, type )
161  {
162  }
163  // ========================================================================
164 
165  // ========================================================================
184  // ========================================================================
185  GSLSpline::GSLSpline( const GSLSpline::Data2D& data, const GaudiMath::Interpolation::Type type )
186  : AbsFunction(), m_spline( data, type )
187  {
188  }
189  // ========================================================================
190 
191  // ========================================================================
193  // ========================================================================
194  GSLSpline::GSLSpline( const SplineBase& right ) : AbsFunction(), m_spline( right ) {}
195  // ========================================================================
196 
197  // ========================================================================
198  double GSLSpline::operator()( double x ) const { return m_spline.eval( x ); }
199  // ========================================================================
200  double GSLSpline::operator()( const Argument& x ) const { return m_spline.eval( x[0] ); }
201  // ========================================================================
202 
203  // ========================================================================
205  // ========================================================================
206  Genfun::Derivative GSLSpline::partial( unsigned int i ) const
207  {
208  if ( i >= 1 ) {
209  const AbsFunction& aux = GaudiMath::Constant( 0.0, 1 );
210  return Genfun::FunctionNoop( &aux );
211  }
212  const AbsFunction& aux = GSLSplineDeriv( *this );
213  return Genfun::FunctionNoop( &aux );
214  }
215  // ========================================================================
216 
217  // ========================================================================
238  // ========================================================================
239  GSLSplineDeriv::GSLSplineDeriv( const GSLSplineDeriv::Data1D& x, const GSLSplineDeriv::Data1D& y,
240  const GaudiMath::Interpolation::Type type )
241  : AbsFunction(), m_spline( x, y, type )
242  {
243  }
244  // ========================================================================
245 
246  // ========================================================================
265  // ========================================================================
266  GSLSplineDeriv::GSLSplineDeriv( const GSLSplineDeriv::Data2D& data, const GaudiMath::Interpolation::Type type )
267  : AbsFunction(), m_spline( data, type )
268  {
269  }
270  // ========================================================================
271 
272  // ========================================================================
274  // ========================================================================
275  GSLSplineDeriv::GSLSplineDeriv( const SplineBase& right ) : AbsFunction(), m_spline( right ) {}
276  // ========================================================================
277 
278  // ========================================================================
280  // ========================================================================
281  GSLSplineDeriv::GSLSplineDeriv( const GSLSplineDeriv& right ) : AbsFunction(), m_spline( right ) {}
282  // ========================================================================
283 
284  // ========================================================================
285  double GSLSplineDeriv::operator()( double x ) const { return m_spline.deriv( x ); }
286  // ========================================================================
287  double GSLSplineDeriv::operator()( const Argument& x ) const { return m_spline.deriv( x[0] ); }
288  // ========================================================================
289 
290  // ========================================================================
292  // ========================================================================
293  Genfun::Derivative GSLSplineDeriv::partial( unsigned int i ) const
294  {
295  if ( i >= 1 ) {
296  const AbsFunction& aux = GaudiMath::Constant( 0.0, 1 );
297  return Genfun::FunctionNoop( &aux );
298  }
299  const AbsFunction& aux = GSLSplineDeriv2( *this );
300  return Genfun::FunctionNoop( &aux );
301  }
302  // ========================================================================
303 
304  // ========================================================================
325  // ========================================================================
326  GSLSplineDeriv2::GSLSplineDeriv2( const GSLSplineDeriv2::Data1D& x, const GSLSplineDeriv2::Data1D& y,
327  const GaudiMath::Interpolation::Type type )
328  : AbsFunction(), m_spline( x, y, type )
329  {
330  }
331  // ========================================================================
332 
333  // ========================================================================
352  // ========================================================================
353  GSLSplineDeriv2::GSLSplineDeriv2( const GSLSplineDeriv2::Data2D& data, const GaudiMath::Interpolation::Type type )
354  : AbsFunction(), m_spline( data, type )
355  {
356  }
357  // ========================================================================
358 
359  // ========================================================================
361  // ========================================================================
362  GSLSplineDeriv2::GSLSplineDeriv2( const SplineBase& right ) : AbsFunction(), m_spline( right ) {}
363  // ========================================================================
364 
365  // ========================================================================
366  double GSLSplineDeriv2::operator()( double x ) const { return m_spline.deriv2( x ); }
367  // ========================================================================
368  double GSLSplineDeriv2::operator()( const Argument& x ) const { return m_spline.deriv2( x[0] ); }
369  // ========================================================================
370 
371  // ========================================================================
373  // ========================================================================
374  Genfun::Derivative GSLSplineDeriv2::partial( unsigned int i ) const
375  {
376  if ( i >= 1 ) {
377  const AbsFunction& aux = GaudiMath::Constant( 0.0, 1 );
378  return Genfun::FunctionNoop( &aux );
379  }
380  const AbsFunction& aux = GaudiMath::Derivative( *this, i );
381  return Genfun::FunctionNoop( &aux );
382  }
383  // ========================================================================
384 
385  // ========================================================================
406  // ========================================================================
407  GSLSplineInteg::GSLSplineInteg( const GSLSplineInteg::Data1D& x, const GSLSplineInteg::Data1D& y,
408  const GaudiMath::Interpolation::Type type, const double low )
409  : AbsFunction(), m_spline( x, y, type ), m_low( low )
410  {
411  }
412  // ========================================================================
413 
414  // ========================================================================
433  // ========================================================================
434  GSLSplineInteg::GSLSplineInteg( const GSLSplineInteg::Data2D& data, const GaudiMath::Interpolation::Type type,
435  const double low )
436  : AbsFunction(), m_spline( data, type ), m_low( low )
437  {
438  }
439  // ========================================================================
440 
441  // ========================================================================
443  // ========================================================================
444  GSLSplineInteg::GSLSplineInteg( const SplineBase& right, const double low )
445  : AbsFunction(), m_spline( right ), m_low( low )
446  {
447  }
448  // ========================================================================
449 
450  // ========================================================================
451  double GSLSplineInteg::operator()( double x ) const { return m_spline.integ( m_low, x ); }
452  // ========================================================================
453  double GSLSplineInteg::operator()( const Argument& x ) const { return m_spline.integ( m_low, x[0] ); }
454  // ========================================================================
455 
456  // ========================================================================
458  // ========================================================================
459  Genfun::Derivative GSLSplineInteg::partial( unsigned int i ) const
460  {
461  if ( i >= 1 ) {
462  const AbsFunction& aux = GaudiMath::Constant( 0.0, 1 );
463  return Genfun::FunctionNoop( &aux );
464  }
465  const AbsFunction& aux = GSLSpline( *this );
466  return Genfun::FunctionNoop( &aux );
467  }
468  // ========================================================================
469  }
470 }
Genfun::GaudiMathImplementation::SplineBase::integ
double integ(const double a, const double b) const
evaluate the integral on [a,b]
Definition: Splines.cpp:127
Genfun::GaudiMathImplementation::GSLSplineDeriv::operator()
double operator()(double a) const override
main method: evaluate the function
Definition: Splines.cpp:285
Genfun::GaudiMathImplementation::SplineBase::SplineBase
SplineBase(const Data1D &x, const Data1D &y, const GaudiMath::Interpolation::Type type)
constructor from vectors and type
Definition: Splines.cpp:31
Genfun::GaudiMathImplementation::SplineBase::eval
double eval(const double x) const
evaluate the function
Definition: Splines.cpp:97
Genfun::GaudiMathImplementation::GSLSplineDeriv2::operator()
double operator()(double a) const override
main method: evaluate the function
Definition: Splines.cpp:366
Genfun::GaudiMathImplementation::GSLSplineInteg::operator()
double operator()(double a) const override
main method: evaluate the function
Definition: Splines.cpp:451
Splines.h
Genfun::GaudiMathImplementation::GSLSplineInteg::partial
Genfun::Derivative partial(unsigned int i) const override
Derivatives.
Definition: Splines.cpp:459
Genfun::GaudiMathImplementation::SplineBase::m_spline
std::unique_ptr< gsl_spline, details::gsl_deleter > m_spline
Definition: Splines.h:138
GaudiMath::Derivative
Genfun::GaudiMathImplementation::NumericalDerivative Derivative
Definition: GaudiMath.h:27
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:140
Genfun::GaudiMathImplementation::GSLSplineDeriv::partial
Genfun::Derivative partial(unsigned int i) const override
Derivatives.
Definition: Splines.cpp:293
Genfun::GaudiMathImplementation::SplineBase::deriv2
double deriv2(const double x) const
evaluate the second derivative
Definition: Splines.cpp:117
Genfun::GaudiMathImplementation::SplineBase::deriv
double deriv(const double x) const
evaluate the first derivative
Definition: Splines.cpp:107
Genfun::GaudiMathImplementation::GSLSplineDeriv2::GSLSplineDeriv2
GSLSplineDeriv2(const Data1D &x, const Data1D &y, const GaudiMath::Interpolation::Type type)
constructor from vectors and type
Definition: Splines.cpp:326
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:198
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:158
Genfun::GaudiMathImplementation::GSLSpline::partial
Genfun::Derivative partial(unsigned int i) const override
Derivatives.
Definition: Splines.cpp:206
Genfun::GaudiMathImplementation::SplineBase::m_accel
std::unique_ptr< gsl_interp_accel, details::gsl_deleter > m_accel
transient
Definition: Splines.h:139
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:407
GaudiMath::Constant
Genfun::GaudiMathImplementation::Constant Constant
Definition: GaudiMath.h:26
GaudiMath.h
Genfun::GaudiMathImplementation::GSLSplineDeriv2::partial
Genfun::Derivative partial(unsigned int i) const override
Derivatives.
Definition: Splines.cpp:374
std::vector::reserve
T reserve(T...args)