#!/usr/bin/env python
# Author : Pierre Schnizer
"""
Wrapper for the interpolations of gsl. This solver wraps all features as described
in Chapter 26 of the gsl documentation.
Difference between spline and interpolation module:
--------------------------------------------------
In the interpolation module the data for the independent and dependent data are
kept as reference in the various objects,whereas the spline module copies these
data into the internal C gsl_spline struct.
"""
import gslwrap
class _acceleration:
"""
This class is meant to be used together with the interpolation and the spline.
It expects the derived class to define the array _xa.
"""
def _accel_alloc (self, *args): return gslwrap.gsl_interp_accel_alloc(*args)
def _accel_find (self, *args): return gslwrap.gsl_interp_accel_find(*args)
def _accel_free (self, *args): return gslwrap.gsl_interp_accel_free(*args)
def _accel_reset (self, *args): return gslwrap.gsl_interp_accel_reset(*args)
def __init__(self):
self._accel = None
assert(self._accel_alloc != None)
assert(self._accel_free != None)
self._accel= self._accel_alloc()
def __del__(self):
if hasattr(self, '_accel'):
if self._accel != None:
self._accel_free(self._accel)
def accel_reset(self):
self._accel_reset(self._accel)
def accel_find(self, x):
"""
This method performs a lookup action on the data array X_ARRAY
of size SIZE, using the given accelerator A. This is how lookups
are performed during evaluation of an interpolation. The function
returns an index i such that `xarray[i] <= x < xarray[i+1]'.
"""
return self._accel_find(self._accel, self._xa, x)
def bsearch(self, x, index_lo, index_high):
"""
input : x, index_lo, index_high
This function returns the index i of the array X_ARRAY such that
`x_array[i] <= x < x_array[i+1]'. The index is searched for in
the range [INDEX_LO,INDEX_HI].
"""
return self._bsearch(self._xa, x, index_lo, index_high)
class _common(_acceleration):
def _alloc (self, *args): return gslwrap.gsl_interp_alloc(*args)
def _free (self, *args): return gslwrap.gsl_interp_free(*args)
def _name (self, *args): return gslwrap.gsl_interp_name(*args)
def _bsearch (self, *args): return gslwrap.gsl_interp_bsearch(*args)
def _eval (self, *args): return gslwrap.gsl_interp_eval(*args)
def _eval_deriv (self, *args): return gslwrap.gsl_interp_eval_deriv(*args)
def _eval_deriv2 (self, *args): return gslwrap.gsl_interp_eval_deriv2(*args)
def _eval_deriv2_e(self, *args): return gslwrap.gsl_interp_eval_deriv2_e(*args)
def _eval_deriv_e (self, *args): return gslwrap.gsl_interp_eval_deriv_e(*args)
def _eval_e (self, *args): return gslwrap.gsl_interp_eval_e(*args)
def _eval_integ (self, *args): return gslwrap.gsl_interp_eval_integ(*args)
def _eval_integ_e (self, *args): return gslwrap.gsl_interp_eval_integ_e(*args)
def _init (self, *args): return gslwrap.gsl_interp_init(*args)
def _min_size (self, *args): return gslwrap.gsl_interp_min_size(*args)
def _name (self, *args): return gslwrap.gsl_interp_name(*args)
_type = None
def __init__(self, size):
self._ptr = None
assert(self._alloc != None)
assert(self._free != None)
assert(self._type != None)
self._ptr = self._alloc(self._type, size)
_acceleration.__init__(self)
self._xa = None
self._xyptr = None
def __del__(self):
if hasattr(self, '_ptr'):
if self._ptr != None:
self._free(self._ptr)
_acceleration.__del__(self)
def init(self, xa, ya):
"""
input : xa, ya
xa ... array of independent values
ya ... array of dependent values
This method initializes this for the
data (xa,ay) where ya and ya are arrays of the size, which was set,
when this object was initialised. The
interpolation object only keeps references to the data arrays
xa and ya and and stores the static state computed from the data.
The ya data array is always assumed to be strictly ordered; the
behavior for other arrangements is not defined.
"""
self._init(self._ptr, (xa,ya))
self.accel_reset()
self._xa = xa
self._ya = ya
self._xyptr = (self._ptr, xa, ya)
def eval(self, x):
"""
input : x
x ... value of the independent variable
output : y
y ... returns the interpolated value at x
"""
return apply(self._eval, self._xyptr + (x, self._accel))
def eval_e(self, x):
"""
input : x
x ... value of the independent variable
output : flag, y
flag error flag
y ... returns the interpolated value at x
"""
return apply(self._eval_e, self._xyptr + (x, self._accel))
def eval_deriv(self, x):
"""
input : x
x ... value of the independent variable
output : y
y ... returns the value of the derivative at x
"""
return apply(self._eval_deriv, self._xyptr + (x, self._accel))
def eval_deriv_e(self, x):
"""
input : x
x ... value of the independent variable
output : flag, y
flag error flag
y ... returns the value of the derivative at x
"""
return apply(self._eval_deriv_e, self._xyptr + (x, self._accel))
def eval_deriv2(self, x):
"""
input : x
x ... value of the independent variable
output : y
y ... the value of the second derivative at x
"""
return apply(self._eval_deriv2, self._xyptr + (x, self._accel))
def eval_deriv2_e(self, x):
"""
input : x
x ... value of the independent variable
output : flag, y
flag error flag
y ... the value of the second derivative at x
"""
return apply(self._eval_deriv2_e, self._xyptr + (x, self._accel))
def eval_integ(self, a, b):
"""
input : a, b,
a ... lower boundary
b ... upper boundary
output : y
y ... the integral of the object over the range [a,b]
"""
return apply(self._eval_integ, self._xyptr + (a,b, self._accel))
def eval_integ_e(self, a, b):
"""
input : a, b,
a ... lower boundary
b ... upper boundary
output : flag, y
flag error flag
y ... the integral of the object over the range [a,b]
"""
return apply(self._eval_integ_e, self._xyptr + (a,b, self._accel))
class _interpolation(_common):
def min_size(self):
"""
This function returns the minimum number of points required by the
interpolation.
"""
return self._min_size(self._ptr)
def name(self):
"""
Returns the name of the interpolation type used
"""
return self._name(self._ptr)
class linear(_interpolation):
"""
Linear interpolation.
"""
_type = gslwrap.cvar.gsl_interp_linear
class polynomial(_interpolation):
"""
Polynomial interpolation. This method should only be used for
interpolating small numbers of points because polynomial interpolation
introduces large oscillations, even for well-behaved datasets. The number
of terms in the interpolating polynomial is equal to the number of points.
"""
_type = gslwrap.cvar.gsl_interp_polynomial
class cspline(_interpolation):
"""
Cubic spline with natural boundary conditions.
"""
_type = gslwrap.cvar.gsl_interp_cspline
class cspline_periodic(_interpolation):
"""
Cubic spline with periodic boundary conditions
"""
_type = gslwrap.cvar.gsl_interp_cspline_periodic
class akima(_interpolation):
"""
Akima spline with natural boundary conditions
"""
_type = gslwrap.cvar.gsl_interp_akima
class akima_periodic(_interpolation):
"""
Akima spline with periodic boundary conditions
"""
_type = gslwrap.cvar.gsl_interp_akima_periodic
syntax highlighted by Code2HTML, v. 0.9.1