interpolate._bsplines

Module Contents

Classes

BSpline(self,t,c,k,extrapolate=True,axis=0) rUnivariate spline in the B-spline basis.

Functions

prod(x) Product of a list of numbers; ~40x faster vs np.prod for Python tuples
_get_dtype(dtype) Return np.complex128 for complex dtypes, np.float64 otherwise.
_as_float_array(x,check_finite=False) Convert the input into a C contiguous float array.
_not_a_knot(x,k) Given data x, construct the knot vector w/ not-a-knot BC.
_augknt(x,k) Construct a knot vector appropriate for the order-k interpolation.
make_interp_spline(x,y,k=3,t=None,bc_type=None,axis=0,check_finite=True) Compute the (coefficients of) interpolating B-spline.
make_lsq_spline(x,y,t,k=3,w=None,axis=0,check_finite=True) rCompute the (coefficients of) an LSQ B-spline.
prod(x)

Product of a list of numbers; ~40x faster vs np.prod for Python tuples

_get_dtype(dtype)

Return np.complex128 for complex dtypes, np.float64 otherwise.

_as_float_array(x, check_finite=False)

Convert the input into a C contiguous float array.

NB: Upcasts half- and single-precision floats to double precision.

class BSpline(t, c, k, extrapolate=True, axis=0)

rUnivariate spline in the B-spline basis.

where are B-spline basis functions of degree k and knots t.

t : ndarray, shape (n+k+1,)
knots
c : ndarray, shape (>=n, …)
spline coefficients
k : int
B-spline order
extrapolate : bool or ‘periodic’, optional
whether to extrapolate beyond the base interval, t[k] .. t[n], or to return nans. If True, extrapolates the first and last polynomial pieces of b-spline functions active on the base interval. If ‘periodic’, periodic extrapolation is used. Default is True.
axis : int, optional
Interpolation axis. Default is zero.
t : ndarray
knot vector
c : ndarray
spline coefficients
k : int
spline degree
extrapolate : bool
If True, extrapolates the first and last polynomial pieces of b-spline functions active on the base interval.
axis : int
Interpolation axis.
tck : tuple
A read-only equivalent of (self.t, self.c, self.k)

__call__ basis_element derivative antiderivative integrate construct_fast

B-spline basis elements are defined via

Implementation details

  • At least k+1 coefficients are required for a spline of degree k, so that n >= k+1. Additional coefficients, c[j] with j > n, are ignored.
  • B-spline basis elements of degree k form a partition of unity on the base interval, t[k] <= x <= t[n].

Translating the recursive definition of B-splines into Python code, we have:

>>> def B(x, k, i, t):
...    if k == 0:
...       return 1.0 if t[i] <= x < t[i+1] else 0.0
...    if t[i+k] == t[i]:
...       c1 = 0.0
...    else:
...       c1 = (x - t[i])/(t[i+k] - t[i]) * B(x, k-1, i, t)
...    if t[i+k+1] == t[i+1]:
...       c2 = 0.0
...    else:
...       c2 = (t[i+k+1] - x)/(t[i+k+1] - t[i+1]) * B(x, k-1, i+1, t)
...    return c1 + c2
>>> def bspline(x, t, c, k):
...    n = len(t) - k - 1
...    assert (n >= k+1) and (len(c) >= n)
...    return sum(c[i] * B(x, k, i, t) for i in range(n))

Note that this is an inefficient (if straightforward) way to evaluate B-splines — this spline class does it in an equivalent, but much more efficient way.

Here we construct a quadratic spline function on the base interval 2 <= x <= 4 and compare with the naive way of evaluating the spline:

>>> from scipy.interpolate import BSpline
>>> k = 2
>>> t = [0, 1, 2, 3, 4, 5, 6]
>>> c = [-1, 2, 0, -1]
>>> spl = BSpline(t, c, k)
>>> spl(2.5)
array(1.375)
>>> bspline(2.5, t, c, k)
1.375

