Gaudi Framework, version v23r5
Home   Generated: Wed Nov 28 2012
 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  std::memcpy(m_x, right.m_x, m_dim);
102  std::memcpy(m_y, right.m_y, m_dim);
103  }
104  // ========================================================================
105 
106  // ========================================================================
108  // ========================================================================
109  SplineBase::~SplineBase()
110  {
111  if ( 0 != m_spline ) { gsl_spline_free ( m_spline ) ; }
112  if ( 0 != m_accel ) { gsl_interp_accel_free ( m_accel ) ; }
113 
114  if ( 0 != m_x ) { delete[] m_x ; }
115  if ( 0 != m_y ) { delete[] m_y ; }
116  }
117  // ========================================================================
118 
119  // ========================================================================
120  void SplineBase::initialize() const
121  {
122  if ( m_init ) { return ; } // RETURN
123 
124  const gsl_interp_type* T = 0 ;
125 
126  switch ( m_type )
127  {
128  case GaudiMath::Interpolation::Linear :
129  T = gsl_interp_linear ; break ;
130  case GaudiMath::Interpolation::Polynomial :
131  T = gsl_interp_polynomial ; break ;
132  case GaudiMath::Interpolation::Cspline :
133  T = gsl_interp_cspline ; break ;
134  case GaudiMath::Interpolation::Cspline_Periodic :
135  T = gsl_interp_cspline_periodic ; break ;
136  case GaudiMath::Interpolation::Akima :
137  T = gsl_interp_akima ; break ;
138  case GaudiMath::Interpolation::Akima_Periodic :
139  T = gsl_interp_akima_periodic ; break ;
140  default :
141  T = gsl_interp_cspline ; break ;
142  };
143 
144  m_spline = gsl_spline_alloc( T , m_dim ) ;
145 
146  gsl_spline_init( m_spline , m_x , m_y , m_dim ) ;
147 
148  m_accel = gsl_interp_accel_alloc() ;
149 
150  m_init = true ;
151 
152  }
153  // ========================================================================
154 
155  // ========================================================================
156  double SplineBase::eval ( const double x ) const
157  {
158  if ( !m_init ) { initialize() ; }
159  return gsl_spline_eval ( m_spline , x , m_accel );
160  }
161  // ========================================================================
162 
163  // ========================================================================
164  double SplineBase::deriv ( const double x ) const
165  {
166  if ( !m_init ) { initialize() ; }
167  return gsl_spline_eval_deriv ( m_spline , x , m_accel );
168  }
169  // ========================================================================
170 
171  // ========================================================================
172  double SplineBase::deriv2 ( const double x ) const
173  {
174  if ( !m_init ) { initialize() ; }
175  return gsl_spline_eval_deriv2 ( m_spline , x , m_accel );
176  }
177  // ========================================================================
178 
179  // ========================================================================
180  double SplineBase::integ ( const double a ,
181  const double b ) const
182  {
183  if ( !m_init ) { initialize() ; }
184  return gsl_spline_eval_integ ( m_spline , a , b , m_accel ) ;
185  }
186  // ========================================================================
187 
188 
189  // ========================================================================
190  FUNCTION_OBJECT_IMP( GSLSpline )
191  // ========================================================================
192 
193  // ========================================================================
214  // ========================================================================
215  GSLSpline::GSLSpline
216  ( const GSLSpline::Data1D& x ,
217  const GSLSpline::Data1D& y ,
218  const GaudiMath::Interpolation::Type type )
219  : AbsFunction ()
220  , m_spline ( x , y , type )
221  {}
222  // ========================================================================
223 
224  // ========================================================================
243  // ========================================================================
244  GSLSpline::GSLSpline
245  ( const GSLSpline::Data2D& data ,
246  const GaudiMath::Interpolation::Type type )
247  : AbsFunction ()
248  , m_spline ( data , type )
249  {}
250  // ========================================================================
251 
252  // ========================================================================
254  // ========================================================================
255  GSLSpline::GSLSpline
256  ( const SplineBase& right )
257  : AbsFunction ()
258  , m_spline ( right )
259  {}
260  // ========================================================================
261 
262  // ========================================================================
264  // ========================================================================
265  GSLSpline::GSLSpline
266  ( const GSLSpline& right )
267  : AbsFunction ()
268  , m_spline ( right )
269  {}
270  // ========================================================================
271 
272  // ========================================================================
274  // ========================================================================
275  GSLSpline::~GSLSpline(){}
276  // ========================================================================
277 
278  // ========================================================================
279  double GSLSpline::operator() ( double x ) const
280  { return m_spline.eval( x ) ; }
281  // ========================================================================
282  double GSLSpline::operator() ( const Argument& x ) const
283  { return m_spline.eval( x[0] ) ; }
284  // ========================================================================
285 
286  // ========================================================================
288  // ========================================================================
289  Genfun::Derivative GSLSpline::partial( unsigned int i ) const
290  {
291  if ( i >= 1 )
292  {
293  const AbsFunction& aux = GaudiMath::Constant( 0.0 , 1 ) ;
294  return Genfun::FunctionNoop( &aux ) ;
295  }
296  const AbsFunction& aux = GSLSplineDeriv( *this ) ;
297  return Genfun::FunctionNoop( &aux ) ;
298  }
299  // ========================================================================
300 
301  // ========================================================================
302  FUNCTION_OBJECT_IMP( GSLSplineDeriv )
303  // ========================================================================
304 
305  // ========================================================================
326  // ========================================================================
327  GSLSplineDeriv::GSLSplineDeriv
328  ( const GSLSplineDeriv::Data1D& x ,
329  const GSLSplineDeriv::Data1D& y ,
330  const GaudiMath::Interpolation::Type type )
331  : AbsFunction ()
332  , m_spline ( x , y , type )
333  {}
334  // ========================================================================
335 
336  // ========================================================================
355  // ========================================================================
356  GSLSplineDeriv::GSLSplineDeriv
357  ( const GSLSplineDeriv::Data2D& data ,
358  const GaudiMath::Interpolation::Type type )
359  : AbsFunction ()
360  , m_spline ( data , type )
361  {}
362  // ========================================================================
363 
364  // ========================================================================
366  // ========================================================================
367  GSLSplineDeriv::GSLSplineDeriv
368  ( const SplineBase& right )
369  : AbsFunction ()
370  , m_spline ( right )
371  {}
372  // ========================================================================
373 
374  // ========================================================================
376  // ========================================================================
377  GSLSplineDeriv::GSLSplineDeriv
378  ( const GSLSplineDeriv& right )
379  : AbsFunction ()
380  , m_spline ( right )
381  {}
382  // ========================================================================
383 
384  // ========================================================================
386  // ========================================================================
387  GSLSplineDeriv::~GSLSplineDeriv(){}
388  // ========================================================================
389 
390  // ========================================================================
391  double GSLSplineDeriv::operator() ( double x ) const
392  { return m_spline.deriv ( x ) ; }
393  // ========================================================================
394  double GSLSplineDeriv::operator() ( const Argument& x ) const
395  { return m_spline.deriv ( x[0] ) ; }
396  // ========================================================================
397 
398  // ========================================================================
400  // ========================================================================
401  Genfun::Derivative GSLSplineDeriv::partial( unsigned int i ) const
402  {
403  if ( i >= 1 )
404  {
405  const AbsFunction& aux = GaudiMath::Constant( 0.0 , 1 ) ;
406  return Genfun::FunctionNoop( &aux ) ;
407  }
408  const AbsFunction& aux = GSLSplineDeriv2( *this ) ;
409  return Genfun::FunctionNoop( &aux ) ;
410  }
411  // ========================================================================
412 
413  // ========================================================================
414  FUNCTION_OBJECT_IMP( GSLSplineDeriv2 )
415  // ========================================================================
416 
417  // ========================================================================
438  // ========================================================================
439  GSLSplineDeriv2::GSLSplineDeriv2
440  ( const GSLSplineDeriv2::Data1D& x ,
441  const GSLSplineDeriv2::Data1D& y ,
442  const GaudiMath::Interpolation::Type type )
443  : AbsFunction ()
444  , m_spline ( x , y , type )
445  {}
446  // ========================================================================
447 
448  // ========================================================================
467  // ========================================================================
468  GSLSplineDeriv2::GSLSplineDeriv2
469  ( const GSLSplineDeriv2::Data2D& data ,
470  const GaudiMath::Interpolation::Type type )
471  : AbsFunction ()
472  , m_spline ( data , type )
473  {}
474  // ========================================================================
475 
476  // ========================================================================
478  // ========================================================================
479  GSLSplineDeriv2::GSLSplineDeriv2
480  ( const SplineBase& right )
481  : AbsFunction ()
482  , m_spline ( right )
483  {}
484  // ========================================================================
485 
486  // ========================================================================
488  // ========================================================================
489  GSLSplineDeriv2::GSLSplineDeriv2
490  ( const GSLSplineDeriv2& right )
491  : AbsFunction ()
492  , m_spline ( right )
493  {}
494  // ========================================================================
495 
496  // ========================================================================
498  // ========================================================================
499  GSLSplineDeriv2::~GSLSplineDeriv2(){}
500  // ========================================================================
501 
502  // ========================================================================
503  double GSLSplineDeriv2::operator() ( double x ) const
504  { return m_spline.deriv2 ( x ) ; }
505  // ========================================================================
506  double GSLSplineDeriv2::operator() ( const Argument& x ) const
507  { return m_spline.deriv2 ( x[0] ) ; }
508  // ========================================================================
509 
510  // ========================================================================
512  // ========================================================================
513  Genfun::Derivative GSLSplineDeriv2::partial( unsigned int i ) const
514  {
515  if ( i >= 1 )
516  {
517  const AbsFunction& aux = GaudiMath::Constant( 0.0 , 1 ) ;
518  return Genfun::FunctionNoop( &aux ) ;
519  }
520  const AbsFunction& aux = GaudiMath::Derivative( *this , i ) ;
521  return Genfun::FunctionNoop( &aux ) ;
522  }
523  // ========================================================================
524 
525 
526  // ========================================================================
527  FUNCTION_OBJECT_IMP( GSLSplineInteg )
528  // ========================================================================
529 
530  // ========================================================================
551  // ========================================================================
552  GSLSplineInteg::GSLSplineInteg
553  ( const GSLSplineInteg::Data1D& x ,
554  const GSLSplineInteg::Data1D& y ,
555  const GaudiMath::Interpolation::Type type ,
556  const double low )
557  : AbsFunction ( )
558  , m_spline ( x , y , type )
559  , m_low ( low )
560  {}
561  // ========================================================================
562 
563  // ========================================================================
582  // ========================================================================
583  GSLSplineInteg::GSLSplineInteg
584  ( const GSLSplineInteg::Data2D& data ,
585  const GaudiMath::Interpolation::Type type ,
586  const double low )
587  : AbsFunction ()
588  , m_spline ( data , type )
589  , m_low ( low )
590  {}
591  // ========================================================================
592 
593  // ========================================================================
595  // ========================================================================
596  GSLSplineInteg::GSLSplineInteg
597  ( const SplineBase& right ,
598  const double low )
599  : AbsFunction ()
600  , m_spline ( right )
601  , m_low ( low )
602  {}
603  // ========================================================================
604 
605  // ========================================================================
607  // ========================================================================
608  GSLSplineInteg::GSLSplineInteg
609  ( const GSLSplineInteg& right )
610  : AbsFunction ()
611  , m_spline ( right )
612  , m_low ( right.m_low )
613  {}
614  // ========================================================================
615 
616  // ========================================================================
618  // ========================================================================
619  GSLSplineInteg::~GSLSplineInteg(){}
620  // ========================================================================
621 
622  // ========================================================================
623  double GSLSplineInteg::operator() ( double x ) const
624  { return m_spline.integ ( m_low , x ) ; }
625  // ========================================================================
626  double GSLSplineInteg::operator() ( const Argument& x ) const
627  { return m_spline.integ ( m_low , x[0] ) ; }
628  // ========================================================================
629 
630  // ========================================================================
632  // ========================================================================
633  Genfun::Derivative GSLSplineInteg::partial( unsigned int i ) const
634  {
635  if ( i >= 1 )
636  {
637  const AbsFunction& aux = GaudiMath::Constant( 0.0 , 1 ) ;
638  return Genfun::FunctionNoop( &aux ) ;
639  }
640  const AbsFunction& aux = GSLSpline( *this ) ;
641  return Genfun::FunctionNoop( &aux ) ;
642  }
643  // ========================================================================
644 
645  }
646 }
Generated at Wed Nov 28 2012 12:17:12 for Gaudi Framework, version v23r5 by Doxygen version 1.8.2 written by Dimitri van Heesch, © 1997-2004