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 ( 0 )
41  , m_accel ( 0 )
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_spline ( 0 )
74  , m_accel ( 0 )
75  , m_type ( type )
76  {
77  double* _x = m_x.get() ;
78  double* _y = m_y.get() ;
79  for( const auto& i : data )
80  {
81  *_x = i.first ; ++_x ;
82  *_y = i.second ; ++_y ;
83  }
84  }
85  // ========================================================================
86 
87  // ========================================================================
89  // ========================================================================
90  SplineBase::SplineBase ( const SplineBase& right )
91  : m_init ( false )
92  , m_dim ( right.m_dim )
93  , m_x ( new double[ right.m_dim ] )
94  , m_y ( new double[ right.m_dim ] )
95  , m_spline ( 0 )
96  , m_accel ( 0 )
97  , m_type ( right.m_type )
98  {
99  std::copy_n(right.m_x.get(),m_dim,m_x.get());
100  std::copy_n(right.m_y.get(),m_dim,m_y.get());
101  }
102  // ========================================================================
103 
104  // ========================================================================
106  // ========================================================================
107  SplineBase::~SplineBase()
108  {
109  if ( 0 != m_spline ) { gsl_spline_free ( m_spline ) ; }
110  if ( 0 != m_accel ) { gsl_interp_accel_free ( m_accel ) ; }
111  }
112  // ========================================================================
113 
114  // ========================================================================
115  void SplineBase::initialize() const
116  {
117  if ( m_init ) { return ; } // RETURN
118 
119  const gsl_interp_type* T = 0 ;
120 
121  switch ( m_type )
122  {
123  case GaudiMath::Interpolation::Linear :
124  T = gsl_interp_linear ; break ;
125  case GaudiMath::Interpolation::Polynomial :
126  T = gsl_interp_polynomial ; break ;
127  case GaudiMath::Interpolation::Cspline :
128  T = gsl_interp_cspline ; break ;
129  case GaudiMath::Interpolation::Cspline_Periodic :
130  T = gsl_interp_cspline_periodic ; break ;
131  case GaudiMath::Interpolation::Akima :
132  T = gsl_interp_akima ; break ;
133  case GaudiMath::Interpolation::Akima_Periodic :
134  T = gsl_interp_akima_periodic ; break ;
135  default :
136  T = gsl_interp_cspline ; break ;
137  };
138 
139  m_spline = gsl_spline_alloc( T , m_dim ) ;
140 
141  gsl_spline_init( m_spline , m_x.get() , m_y.get() , m_dim ) ;
142 
143  m_accel = gsl_interp_accel_alloc() ;
144 
145  m_init = true ;
146 
147  }
148  // ========================================================================
149 
150  // ========================================================================
151  double SplineBase::eval ( const double x ) const
152  {
153  if ( !m_init ) { initialize() ; }
154  return gsl_spline_eval ( m_spline , x , m_accel );
155  }
156  // ========================================================================
157 
158  // ========================================================================
159  double SplineBase::deriv ( const double x ) const
160  {
161  if ( !m_init ) { initialize() ; }
162  return gsl_spline_eval_deriv ( m_spline , x , m_accel );
163  }
164  // ========================================================================
165 
166  // ========================================================================
167  double SplineBase::deriv2 ( const double x ) const
168  {
169  if ( !m_init ) { initialize() ; }
170  return gsl_spline_eval_deriv2 ( m_spline , x , m_accel );
171  }
172  // ========================================================================
173 
174  // ========================================================================
175  double SplineBase::integ ( const double a ,
176  const double b ) const
177  {
178  if ( !m_init ) { initialize() ; }
179  return gsl_spline_eval_integ ( m_spline , a , b , m_accel ) ;
180  }
181  // ========================================================================
182 
183 
184  // ========================================================================
185  FUNCTION_OBJECT_IMP( GSLSpline )
186  // ========================================================================
187 
188  // ========================================================================
209  // ========================================================================
210  GSLSpline::GSLSpline
211  ( const GSLSpline::Data1D& x ,
212  const GSLSpline::Data1D& y ,
213  const GaudiMath::Interpolation::Type type )
214  : AbsFunction ()
215  , m_spline ( x , y , type )
216  {}
217  // ========================================================================
218 
219  // ========================================================================
238  // ========================================================================
239  GSLSpline::GSLSpline
240  ( const GSLSpline::Data2D& data ,
241  const GaudiMath::Interpolation::Type type )
242  : AbsFunction ()
243  , m_spline ( data , type )
244  {}
245  // ========================================================================
246 
247  // ========================================================================
249  // ========================================================================
250  GSLSpline::GSLSpline
251  ( const SplineBase& right )
252  : AbsFunction ()
253  , m_spline ( right )
254  {}
255  // ========================================================================
256 
257  // ========================================================================
259  // ========================================================================
260  GSLSpline::GSLSpline
261  ( const GSLSpline& right )
262  : AbsFunction ()
263  , m_spline ( right )
264  {}
265  // ========================================================================
266 
267  // ========================================================================
269  // ========================================================================
270  GSLSpline::~GSLSpline(){}
271  // ========================================================================
272 
273  // ========================================================================
274  double GSLSpline::operator() ( double x ) const
275  { return m_spline.eval( x ) ; }
276  // ========================================================================
277  double GSLSpline::operator() ( const Argument& x ) const
278  { return m_spline.eval( x[0] ) ; }
279  // ========================================================================
280 
281  // ========================================================================
283  // ========================================================================
284  Genfun::Derivative GSLSpline::partial( unsigned int i ) const
285  {
286  if ( i >= 1 )
287  {
288  const AbsFunction& aux = GaudiMath::Constant( 0.0 , 1 ) ;
289  return Genfun::FunctionNoop( &aux ) ;
290  }
291  const AbsFunction& aux = GSLSplineDeriv( *this ) ;
292  return Genfun::FunctionNoop( &aux ) ;
293  }
294  // ========================================================================
295 
296  // ========================================================================
297  FUNCTION_OBJECT_IMP( GSLSplineDeriv )
298  // ========================================================================
299 
300  // ========================================================================
321  // ========================================================================
322  GSLSplineDeriv::GSLSplineDeriv
323  ( const GSLSplineDeriv::Data1D& x ,
324  const GSLSplineDeriv::Data1D& y ,
325  const GaudiMath::Interpolation::Type type )
326  : AbsFunction ()
327  , m_spline ( x , y , type )
328  {}
329  // ========================================================================
330 
331  // ========================================================================
350  // ========================================================================
351  GSLSplineDeriv::GSLSplineDeriv
352  ( const GSLSplineDeriv::Data2D& data ,
353  const GaudiMath::Interpolation::Type type )
354  : AbsFunction ()
355  , m_spline ( data , type )
356  {}
357  // ========================================================================
358 
359  // ========================================================================
361  // ========================================================================
362  GSLSplineDeriv::GSLSplineDeriv
363  ( const SplineBase& right )
364  : AbsFunction ()
365  , m_spline ( right )
366  {}
367  // ========================================================================
368 
369  // ========================================================================
371  // ========================================================================
372  GSLSplineDeriv::GSLSplineDeriv
373  ( const GSLSplineDeriv& right )
374  : AbsFunction ()
375  , m_spline ( right )
376  {}
377  // ========================================================================
378 
379  // ========================================================================
381  // ========================================================================
382  GSLSplineDeriv::~GSLSplineDeriv(){}
383  // ========================================================================
384 
385  // ========================================================================
386  double GSLSplineDeriv::operator() ( double x ) const
387  { return m_spline.deriv ( x ) ; }
388  // ========================================================================
389  double GSLSplineDeriv::operator() ( const Argument& x ) const
390  { return m_spline.deriv ( x[0] ) ; }
391  // ========================================================================
392 
393  // ========================================================================
395  // ========================================================================
396  Genfun::Derivative GSLSplineDeriv::partial( unsigned int i ) const
397  {
398  if ( i >= 1 )
399  {
400  const AbsFunction& aux = GaudiMath::Constant( 0.0 , 1 ) ;
401  return Genfun::FunctionNoop( &aux ) ;
402  }
403  const AbsFunction& aux = GSLSplineDeriv2( *this ) ;
404  return Genfun::FunctionNoop( &aux ) ;
405  }
406  // ========================================================================
407 
408  // ========================================================================
409  FUNCTION_OBJECT_IMP( GSLSplineDeriv2 )
410  // ========================================================================
411 
412  // ========================================================================
433  // ========================================================================
434  GSLSplineDeriv2::GSLSplineDeriv2
435  ( const GSLSplineDeriv2::Data1D& x ,
436  const GSLSplineDeriv2::Data1D& y ,
437  const GaudiMath::Interpolation::Type type )
438  : AbsFunction ()
439  , m_spline ( x , y , type )
440  {}
441  // ========================================================================
442 
443  // ========================================================================
462  // ========================================================================
463  GSLSplineDeriv2::GSLSplineDeriv2
464  ( const GSLSplineDeriv2::Data2D& data ,
465  const GaudiMath::Interpolation::Type type )
466  : AbsFunction ()
467  , m_spline ( data , type )
468  {}
469  // ========================================================================
470 
471  // ========================================================================
473  // ========================================================================
474  GSLSplineDeriv2::GSLSplineDeriv2
475  ( const SplineBase& right )
476  : AbsFunction ()
477  , m_spline ( right )
478  {}
479  // ========================================================================
480 
481  // ========================================================================
483  // ========================================================================
484  GSLSplineDeriv2::GSLSplineDeriv2
485  ( const GSLSplineDeriv2& right )
486  : AbsFunction ()
487  , m_spline ( right )
488  {}
489  // ========================================================================
490 
491  // ========================================================================
493  // ========================================================================
494  GSLSplineDeriv2::~GSLSplineDeriv2(){}
495  // ========================================================================
496 
497  // ========================================================================
498  double GSLSplineDeriv2::operator() ( double x ) const
499  { return m_spline.deriv2 ( x ) ; }
500  // ========================================================================
501  double GSLSplineDeriv2::operator() ( const Argument& x ) const
502  { return m_spline.deriv2 ( x[0] ) ; }
503  // ========================================================================
504 
505  // ========================================================================
507  // ========================================================================
508  Genfun::Derivative GSLSplineDeriv2::partial( unsigned int i ) const
509  {
510  if ( i >= 1 )
511  {
512  const AbsFunction& aux = GaudiMath::Constant( 0.0 , 1 ) ;
513  return Genfun::FunctionNoop( &aux ) ;
514  }
515  const AbsFunction& aux = GaudiMath::Derivative( *this , i ) ;
516  return Genfun::FunctionNoop( &aux ) ;
517  }
518  // ========================================================================
519 
520 
521  // ========================================================================
522  FUNCTION_OBJECT_IMP( GSLSplineInteg )
523  // ========================================================================
524 
525  // ========================================================================
546  // ========================================================================
547  GSLSplineInteg::GSLSplineInteg
548  ( const GSLSplineInteg::Data1D& x ,
549  const GSLSplineInteg::Data1D& y ,
550  const GaudiMath::Interpolation::Type type ,
551  const double low )
552  : AbsFunction ( )
553  , m_spline ( x , y , type )
554  , m_low ( low )
555  {}
556  // ========================================================================
557 
558  // ========================================================================
577  // ========================================================================
578  GSLSplineInteg::GSLSplineInteg
579  ( const GSLSplineInteg::Data2D& data ,
580  const GaudiMath::Interpolation::Type type ,
581  const double low )
582  : AbsFunction ()
583  , m_spline ( data , type )
584  , m_low ( low )
585  {}
586  // ========================================================================
587 
588  // ========================================================================
590  // ========================================================================
591  GSLSplineInteg::GSLSplineInteg
592  ( const SplineBase& right ,
593  const double low )
594  : AbsFunction ()
595  , m_spline ( right )
596  , m_low ( low )
597  {}
598  // ========================================================================
599 
600  // ========================================================================
602  // ========================================================================
603  GSLSplineInteg::GSLSplineInteg
604  ( const GSLSplineInteg& right )
605  : AbsFunction ()
606  , m_spline ( right )
607  , m_low ( right.m_low )
608  {}
609  // ========================================================================
610 
611  // ========================================================================
613  // ========================================================================
614  GSLSplineInteg::~GSLSplineInteg(){}
615  // ========================================================================
616 
617  // ========================================================================
618  double GSLSplineInteg::operator() ( double x ) const
619  { return m_spline.integ ( m_low , x ) ; }
620  // ========================================================================
621  double GSLSplineInteg::operator() ( const Argument& x ) const
622  { return m_spline.integ ( m_low , x[0] ) ; }
623  // ========================================================================
624 
625  // ========================================================================
627  // ========================================================================
628  Genfun::Derivative GSLSplineInteg::partial( unsigned int i ) const
629  {
630  if ( i >= 1 )
631  {
632  const AbsFunction& aux = GaudiMath::Constant( 0.0 , 1 ) ;
633  return Genfun::FunctionNoop( &aux ) ;
634  }
635  const AbsFunction& aux = GSLSpline( *this ) ;
636  return Genfun::FunctionNoop( &aux ) ;
637  }
638  // ========================================================================
639 
640  }
641 }
Genfun::GaudiMathImplementation::SplineBase::integ
double integ(const double a, const double b) const
evaluate the integral on [a,b]
Definition: Splines.cpp:175
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:396
Genfun::GaudiMathImplementation::SplineBase::eval
double eval(const double x) const
evaluate the function
Definition: Splines.cpp:151
Genfun::GaudiMathImplementation::GSLSpline::operator()
virtual double operator()(double a) const
main methgod: evaluate teh function
Definition: Splines.cpp:274
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:628
Genfun::GaudiMathImplementation::GSLSplineDeriv::operator()
double operator()(double a) const override
main method: evaluate the function
Definition: Splines.cpp:386
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:167
Genfun::GaudiMathImplementation::SplineBase::deriv
double deriv(const double x) const
evaluate the first derivative
Definition: Splines.cpp:159
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:498
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:284
Genfun::GaudiMathImplementation::GSLSplineDeriv2::partial
Genfun::Derivative partial(unsigned int i) const override
Derivatives.
Definition: Splines.cpp:508
Genfun::GaudiMathImplementation::GSLSplineInteg::operator()
double operator()(double a) const override
main method: evaluate the function
Definition: Splines.cpp:618
ana.i
list i
Definition: ana.py:128
Tuples::Type
Type
the list of available types for ntuples
Definition: TupleObj.h:73
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