Note that outside of the base interval results differ. This is because BSpline extrapolates the first and last polynomial pieces of b-spline functions active on the base interval.

>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots()
>>> xx = np.linspace(1.5, 4.5, 50)
>>> ax.plot(xx, [bspline(x, t, c ,k) for x in xx], 'r-', lw=3, label='naive')
>>> ax.plot(xx, spl(xx), 'b-', lw=4, alpha=0.7, label='BSpline')
>>> ax.grid(True)
>>> ax.legend(loc='best')
>>> plt.show()
[1]Tom Lyche and Knut Morken, Spline methods, http://www.uio.no/studier/emner/matnat/ifi/INF-MAT5340/v05/undervisningsmateriale/
[2]Carl de Boor, A practical guide to splines, Springer, 2001.
__init__(t, c, k, extrapolate=True, axis=0)
construct_fast(t, c, k, extrapolate=True, axis=0)

Construct a spline without making checks.

Accepts same parameters as the regular constructor. Input arrays t and c must of correct shape and dtype.

tck()

Equvalent to (self.t, self.c, self.k) (read-only).

basis_element(t, extrapolate=True)

Return a B-spline basis element B(x | t[0], ..., t[k+1]).

t : ndarray, shape (k+1,)
internal knots
extrapolate : bool or ‘periodic’, optional
whether to extrapolate beyond the base interval, t[0] .. t[k+1], or to return nans. If ‘periodic’, periodic extrapolation is used. Default is True.
basis_element : callable
A callable representing a B-spline basis element for the knot vector t.

The order of the b-spline, k, is inferred from the length of t as len(t)-2. The knot vector is constructed by appending and prepending k+1 elements to internal knots t.

Construct a cubic b-spline:

>>> from scipy.interpolate import BSpline
>>> b = BSpline.basis_element([0, 1, 2, 3, 4])
>>> k = b.k
>>> b.t[k:-k]
array([ 0.,  1.,  2.,  3.,  4.])
>>> k
3

Construct a second order b-spline on [0, 1, 1, 2], and compare to its explicit form:

>>> t = [-1, 0, 1, 1, 2]
>>> b = BSpline.basis_element(t[1:])
>>> def f(x):
...     return np.where(x < 1, x*x, (2. - x)**2)
>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots()
>>> x = np.linspace(0, 2, 51)
>>> ax.plot(x, b(x), 'g', lw=3)
>>> ax.plot(x, f(x), 'r', lw=8, alpha=0.4)
>>> ax.grid(True)
>>> plt.show()
__call__(x, nu=0, extrapolate=None)

Evaluate a spline function.

x : array_like
points to evaluate the spline at.
nu: int, optional
derivative to evaluate (default is 0).
extrapolate : bool or ‘periodic’, optional
whether to extrapolate based on the first and last intervals or return nans. If ‘periodic’, periodic extrapolation is used. Default is self.extrapolate.
y : array_like
Shape is determined by replacing the interpolation axis in the coefficient array with the shape of x.
_evaluate(xp, nu, extrapolate, out)
_ensure_c_contiguous()

c and t may be modified by the user. The Cython code expects that they are C contiguous.

derivative(nu=1)

Return a b-spline representing the derivative.

nu : int, optional
Derivative order. Default is 1.
b : BSpline object
A new instance representing the derivative.

splder, splantider

antiderivative(nu=1)

Return a b-spline representing the antiderivative.

nu : int, optional
Antiderivative order. Default is 1.
b : BSpline object
A new instance representing the antiderivative.

If antiderivative is computed and self.extrapolate='periodic', it will be set to False for the returned instance. This is done because the antiderivative is no longer periodic and its correct evaluation outside of the initially given x interval is difficult.

splder, splantider

integrate(a, b, extrapolate=None)

Compute a definite integral of the spline.

