All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
Splines.cpp
Go to the documentation of this file.
1 // $Id: Splines.cpp,v 1.2 2005/11/25 10:27:03 mato Exp $
2 // ============================================================================
3 // Include files
4 // ============================================================================
5 // GSL
6 // ============================================================================
7 #include "gsl/gsl_interp.h"
8 #include "gsl/gsl_spline.h"
9 // ============================================================================
10 // GaudiMath
11 // ============================================================================
12 #include "GaudiMath/Splines.h"
13 #include "GaudiMath/GaudiMath.h"
14 // ============================================================================
15 #include <cstring>
16 
17 
24 namespace Genfun
25 {
26  namespace GaudiMathImplementation
27  {
33  SplineBase::SplineBase
34  ( const SplineBase::Data1D& x ,
35  const SplineBase::Data1D& y ,
36  const GaudiMath::Interpolation::Type type )
37  : m_init ( false )
38  , m_dim ( x.size() )
39  , m_x ( new double[ x.size() ] )
40  , m_y ( new double[ y.size() ] )
41  , m_spline ( 0 )
42  , m_accel ( 0 )
43  , m_type ( type )
44  {
45 #ifdef WIN32
46 // Disable the warning
47 // C4996: 'std::copy': Function call with parameters that may be unsafe
48 // The parameters are checked
49 #pragma warning(push)
50 #pragma warning(disable:4996)
51 #endif
52  std::copy( x.begin() , x.end() , m_x ) ;
53  std::copy( y.begin() , y.end() , m_y ) ;
54 #ifdef WIN32
55 #pragma warning(pop)
56 #endif
57 
58  }
59  // ========================================================================
60 
61  // ========================================================================
66  // ========================================================================
67  SplineBase::SplineBase
68  ( const SplineBase::Data2D& data ,
69  const GaudiMath::Interpolation::Type type )
70  : m_init ( false )
71  , m_dim ( data.size() )
72  , m_x ( new double[ data.size() ] )
73  , m_y ( new double[ data.size() ] )
74  , m_spline ( 0 )
75  , m_accel ( 0 )
76  , m_type ( type )
77  {
78  double* _x = m_x ;
79  double* _y = m_y ;
80  for( Data2D::const_iterator it =
81  data.begin() ; data.end() != it ; ++it )
82  {
83  *_x = it -> first ; ++_x ;
84  *_y = it -> second ; ++_y ;
85  }
86  }
87  // ========================================================================
88 
89  // ========================================================================
91  // ========================================================================
92  SplineBase::SplineBase ( const SplineBase& right )
93  : m_init ( false )
94  , m_dim ( right.m_dim )
95  , m_x ( new double[ right.m_dim ] )
96  , m_y ( new double[ right.m_dim ] )
97  , m_spline ( 0 )
98  , m_accel ( 0 )
99  , m_type ( right.m_type )
100  {
101  const std::size_t size = sizeof(double) * m_dim;
102  std::memcpy(m_x, right.m_x, size);
103  std::memcpy(m_y, right.m_y, size);
104  }
105  // ========================================================================
106 
107  // ========================================================================
109  // ========================================================================
110  SplineBase::~SplineBase()
111  {
112  if ( 0 != m_spline ) { gsl_spline_free ( m_spline ) ; }
113  if ( 0 != m_accel ) { gsl_interp_accel_free ( m_accel ) ; }
114 
115  if ( 0 != m_x ) { delete[] m_x ; }
116  if ( 0 != m_y ) { delete[] m_y ; }
117  }
118  // ========================================================================
119 
120  // ========================================================================
121  void SplineBase::initialize() const
122  {
123  if ( m_init ) { return ; } // RETURN
124 
125  const gsl_interp_type* T = 0 ;
126 
127  switch ( m_type )
128  {
129  case GaudiMath::Interpolation::Linear :
130  T = gsl_interp_linear ; break ;
131  case GaudiMath::Interpolation::Polynomial :
132  T = gsl_interp_polynomial ; break ;
133  case GaudiMath::Interpolation::Cspline :
134  T = gsl_interp_cspline ; break ;
135  case GaudiMath::Interpolation::Cspline_Periodic :
136  T = gsl_interp_cspline_periodic ; break ;
137  case GaudiMath::Interpolation::Akima :
138  T = gsl_interp_akima ; break ;
139  case GaudiMath::Interpolation::Akima_Periodic :
140  T = gsl_interp_akima_periodic ; break ;
141  default :
142  T = gsl_interp_cspline ; break ;
143  };
144 
145  m_spline = gsl_spline_alloc( T , m_dim ) ;
146 
147  gsl_spline_init( m_spline , m_x , m_y , m_dim ) ;
148 
149  m_accel = gsl_interp_accel_alloc() ;
150 
151  m_init = true ;
152 
153  }
154  // ========================================================================
155 
156  // ========================================================================
157  double SplineBase::eval ( const double x ) const
158  {
159  if ( !m_init ) { initialize() ; }
160  return gsl_spline_eval ( m_spline , x , m_accel );
161  }
162  // ========================================================================
163 
164  // ========================================================================
165  double SplineBase::deriv ( const double x ) const
166  {
167  if ( !m_init ) { initialize() ; }
168  return gsl_spline_eval_deriv ( m_spline , x , m_accel );
169  }
170  // ========================================================================
171 
172  // ========================================================================
173  double SplineBase::deriv2 ( const double x ) const
174  {
175  if ( !m_init ) { initialize() ; }
176  return gsl_spline_eval_deriv2 ( m_spline , x , m_accel );
177  }
178  // ========================================================================
179 
180  // ========================================================================
181  double SplineBase::integ ( const double a ,
182  const double b ) const
183  {
184  if ( !m_init ) { initialize() ; }
185  return gsl_spline_eval_integ ( m_spline , a , b , m_accel ) ;
186  }
187  // ========================================================================
188 
189 
190  // ========================================================================
191  FUNCTION_OBJECT_IMP( GSLSpline )
192  // ========================================================================
193 
194  // ========================================================================
215  // ========================================================================
216  GSLSpline::GSLSpline
217  ( const GSLSpline::Data1D& x ,
218  const GSLSpline::Data1D& y ,
219  const GaudiMath::Interpolation::Type type )
220  : AbsFunction ()
221  , m_spline ( x , y , type )
222  {}
223  // ========================================================================
224 
225  // ========================================================================
244  // ========================================================================
245  GSLSpline::GSLSpline
246  ( const GSLSpline::Data2D& data ,
247  const GaudiMath::Interpolation::Type type )
248  : AbsFunction ()
249  , m_spline ( data , type )
250  {}
251  // ========================================================================
252 
253  // ========================================================================
255  // ========================================================================
256  GSLSpline::GSLSpline
257  ( const SplineBase& right )
258  : AbsFunction ()
259  , m_spline ( right )
260  {}
261  // ========================================================================
262 
263  // ========================================================================
265  // ========================================================================
266  GSLSpline::GSLSpline
267  ( const GSLSpline& right )
268  : AbsFunction ()
269  , m_spline ( right )
270  {}
271  // ========================================================================
272 
273  // ========================================================================
275  // ========================================================================
276  GSLSpline::~GSLSpline(){}
277  // ========================================================================
278 
279  // ========================================================================
280  double GSLSpline::operator() ( double x ) const
281  { return m_spline.eval( x ) ; }
282  // ========================================================================
283  double GSLSpline::operator() ( const Argument& x ) const
284  { return m_spline.eval( x[0] ) ; }
285  // ========================================================================
286 
287  // ========================================================================
289  // ========================================================================
290  Genfun::Derivative GSLSpline::partial( unsigned int i ) const
291  {
292  if ( i >= 1 )
293  {
294  const AbsFunction& aux = GaudiMath::Constant( 0.0 , 1 ) ;
295  return Genfun::FunctionNoop( &aux ) ;
296  }
297  const AbsFunction& aux = GSLSplineDeriv( *this ) ;
298  return Genfun::FunctionNoop( &aux ) ;
299  }
300  // ========================================================================
301 
302  // ========================================================================
303  FUNCTION_OBJECT_IMP( GSLSplineDeriv )
304  // ========================================================================
305 
306  // ========================================================================
327  // ========================================================================
328  GSLSplineDeriv::GSLSplineDeriv
329  ( const GSLSplineDeriv::Data1D& x ,
330  const GSLSplineDeriv::Data1D& y ,
331  const GaudiMath::Interpolation::Type type )
332  : AbsFunction ()
333  , m_spline ( x , y , type )
334  {}
335  // ========================================================================
336 
337  // ========================================================================
356  // ========================================================================
357  GSLSplineDeriv::GSLSplineDeriv
358  ( const GSLSplineDeriv::Data2D& data ,
359  const GaudiMath::Interpolation::Type type )
360  : AbsFunction ()
361  , m_spline ( data , type )
362  {}
363  // ========================================================================
364 
365  // ========================================================================
367  // ========================================================================
368  GSLSplineDeriv::GSLSplineDeriv
369  ( const SplineBase& right )
370  : AbsFunction ()
371  , m_spline ( right )
372  {}
373  // ========================================================================
374 
375  // ========================================================================
377  // ========================================================================
378  GSLSplineDeriv::GSLSplineDeriv
379  ( const GSLSplineDeriv& right )
380  : AbsFunction ()
381  , m_spline ( right )
382  {}
383  // ========================================================================
384 
385  // ========================================================================
387  // ========================================================================
388  GSLSplineDeriv::~GSLSplineDeriv(){}
389  // ========================================================================
390 
391  // ========================================================================
392  double GSLSplineDeriv::operator() ( double x ) const
393  { return m_spline.deriv ( x ) ; }
394  // ========================================================================
395  double GSLSplineDeriv::operator() ( const Argument& x ) const
396  { return m_spline.deriv ( x[0] ) ; }
397  // ========================================================================
398 
399  // ========================================================================
401  // ========================================================================
402  Genfun::Derivative GSLSplineDeriv::partial( unsigned int i ) const
403  {
404  if ( i >= 1 )
405  {
406  const AbsFunction& aux = GaudiMath::Constant( 0.0 , 1 ) ;
407  return Genfun::FunctionNoop( &aux ) ;
408  }
409  const AbsFunction& aux = GSLSplineDeriv2( *this ) ;
410  return Genfun::FunctionNoop( &aux ) ;
411  }
412  // ========================================================================
413 
414  // ========================================================================
415  FUNCTION_OBJECT_IMP( GSLSplineDeriv2 )
416  // ========================================================================
417 
418  // ========================================================================
439  // ========================================================================
440  GSLSplineDeriv2::GSLSplineDeriv2
441  ( const GSLSplineDeriv2::Data1D& x ,
442  const GSLSplineDeriv2::Data1D& y ,
443  const GaudiMath::Interpolation::Type type )
444  : AbsFunction ()
445  , m_spline ( x , y , type )
446  {}
447  // ========================================================================
448 
449  // ========================================================================
468  // ========================================================================
469  GSLSplineDeriv2::GSLSplineDeriv2
470  ( const GSLSplineDeriv2::Data2D& data ,
471  const GaudiMath::Interpolation::Type type )
472  : AbsFunction ()
473  , m_spline ( data , type )
474  {}
475  // ========================================================================
476 
477  // ========================================================================
479  // ========================================================================
480  GSLSplineDeriv2::GSLSplineDeriv2
481  ( const SplineBase& right )
482  : AbsFunction ()
483  , m_spline ( right )
484  {}
485  // ========================================================================
486 
487  // ========================================================================
489  // ========================================================================
490  GSLSplineDeriv2::GSLSplineDeriv2
491  ( const GSLSplineDeriv2& right )
492  : AbsFunction ()
493  , m_spline ( right )
494  {}
495  // ========================================================================
496 
497  // ========================================================================
499  // ========================================================================
500  GSLSplineDeriv2::~GSLSplineDeriv2(){}
501  // ========================================================================
502 
503  // ========================================================================
504  double GSLSplineDeriv2::operator() ( double x ) const
505  { return m_spline.deriv2 ( x ) ; }
506  // ========================================================================
507  double GSLSplineDeriv2::operator() ( const Argument& x ) const
508  { return m_spline.deriv2 ( x[0] ) ; }
509  // ========================================================================
510 
511  // ========================================================================
513  // ========================================================================
514  Genfun::Derivative GSLSplineDeriv2::partial( unsigned int i ) const
515  {
516  if ( i >= 1 )
517  {
518  const AbsFunction& aux = GaudiMath::Constant( 0.0 , 1 ) ;
519  return Genfun::FunctionNoop( &aux ) ;
520  }
521  const AbsFunction& aux = GaudiMath::Derivative( *this , i ) ;
522  return Genfun::FunctionNoop( &aux ) ;
523  }
524  // ========================================================================
525 
526 
527  // ========================================================================
528  FUNCTION_OBJECT_IMP( GSLSplineInteg )
529  // ========================================================================
530 
531  // ========================================================================
552  // ========================================================================
553  GSLSplineInteg::GSLSplineInteg
554  ( const GSLSplineInteg::Data1D& x ,
555  const GSLSplineInteg::Data1D& y ,
556  const GaudiMath::Interpolation::Type type ,
557  const double low )
558  : AbsFunction ( )
559  , m_spline ( x , y , type )
560  , m_low ( low )
561  {}
562  // ========================================================================
563 
564  // ========================================================================
583  // ========================================================================
584  GSLSplineInteg::GSLSplineInteg
585  ( const GSLSplineInteg::Data2D& data ,
586  const GaudiMath::Interpolation::Type type ,
587  const double low )
588  : AbsFunction ()
589  , m_spline ( data , type )
590  , m_low ( low )
591  {}
592  // ========================================================================
593 
594  // ========================================================================
596  // ========================================================================
597  GSLSplineInteg::GSLSplineInteg
598  ( const SplineBase& right ,
599  const double low )
600  : AbsFunction ()
601  , m_spline ( right )
602  , m_low ( low )
603  {}
604  // ========================================================================
605 
606  // ========================================================================
608  // ========================================================================
609  GSLSplineInteg::GSLSplineInteg
610  ( const GSLSplineInteg& right )
611  : AbsFunction ()
612  , m_spline ( right )
613  , m_low ( right.m_low )
614  {}
615  // ========================================================================
616 
617  // ========================================================================
619  // ========================================================================
620  GSLSplineInteg::~GSLSplineInteg(){}
621  // ========================================================================
622 
623  // ========================================================================
624  double GSLSplineInteg::operator() ( double x ) const
625  { return m_spline.integ ( m_low , x ) ; }
626  // ========================================================================
627  double GSLSplineInteg::operator() ( const Argument& x ) const
628  { return m_spline.integ ( m_low , x[0] ) ; }
629  // ========================================================================
630 
631  // ========================================================================
633  // ========================================================================
634  Genfun::Derivative GSLSplineInteg::partial( unsigned int i ) const
635  {
636  if ( i >= 1 )
637  {
638  const AbsFunction& aux = GaudiMath::Constant( 0.0 , 1 ) ;
639  return Genfun::FunctionNoop( &aux ) ;
640  }
641  const AbsFunction& aux = GSLSpline( *this ) ;
642  return Genfun::FunctionNoop( &aux ) ;
643  }
644  // ========================================================================
645 
646  }
647 }
Genfun::GaudiMathImplementation::SplineBase::integ
double integ(const double a, const double b) const
evaluate the integral on [a,b]
Definition: Splines.cpp:181
Genfun::GaudiMathImplementation::GSLSplineDeriv2::GSLSplineDeriv2
GSLSplineDeriv2()
default construtor is desabled ;
Genfun::GaudiMathImplementation::GSLSpline::GSLSpline
GSLSpline()
default construtor is desabled ;
Genfun::GaudiMathImplementation::GSLSplineInteg::operator()
virtual double operator()(double a) const
main method: evaluate the function
Definition: Splines.cpp:624
Genfun::GaudiMathImplementation::SplineBase::eval
double eval(const double x) const
evaluate the function
Definition: Splines.cpp:157
Genfun::GaudiMathImplementation::GSLSpline::operator()
virtual double operator()(double a) const
main methgod: evaluate teh function
Definition: Splines.cpp:280
Splines.h
GaudiMath::Derivative
Genfun::GaudiMathImplementation::NumericalDerivative Derivative
Definition: GaudiMath.h:31
Genfun::GaudiMathImplementation::GSLSplineDeriv::partial
virtual Genfun::Derivative partial(unsigned int i) const
Derivatives.
Definition: Splines.cpp:402
Genfun::GaudiMathImplementation::GSLSplineInteg::partial
virtual Genfun::Derivative partial(unsigned int i) const
Derivatives.
Definition: Splines.cpp:634
Genfun::GaudiMathImplementation::GSLSplineDeriv::operator()
virtual double operator()(double a) const
main method: evaluate the function
Definition: Splines.cpp:392
gaudirun.type
string type
Definition: gaudirun.py:126
Genfun::GaudiMathImplementation::SplineBase::m_type
GaudiMath::Interpolation::Type m_type
Definition: Splines.h:137
Genfun::GaudiMathImplementation::GSLSplineDeriv2::partial
virtual Genfun::Derivative partial(unsigned int i) const
Derivatives.
Definition: Splines.cpp:514
Genfun::GaudiMathImplementation::SplineBase::deriv2
double deriv2(const double x) const
evaluate the second derivative
Definition: Splines.cpp:173
Genfun::GaudiMathImplementation::SplineBase::deriv
double deriv(const double x) const
evaluate the first derivative
Definition: Splines.cpp:165
Genfun::GaudiMathImplementation::GSLSplineDeriv::GSLSplineDeriv
GSLSplineDeriv()
default construtor is desabled ;
Genfun::GaudiMathImplementation::GSLSplineInteg::GSLSplineInteg
GSLSplineInteg()
default construtor is desabled ;
Genfun::GaudiMathImplementation::GSLSplineDeriv2::operator()
virtual double operator()(double a) const
main method: evaluate the function
Definition: Splines.cpp:504
Genfun::GaudiMathImplementation::GSLSpline::partial
virtual Genfun::Derivative partial(unsigned int i) const
Derivatives.
Definition: Splines.cpp:290
ana.i
list i
Definition: ana.py:128
Tuples::Type
Type
the list of available types for ntuples
Definition: TupleObj.h:63
GaudiMath::Constant
Genfun::GaudiMathImplementation::Constant Constant
Definition: GaudiMath.h:29
GaudiMath.h
Genfun::GaudiMathImplementation::SplineBase::Data2D
std::vector< std::pair< double, double > > Data2D
Definition: Splines.h:40