00001
00002
00003
00004
00005
00006
00007 #include "gsl/gsl_interp.h"
00008 #include "gsl/gsl_spline.h"
00009
00010
00011
00012 #include "GaudiMath/Splines.h"
00013 #include "GaudiMath/GaudiMath.h"
00014
00015
00022 namespace Genfun
00023 {
00024 namespace GaudiMathImplementation
00025 {
00031 SplineBase::SplineBase
00032 ( const SplineBase::Data1D& x ,
00033 const SplineBase::Data1D& y ,
00034 const GaudiMath::Interpolation::Type type )
00035 : m_init ( false )
00036 , m_dim ( x.size() )
00037 , m_x ( new double[ x.size() ] )
00038 , m_y ( new double[ y.size() ] )
00039 , m_spline ( 0 )
00040 , m_accel ( 0 )
00041 , m_type ( type )
00042 {
00043 std::copy( x.begin() , x.end() , m_x ) ;
00044 std::copy( y.begin() , y.end() , m_y ) ;
00045 }
00046
00047
00048
00053
00054 SplineBase::SplineBase
00055 ( const SplineBase::Data2D& data ,
00056 const GaudiMath::Interpolation::Type type )
00057 : m_init ( false )
00058 , m_dim ( data.size() )
00059 , m_x ( new double[ data.size() ] )
00060 , m_y ( new double[ data.size() ] )
00061 , m_spline ( 0 )
00062 , m_accel ( 0 )
00063 , m_type ( type )
00064 {
00065 double* _x = m_x ;
00066 double* _y = m_y ;
00067 for( Data2D::const_iterator it =
00068 data.begin() ; data.end() != it ; ++it )
00069 {
00070 *_x = it -> first ; ++_x ;
00071 *_y = it -> second ; ++_y ;
00072 }
00073 }
00074
00075
00076
00078
00079 SplineBase::SplineBase ( const SplineBase& right )
00080 : m_init ( false )
00081 , m_dim ( right.m_dim )
00082 , m_x ( new double[ right.m_dim ] )
00083 , m_y ( new double[ right.m_dim ] )
00084 , m_spline ( 0 )
00085 , m_accel ( 0 )
00086 , m_type ( right.m_type )
00087 {
00088 std::copy( right.m_x , right.m_x + right.m_dim , m_x ) ;
00089 std::copy( right.m_y , right.m_y + right.m_dim , m_y ) ;
00090 }
00091
00092
00093
00095
00096 SplineBase::~SplineBase()
00097 {
00098 if ( 0 != m_spline ) { gsl_spline_free ( m_spline ) ; }
00099 if ( 0 != m_accel ) { gsl_interp_accel_free ( m_accel ) ; }
00100
00101 if ( 0 != m_x ) { delete[] m_x ; }
00102 if ( 0 != m_y ) { delete[] m_y ; }
00103 }
00104
00105
00106
00107 void SplineBase::initialize() const
00108 {
00109 if ( m_init ) { return ; }
00110
00111 const gsl_interp_type* T = 0 ;
00112
00113 switch ( m_type )
00114 {
00115 case GaudiMath::Interpolation::Linear :
00116 T = gsl_interp_linear ; break ;
00117 case GaudiMath::Interpolation::Polynomial :
00118 T = gsl_interp_polynomial ; break ;
00119 case GaudiMath::Interpolation::Cspline :
00120 T = gsl_interp_cspline ; break ;
00121 case GaudiMath::Interpolation::Cspline_Periodic :
00122 T = gsl_interp_cspline_periodic ; break ;
00123 case GaudiMath::Interpolation::Akima :
00124 T = gsl_interp_akima ; break ;
00125 case GaudiMath::Interpolation::Akima_Periodic :
00126 T = gsl_interp_akima_periodic ; break ;
00127 default :
00128 T = gsl_interp_cspline ; break ;
00129 };
00130
00131 m_spline = gsl_spline_alloc( T , m_dim ) ;
00132
00133 gsl_spline_init( m_spline , m_x , m_y , m_dim ) ;
00134
00135 m_accel = gsl_interp_accel_alloc() ;
00136
00137 m_init = true ;
00138
00139 }
00140
00141
00142
00143 double SplineBase::eval ( const double x ) const
00144 {
00145 if ( !m_init ) { initialize() ; }
00146 return gsl_spline_eval ( m_spline , x , m_accel );
00147 }
00148
00149
00150
00151 double SplineBase::deriv ( const double x ) const
00152 {
00153 if ( !m_init ) { initialize() ; }
00154 return gsl_spline_eval_deriv ( m_spline , x , m_accel );
00155 }
00156
00157
00158
00159 double SplineBase::deriv2 ( const double x ) const
00160 {
00161 if ( !m_init ) { initialize() ; }
00162 return gsl_spline_eval_deriv2 ( m_spline , x , m_accel );
00163 }
00164
00165
00166
00167 double SplineBase::integ ( const double a ,
00168 const double b ) const
00169 {
00170 if ( !m_init ) { initialize() ; }
00171 return gsl_spline_eval_integ ( m_spline , a , b , m_accel ) ;
00172 }
00173
00174
00175
00176
00177 FUNCTION_OBJECT_IMP( GSLSpline )
00178
00179
00180
00201
00202 GSLSpline::GSLSpline
00203 ( const GSLSpline::Data1D& x ,
00204 const GSLSpline::Data1D& y ,
00205 const GaudiMath::Interpolation::Type type )
00206 : AbsFunction ()
00207 , m_spline ( x , y , type )
00208 {}
00209
00210
00211
00230
00231 GSLSpline::GSLSpline
00232 ( const GSLSpline::Data2D& data ,
00233 const GaudiMath::Interpolation::Type type )
00234 : AbsFunction ()
00235 , m_spline ( data , type )
00236 {}
00237
00238
00239
00241
00242 GSLSpline::GSLSpline
00243 ( const SplineBase& right )
00244 : AbsFunction ()
00245 , m_spline ( right )
00246 {}
00247
00248
00249
00251
00252 GSLSpline::GSLSpline
00253 ( const GSLSpline& right )
00254 : AbsFunction ()
00255 , m_spline ( right )
00256 {}
00257
00258
00259
00261
00262 GSLSpline::~GSLSpline(){}
00263
00264
00265
00266 double GSLSpline::operator() ( double x ) const
00267 { return m_spline.eval( x ) ; }
00268
00269 double GSLSpline::operator() ( const Argument& x ) const
00270 { return m_spline.eval( x[0] ) ; }
00271
00272
00273
00275
00276 Genfun::Derivative GSLSpline::partial( unsigned int i ) const
00277 {
00278 if ( i >= 1 )
00279 {
00280 const AbsFunction& aux = GaudiMath::Constant( 0.0 , 1 ) ;
00281 return Genfun::FunctionNoop( &aux ) ;
00282 }
00283 const AbsFunction& aux = GSLSplineDeriv( *this ) ;
00284 return Genfun::FunctionNoop( &aux ) ;
00285 }
00286
00287
00288
00289 FUNCTION_OBJECT_IMP( GSLSplineDeriv )
00290
00291
00292
00313
00314 GSLSplineDeriv::GSLSplineDeriv
00315 ( const GSLSplineDeriv::Data1D& x ,
00316 const GSLSplineDeriv::Data1D& y ,
00317 const GaudiMath::Interpolation::Type type )
00318 : AbsFunction ()
00319 , m_spline ( x , y , type )
00320 {}
00321
00322
00323
00342
00343 GSLSplineDeriv::GSLSplineDeriv
00344 ( const GSLSplineDeriv::Data2D& data ,
00345 const GaudiMath::Interpolation::Type type )
00346 : AbsFunction ()
00347 , m_spline ( data , type )
00348 {}
00349
00350
00351
00353
00354 GSLSplineDeriv::GSLSplineDeriv
00355 ( const SplineBase& right )
00356 : AbsFunction ()
00357 , m_spline ( right )
00358 {}
00359
00360
00361
00363
00364 GSLSplineDeriv::GSLSplineDeriv
00365 ( const GSLSplineDeriv& right )
00366 : AbsFunction ()
00367 , m_spline ( right )
00368 {}
00369
00370
00371
00373
00374 GSLSplineDeriv::~GSLSplineDeriv(){}
00375
00376
00377
00378 double GSLSplineDeriv::operator() ( double x ) const
00379 { return m_spline.deriv ( x ) ; }
00380
00381 double GSLSplineDeriv::operator() ( const Argument& x ) const
00382 { return m_spline.deriv ( x[0] ) ; }
00383
00384
00385
00387
00388 Genfun::Derivative GSLSplineDeriv::partial( unsigned int i ) const
00389 {
00390 if ( i >= 1 )
00391 {
00392 const AbsFunction& aux = GaudiMath::Constant( 0.0 , 1 ) ;
00393 return Genfun::FunctionNoop( &aux ) ;
00394 }
00395 const AbsFunction& aux = GSLSplineDeriv2( *this ) ;
00396 return Genfun::FunctionNoop( &aux ) ;
00397 }
00398
00399
00400
00401 FUNCTION_OBJECT_IMP( GSLSplineDeriv2 )
00402
00403
00404
00425
00426 GSLSplineDeriv2::GSLSplineDeriv2
00427 ( const GSLSplineDeriv2::Data1D& x ,
00428 const GSLSplineDeriv2::Data1D& y ,
00429 const GaudiMath::Interpolation::Type type )
00430 : AbsFunction ()
00431 , m_spline ( x , y , type )
00432 {}
00433
00434
00435
00454
00455 GSLSplineDeriv2::GSLSplineDeriv2
00456 ( const GSLSplineDeriv2::Data2D& data ,
00457 const GaudiMath::Interpolation::Type type )
00458 : AbsFunction ()
00459 , m_spline ( data , type )
00460 {}
00461
00462
00463
00465
00466 GSLSplineDeriv2::GSLSplineDeriv2
00467 ( const SplineBase& right )
00468 : AbsFunction ()
00469 , m_spline ( right )
00470 {}
00471
00472
00473
00475
00476 GSLSplineDeriv2::GSLSplineDeriv2
00477 ( const GSLSplineDeriv2& right )
00478 : AbsFunction ()
00479 , m_spline ( right )
00480 {}
00481
00482
00483
00485
00486 GSLSplineDeriv2::~GSLSplineDeriv2(){}
00487
00488
00489
00490 double GSLSplineDeriv2::operator() ( double x ) const
00491 { return m_spline.deriv2 ( x ) ; }
00492
00493 double GSLSplineDeriv2::operator() ( const Argument& x ) const
00494 { return m_spline.deriv2 ( x[0] ) ; }
00495
00496
00497
00499
00500 Genfun::Derivative GSLSplineDeriv2::partial( unsigned int i ) const
00501 {
00502 if ( i >= 1 )
00503 {
00504 const AbsFunction& aux = GaudiMath::Constant( 0.0 , 1 ) ;
00505 return Genfun::FunctionNoop( &aux ) ;
00506 }
00507 const AbsFunction& aux = GaudiMath::Derivative( *this , i ) ;
00508 return Genfun::FunctionNoop( &aux ) ;
00509 }
00510
00511
00512
00513
00514 FUNCTION_OBJECT_IMP( GSLSplineInteg )
00515
00516
00517
00538
00539 GSLSplineInteg::GSLSplineInteg
00540 ( const GSLSplineInteg::Data1D& x ,
00541 const GSLSplineInteg::Data1D& y ,
00542 const GaudiMath::Interpolation::Type type ,
00543 const double low )
00544 : AbsFunction ( )
00545 , m_spline ( x , y , type )
00546 , m_low ( low )
00547 {}
00548
00549
00550
00569
00570 GSLSplineInteg::GSLSplineInteg
00571 ( const GSLSplineInteg::Data2D& data ,
00572 const GaudiMath::Interpolation::Type type ,
00573 const double low )
00574 : AbsFunction ()
00575 , m_spline ( data , type )
00576 , m_low ( low )
00577 {}
00578
00579
00580
00582
00583 GSLSplineInteg::GSLSplineInteg
00584 ( const SplineBase& right ,
00585 const double low )
00586 : AbsFunction ()
00587 , m_spline ( right )
00588 , m_low ( low )
00589 {}
00590
00591
00592
00594
00595 GSLSplineInteg::GSLSplineInteg
00596 ( const GSLSplineInteg& right )
00597 : AbsFunction ()
00598 , m_spline ( right )
00599 , m_low ( right.m_low )
00600 {}
00601
00602
00603
00605
00606 GSLSplineInteg::~GSLSplineInteg(){}
00607
00608
00609
00610 double GSLSplineInteg::operator() ( double x ) const
00611 { return m_spline.integ ( m_low , x ) ; }
00612
00613 double GSLSplineInteg::operator() ( const Argument& x ) const
00614 { return m_spline.integ ( m_low , x[0] ) ; }
00615
00616
00617
00619
00620 Genfun::Derivative GSLSplineInteg::partial( unsigned int i ) const
00621 {
00622 if ( i >= 1 )
00623 {
00624 const AbsFunction& aux = GaudiMath::Constant( 0.0 , 1 ) ;
00625 return Genfun::FunctionNoop( &aux ) ;
00626 }
00627 const AbsFunction& aux = GSLSpline( *this ) ;
00628 return Genfun::FunctionNoop( &aux ) ;
00629 }
00630
00631
00632 }
00633 }