a : float
Lower limit of integration.
b : float
Upper limit of integration.
extrapolate : bool or ‘periodic’, optional
whether to extrapolate beyond the base interval, t[k] .. t[-k-1], or take the spline to be zero outside of the base interval. If ‘periodic’, periodic extrapolation is used. If None (default), use self.extrapolate.
I : array_like
Definite integral of the spline over the interval [a, b].

Construct the linear spline x if x < 1 else 2 - x on the base interval , and integrate it

>>> from scipy.interpolate import BSpline
>>> b = BSpline.basis_element([0, 1, 2])
>>> b.integrate(0, 1)
array(0.5)

If the integration limits are outside of the base interval, the result is controlled by the extrapolate parameter

>>> b.integrate(-1, 1)
array(0.0)
>>> b.integrate(-1, 1, extrapolate=False)
array(0.5)
>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots()
>>> ax.grid(True)
>>> ax.axvline(0, c='r', lw=5, alpha=0.5)  # base interval
>>> ax.axvline(2, c='r', lw=5, alpha=0.5)
>>> xx = [-1, 1, 2]
>>> ax.plot(xx, b(xx))
>>> plt.show()
_not_a_knot(x, k)

Given data x, construct the knot vector w/ not-a-knot BC. cf de Boor, XIII(12).

_augknt(x, k)

Construct a knot vector appropriate for the order-k interpolation.

make_interp_spline(x, y, k=3, t=None, bc_type=None, axis=0, check_finite=True)

Compute the (coefficients of) interpolating B-spline.

x : array_like, shape (n,)
Abscissas.
y : array_like, shape (n, …)
Ordinates.
k : int, optional
B-spline degree. Default is cubic, k=3.
t : array_like, shape (nt + k + 1,), optional.
Knots. The number of knots needs to agree with the number of datapoints and the number of derivatives at the edges. Specifically, nt - n must equal len(deriv_l) + len(deriv_r).
bc_type : 2-tuple or None
Boundary conditions. Default is None, which means choosing the boundary conditions automatically. Otherwise, it must be a length-two tuple where the first element sets the boundary conditions at x[0] and the second element sets the boundary conditions at x[-1]. Each of these must be an iterable of pairs (order, value) which gives the values of derivatives of specified orders at the given edge of the interpolation interval.
axis : int, optional
Interpolation axis. Default is 0.
check_finite : bool, optional
Whether to check that the input arrays contain only finite numbers. Disabling may give a performance gain, but may result in problems (crashes, non-termination) if the inputs do contain infinities or NaNs. Default is True.

b : a BSpline object of the degree k and with knots t.

Use cubic interpolation on Chebyshev nodes:

>>> def cheb_nodes(N):
...     jj = 2.*np.arange(N) + 1
...     x = np.cos(np.pi * jj / 2 / N)[::-1]
...     return x
>>> x = cheb_nodes(20)
>>> y = np.sqrt(1 - x**2)
>>> from scipy.interpolate import BSpline, make_interp_spline
>>> b = make_interp_spline(x, y)
>>> np.allclose(b(x), y)
True

Note that the default is a cubic spline with a not-a-knot boundary condition

>>> b.k
3

Here we use a ‘natural’ spline, with zero 2nd derivatives at edges:

>>> l, r = [(2, 0)], [(2, 0)]
>>> b_n = make_interp_spline(x, y, bc_type=(l, r))
>>> np.allclose(b_n(x), y)
True
>>> x0, x1 = x[0], x[-1]
>>> np.allclose([b_n(x0, 2), b_n(x1, 2)], [0, 0])
True

Interpolation of parametric curves is also supported. As an example, we compute a discretization of a snail curve in polar coordinates

>>> phi = np.linspace(0, 2.*np.pi, 40)
>>> r = 0.3 + np.cos(phi)
>>> x, y = r*np.cos(phi), r*np.sin(phi)  # convert to Cartesian coordinates

