Gaudi Framework, version v20r4
Generated: 8 Jan 2009
|
|
|
Generated: 8 Jan 2009 |
00001 // $Id: Splines.cpp,v 1.2 2005/11/25 10:27:03 mato Exp $ 00002 // ============================================================================ 00003 // Include files 00004 // ============================================================================ 00005 // GSL 00006 // ============================================================================ 00007 #include "gsl/gsl_interp.h" 00008 #include "gsl/gsl_spline.h" 00009 // ============================================================================ 00010 // GaudiMath 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 ; } // 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 };