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/Splines.h"
12 #include "GaudiMath/GaudiMath.h"
13 // ============================================================================
14 #include <cstring>
15 
16 
23 namespace Genfun
24 {
25  namespace GaudiMathImplementation
26  {
32  SplineBase::SplineBase
33  ( const SplineBase::Data1D& x ,
34  const SplineBase::Data1D& y ,
35  const GaudiMath::Interpolation::Type type )
36  : m_init ( false )
37  , m_dim ( x.size() )
38  , m_x ( new double[ x.size() ] )
39  , m_y ( new double[ y.size() ] )
40  , m_spline ( nullptr )
41  , m_accel ( nullptr )
42  , m_type ( type )
43  {
44 #ifdef WIN32
45 // Disable the warning
46 // C4996: 'std::copy': Function call with parameters that may be unsafe
47 // The parameters are checked
48 #pragma warning(push)
49 #pragma warning(disable:4996)
50 #endif
51  std::copy( x.begin() , x.end() , m_x.get() ) ;
52  std::copy( y.begin() , y.end() , m_y.get() ) ;
53 #ifdef WIN32
54 #pragma warning(pop)
55 #endif
56 
57  }
58  // ========================================================================
59 
60  // ========================================================================
65  // ========================================================================
66  SplineBase::SplineBase
67  ( const SplineBase::Data2D& data ,
68  const GaudiMath::Interpolation::Type type )
69  : m_init ( false )
70  , m_dim ( data.size() )
71  , m_x ( new double[ data.size() ] )
72  , m_y ( new double[ data.size() ] )
73  , m_type ( type )
74  {
75  double* _x = m_x.get() ;
76  double* _y = m_y.get() ;
77  for( const auto& i : data )
78  {
79  *_x = i.first ; ++_x ;
80  *_y = i.second ; ++_y ;
81  }
82  }
83  // ========================================================================
84 
85  // ========================================================================
87  // ========================================================================
88  SplineBase::SplineBase ( const SplineBase& right )
89  : m_init ( false )
90  , m_dim ( right.m_dim )
91  , m_x ( new double[ right.m_dim ] )
92  , m_y ( new double[ right.m_dim ] )
93  , m_type ( right.m_type )
94  {
95  std::copy_n(right.m_x.get(),m_dim,m_x.get());
96  std::copy_n(right.m_y.get(),m_dim,m_y.get());
97  }
98  // ========================================================================
99 
100  // ========================================================================
102  // ========================================================================
103  SplineBase::~SplineBase()
104  {
105  if ( m_spline ) { gsl_spline_free ( m_spline ) ; }
106  if ( m_accel ) { gsl_interp_accel_free ( m_accel ) ; }
107  }
108  // ========================================================================
109 
110  // ========================================================================
111  void SplineBase::initialize() const
112  {
113  if ( m_init ) { return ; } // RETURN
114 
115  const gsl_interp_type* T = nullptr ;
116 
117  switch ( m_type )
118  {
119  case GaudiMath::Interpolation::Linear :
120  T = gsl_interp_linear ; break ;
121  case GaudiMath::Interpolation::Polynomial :
122  T = gsl_interp_polynomial ; break ;
123  case GaudiMath::Interpolation::Cspline :
124  T = gsl_interp_cspline ; break ;
125  case GaudiMath::Interpolation::Cspline_Periodic :
126  T = gsl_interp_cspline_periodic ; break ;
127  case GaudiMath::Interpolation::Akima :
128  T = gsl_interp_akima ; break ;
129  case GaudiMath::Interpolation::Akima_Periodic :
130  T = gsl_interp_akima_periodic ; break ;
131  default :
132  T = gsl_interp_cspline ; break ;
133  };
134 
135  m_spline = gsl_spline_alloc( T , m_dim ) ;
136 
137  gsl_spline_init( m_spline , m_x.get() , m_y.get() , m_dim ) ;
138 
139  m_accel = gsl_interp_accel_alloc() ;
140 
141  m_init = true ;
142 
143  }
144  // ========================================================================
145 
146  // ========================================================================
147  double SplineBase::eval ( const double x ) const
148  {
149  if ( !m_init ) { initialize() ; }
150  return gsl_spline_eval ( m_spline , x , m_accel );
151  }
152  // ========================================================================
153 
154  // ========================================================================
155  double SplineBase::deriv ( const double x ) const
156  {
157  if ( !m_init ) { initialize() ; }
158  return gsl_spline_eval_deriv ( m_spline , x , m_accel );
159  }
160  // ========================================================================
161 
162  // ========================================================================
163  double SplineBase::deriv2 ( const double x ) const
164  {
165  if ( !m_init ) { initialize() ; }
166  return gsl_spline_eval_deriv2 ( m_spline , x , m_accel );
167  }
168  // ========================================================================
169 
170  // ========================================================================
171  double SplineBase::integ ( const double a ,
172  const double b ) const
173  {
174  if ( !m_init ) { initialize() ; }
175  return gsl_spline_eval_integ ( m_spline , a , b , m_accel ) ;
176  }
177  // ========================================================================
178 
179 
180  // ========================================================================
181  FUNCTION_OBJECT_IMP( GSLSpline )
182  // ========================================================================
183 
184  // ========================================================================
205  // ========================================================================
206  GSLSpline::GSLSpline
207  ( const GSLSpline::Data1D& x ,
208  const GSLSpline::Data1D& y ,
209  const GaudiMath::Interpolation::Type type )
210  : AbsFunction ()
211  , m_spline ( x , y , type )
212  {}
213  // ========================================================================
214 
215  // ========================================================================
234  // ========================================================================
235  GSLSpline::GSLSpline
236  ( const GSLSpline::Data2D& data ,
237  const GaudiMath::Interpolation::Type type )
238  : AbsFunction ()
239  , m_spline ( data , type )
240  {}
241  // ========================================================================
242 
243  // ========================================================================
245  // ========================================================================
246  GSLSpline::GSLSpline
247  ( const SplineBase& right )
248  : AbsFunction ()
249  , m_spline ( right )
250  {}
251  // ========================================================================
252 
253  // ========================================================================
255  // ========================================================================
256  GSLSpline::GSLSpline
257  ( const GSLSpline& right )
258  : AbsFunction ()
259  , m_spline ( right )
260  {}
261  // ========================================================================
262 
263  // ========================================================================
265  // ========================================================================
266  GSLSpline::~GSLSpline() = default;
267  // ========================================================================
268 
269  // ========================================================================
270  double GSLSpline::operator() ( double x ) const
271  { return m_spline.eval( x ) ; }
272  // ========================================================================
273  double GSLSpline::operator() ( const Argument& x ) const
274  { return m_spline.eval( x[0] ) ; }
275  // ========================================================================
276 
277  // ========================================================================
279  // ========================================================================
280  Genfun::Derivative GSLSpline::partial( unsigned int i ) const
281  {
282  if ( i >= 1 )
283  {
284  const AbsFunction& aux = GaudiMath::Constant( 0.0 , 1 ) ;
285  return Genfun::FunctionNoop( &aux ) ;
286  }
287  const AbsFunction& aux = GSLSplineDeriv( *this ) ;
288  return Genfun::FunctionNoop( &aux ) ;
289  }
290  // ========================================================================
291 
292  // ========================================================================
293  FUNCTION_OBJECT_IMP( GSLSplineDeriv )
294  // ========================================================================
295 
296  // ========================================================================
317  // ========================================================================
318  GSLSplineDeriv::GSLSplineDeriv
319  ( const GSLSplineDeriv::Data1D& x ,
320  const GSLSplineDeriv::Data1D& y ,
321  const GaudiMath::Interpolation::Type type )
322  : AbsFunction ()
323  , m_spline ( x , y , type )
324  {}
325  // ========================================================================
326 
327  // ========================================================================
346  // ========================================================================
347  GSLSplineDeriv::GSLSplineDeriv
348  ( const GSLSplineDeriv::Data2D& data ,
349  const GaudiMath::Interpolation::Type type )
350  : AbsFunction ()
351  , m_spline ( data , type )
352  {}
353  // ========================================================================
354 
355  // ========================================================================
357  // ========================================================================
358  GSLSplineDeriv::GSLSplineDeriv
359  ( const SplineBase& right )
360  : AbsFunction ()
361  , m_spline ( right )
362  {}
363  // ========================================================================
364 
365  // ========================================================================
367  // ========================================================================
368  GSLSplineDeriv::GSLSplineDeriv
369  ( const GSLSplineDeriv& right )
370  : AbsFunction ()
371  , m_spline ( right )
372  {}
373  // ========================================================================
374 
375  // ========================================================================
377  // ========================================================================
378  GSLSplineDeriv::~GSLSplineDeriv(){}
379  // ========================================================================
380 
381  // ========================================================================
382  double GSLSplineDeriv::operator() ( double x ) const
383  { return m_spline.deriv ( x ) ; }
384  // ========================================================================
385  double GSLSplineDeriv::operator() ( const Argument& x ) const
386  { return m_spline.deriv ( x[0] ) ; }
387  // ========================================================================
388 
389  // ========================================================================
391  // ========================================================================
392  Genfun::Derivative GSLSplineDeriv::partial( unsigned int i ) const
393  {
394  if ( i >= 1 )
395  {
396  const AbsFunction& aux = GaudiMath::Constant( 0.0 , 1 ) ;
397  return Genfun::FunctionNoop( &aux ) ;
398  }
399  const AbsFunction& aux = GSLSplineDeriv2( *this ) ;
400  return Genfun::FunctionNoop( &aux ) ;
401  }
402  // ========================================================================
403 
404  // ========================================================================
405  FUNCTION_OBJECT_IMP( GSLSplineDeriv2 )
406  // ========================================================================
407 
408  // ========================================================================
429  // ========================================================================
430  GSLSplineDeriv2::GSLSplineDeriv2
431  ( const GSLSplineDeriv2::Data1D& x ,
432  const GSLSplineDeriv2::Data1D& y ,
433  const GaudiMath::Interpolation::Type type )
434  : AbsFunction ()
435  , m_spline ( x , y , type )
436  {}
437  // ========================================================================
438 
439  // ========================================================================
458  // ========================================================================
459  GSLSplineDeriv2::GSLSplineDeriv2
460  ( const GSLSplineDeriv2::Data2D& data ,
461  const GaudiMath::Interpolation::Type type )
462  : AbsFunction ()
463  , m_spline ( data , type )
464  {}
465  // ========================================================================
466 
467  // ========================================================================
469  // ========================================================================
470  GSLSplineDeriv2::GSLSplineDeriv2
471  ( const SplineBase& right )
472  : AbsFunction ()
473  , m_spline ( right )
474  {}
475  // ========================================================================
476 
477  // ========================================================================
479  // ========================================================================
480  GSLSplineDeriv2::GSLSplineDeriv2
481  ( const GSLSplineDeriv2& right )
482  : AbsFunction ()
483  , m_spline ( right )
484  {}
485  // ========================================================================
486 
487  // ========================================================================
489  // ========================================================================
490  GSLSplineDeriv2::~GSLSplineDeriv2(){}
491  // ========================================================================
492 
493  // ========================================================================
494  double GSLSplineDeriv2::operator() ( double x ) const
495  { return m_spline.deriv2 ( x ) ; }
496  // ========================================================================
497  double GSLSplineDeriv2::operator() ( const Argument& x ) const
498  { return m_spline.deriv2 ( x[0] ) ; }
499  // ========================================================================
500 
501  // ========================================================================
503  // ========================================================================
504  Genfun::Derivative GSLSplineDeriv2::partial( unsigned int i ) const
505  {
506  if ( i >= 1 )
507  {
508  const AbsFunction& aux = GaudiMath::Constant( 0.0 , 1 ) ;
509  return Genfun::FunctionNoop( &aux ) ;
510  }
511  const AbsFunction& aux = GaudiMath::Derivative( *this , i ) ;
512  return Genfun::FunctionNoop( &aux ) ;
513  }
514  // ========================================================================
515 
516 
517  // ========================================================================
518  FUNCTION_OBJECT_IMP( GSLSplineInteg )
519  // ========================================================================
520 
521  // ========================================================================
542  // ========================================================================
543  GSLSplineInteg::GSLSplineInteg
544  ( const GSLSplineInteg::Data1D& x ,
545  const GSLSplineInteg::Data1D& y ,
546  const GaudiMath::Interpolation::Type type ,
547  const double low )
548  : AbsFunction ( )
549  , m_spline ( x , y , type )
550  , m_low ( low )
551  {}
552  // ========================================================================
553 
554  // ========================================================================
573  // ========================================================================
574  GSLSplineInteg::GSLSplineInteg
575  ( const GSLSplineInteg::Data2D& data ,
576  const GaudiMath::Interpolation::Type type ,
577  const double low )
578  : AbsFunction ()
579  , m_spline ( data , type )
580  , m_low ( low )
581  {}
582  // ========================================================================
583 
584  // ========================================================================
586  // ========================================================================
587  GSLSplineInteg::GSLSplineInteg
588  ( const SplineBase& right ,
589  const double low )
590  : AbsFunction ()
591  , m_spline ( right )
592  , m_low ( low )
593  {}
594  // ========================================================================
595 
596  // ========================================================================
598  // ========================================================================
599  GSLSplineInteg::GSLSplineInteg
600  ( const GSLSplineInteg& right )
601  : AbsFunction ()
602  , m_spline ( right )
603  , m_low ( right.m_low )
604  {}
605  // ========================================================================
606 
607  // ========================================================================
609  // ========================================================================
610  GSLSplineInteg::~GSLSplineInteg(){}
611  // ========================================================================
612 
613  // ========================================================================
614  double GSLSplineInteg::operator() ( double x ) const
615  { return m_spline.integ ( m_low , x ) ; }
616  // ========================================================================
617  double GSLSplineInteg::operator() ( const Argument& x ) const
618  { return m_spline.integ ( m_low , x[0] ) ; }
619  // ========================================================================
620 
621  // ========================================================================
623  // ========================================================================
624  Genfun::Derivative GSLSplineInteg::partial( unsigned int i ) const
625  {
626  if ( i >= 1 )
627  {
628  const AbsFunction& aux = GaudiMath::Constant( 0.0 , 1 ) ;
629  return Genfun::FunctionNoop( &aux ) ;
630  }
631  const AbsFunction& aux = GSLSpline( *this ) ;
632  return Genfun::FunctionNoop( &aux ) ;
633  }
634  // ========================================================================
635 
636  }
637 }
Genfun::GaudiMathImplementation::SplineBase::integ
double integ(const double a, const double b) const
evaluate the integral on [a,b]
Definition: Splines.cpp:171
Genfun::GaudiMathImplementation::GSLSpline::GSLSpline
GSLSpline()=delete
default construtor is disabled ;
Genfun::GaudiMathImplementation::GSLSplineDeriv2::GSLSplineDeriv2
GSLSplineDeriv2()=delete
default construtor is disabled ;
Genfun::GaudiMathImplementation::GSLSplineDeriv::partial
Genfun::Derivative partial(unsigned int i) const override
Derivatives.
Definition: Splines.cpp:392
Genfun::GaudiMathImplementation::SplineBase::eval
double eval(const double x) const
evaluate the function
Definition: Splines.cpp:147
Genfun::GaudiMathImplementation::GSLSpline::operator()
virtual double operator()(double a) const
main methgod: evaluate teh function
Definition: Splines.cpp:270
Genfun::GaudiMathImplementation::SplineBase::m_y
std::unique_ptr< double[]> m_y
Definition: Splines.h:127
Genfun::GaudiMathImplementation::GSLSplineInteg::partial
Genfun::Derivative partial(unsigned int i) const override
Derivatives.
Definition: Splines.cpp:624
Genfun::GaudiMathImplementation::GSLSplineDeriv::operator()
double operator()(double a) const override
main method: evaluate the function
Definition: Splines.cpp:382
GaudiMath::Derivative
Genfun::GaudiMathImplementation::NumericalDerivative Derivative
Definition: GaudiMath.h:29
Genfun::GaudiMathImplementation::GSLSplineInteg::GSLSplineInteg
GSLSplineInteg()=delete
default construtor is disabled ;
Genfun::GaudiMathImplementation::SplineBase::m_type
GaudiMath::Interpolation::Type m_type
Definition: Splines.h:130
Genfun::GaudiMathImplementation::SplineBase::m_x
std::unique_ptr< double[]> m_x
Definition: Splines.h:126
Genfun::GaudiMathImplementation::SplineBase::deriv2
double deriv2(const double x) const
evaluate the second derivative
Definition: Splines.cpp:163
Genfun::GaudiMathImplementation::SplineBase::deriv
double deriv(const double x) const
evaluate the first derivative
Definition: Splines.cpp:155
Genfun::GaudiMathImplementation::GSLSplineDeriv::GSLSplineDeriv
GSLSplineDeriv()
default construtor is desabled ;
Genfun::GaudiMathImplementation::GSLSplineDeriv2::operator()
double operator()(double a) const override
main method: evaluate the function
Definition: Splines.cpp:494
GaudiMath
GaudiMath.h GaudiMath/GaudiMath.h.
Definition: Adapters.h:13
Genfun
CLHEP.
Definition: IEqSolver.h:13
Genfun::GaudiMathImplementation::GSLSpline::partial
virtual Genfun::Derivative partial(unsigned int i) const
Derivatives.
Definition: Splines.cpp:280
Genfun::GaudiMathImplementation::GSLSplineDeriv2::partial
Genfun::Derivative partial(unsigned int i) const override
Derivatives.
Definition: Splines.cpp:504
Genfun::GaudiMathImplementation::GSLSplineInteg::operator()
double operator()(double a) const override
main method: evaluate the function
Definition: Splines.cpp:614
ana.i
list i
Definition: ana.py:128
Tuples::Type
Type
the list of available types for ntuples
Definition: TupleObj.h:79
GaudiMath::Constant
Genfun::GaudiMathImplementation::Constant Constant
Definition: GaudiMath.h:27
Genfun::GaudiMathImplementation::SplineBase::Data2D
std::vector< std::pair< double, double > > Data2D
Definition: Splines.h:39
gaudirun.type
string type
Definition: gaudirun.py:151