Build an interpolating curve, parameterizing it by the angle

>>> from scipy.interpolate import make_interp_spline
>>> spl = make_interp_spline(phi, np.c_[x, y])

Evaluate the interpolant on a finer grid (note that we transpose the result to unpack it into a pair of x- and y-arrays)

>>> phi_new = np.linspace(0, 2.*np.pi, 100)
>>> x_new, y_new = spl(phi_new).T

Plot the result

>>> import matplotlib.pyplot as plt
>>> plt.plot(x, y, 'o')
>>> plt.plot(x_new, y_new, '-')
>>> plt.show()

BSpline : base class representing the B-spline objects CubicSpline : a cubic spline in the polynomial basis make_lsq_spline : a similar factory function for spline fitting UnivariateSpline : a wrapper over FITPACK spline fitting routines splrep : a wrapper over FITPACK spline fitting routines

make_lsq_spline(x, y, t, k=3, w=None, axis=0, check_finite=True)

rCompute the (coefficients of) an LSQ B-spline.

The result is a linear combination

of the B-spline basis elements, , which minimizes

x : array_like, shape (m,)
Abscissas.
y : array_like, shape (m, …)
Ordinates.
t : array_like, shape (n + k + 1,).
Knots. Knots and data points must satisfy Schoenberg-Whitney conditions.
k : int, optional
B-spline degree. Default is cubic, k=3.
w : array_like, shape (n,), optional
Weights for spline fitting. Must be positive. If None, then weights are all equal. Default is None.
axis : int, optional
Interpolation axis. Default is zero.
check_finite : bool, optional
Whether to check that the input arrays contain only finite numbers. Disabling may give a performance gain, but may result in problems (crashes, non-termination) if the inputs do contain infinities or NaNs. Default is True.

b : a BSpline object of the degree k with knots t.

The number of data points must be larger than the spline degree k.

Knots t must satisfy the Schoenberg-Whitney conditions, i.e., there must be a subset of data points x[j] such that t[j] < x[j] < t[j+k+1], for j=0, 1,...,n-k-2.

Generate some noisy data:

>>> x = np.linspace(-3, 3, 50)
>>> y = np.exp(-x**2) + 0.1 * np.random.randn(50)

Now fit a smoothing cubic spline with a pre-defined internal knots. Here we make the knot vector (k+1)-regular by adding boundary knots:

>>> from scipy.interpolate import make_lsq_spline, BSpline
>>> t = [-1, 0, 1]
>>> k = 3
>>> t = np.r_[(x[0],)*(k+1),
...           t,
...           (x[-1],)*(k+1)]
>>> spl = make_lsq_spline(x, y, t, k)

For comparison, we also construct an interpolating spline for the same set of data:

>>> from scipy.interpolate import make_interp_spline
>>> spl_i = make_interp_spline(x, y)

Plot both:

>>> import matplotlib.pyplot as plt
>>> xs = np.linspace(-3, 3, 100)
>>> plt.plot(x, y, 'ro', ms=5)
>>> plt.plot(xs, spl(xs), 'g-', lw=3, label='LSQ spline')
>>> plt.plot(xs, spl_i(xs), 'b-', lw=3, alpha=0.7, label='interp spline')
>>> plt.legend(loc='best')
>>> plt.show()

NaN handling: If the input arrays contain nan values, the result is not useful since the underlying spline fitting routines cannot deal with nan. A workaround is to use zero weights for not-a-number data points:

>>> y[8] = np.nan
>>> w = np.isnan(y)
>>> y[w] = 0.
>>> tck = make_lsq_spline(x, y, t, w=~w)

Notice the need to replace a nan by a numerical value (precise value does not matter as long as the corresponding weight is zero.)

BSpline : base class representing the B-spline objects make_interp_spline : a similar factory function for interpolating splines LSQUnivariateSpline : a FITPACK-based spline fitting routine splrep : a FITPACK-based fitting routine