[box type=”note” align=”” class=”” width=””]This article is an excerpt from a book co-authored by L. Felipe Martins, Ruben Oliva Ramos and V Kishore Ayyadevara titled SciPy Recipes. This book provides numerous recipes in mastering common tasks related to SciPy and associated libraries such as NumPy, pandas, and matplotlib.[/box]
In today’s tutorial, we will see how to compute and solve polynomial, univariate interpolations using SciPy with detailed process and instructions.
In this recipe, we will look at how to compute data polynomial interpolation by applying some important methods which are discussed in detail in the coming How to do it… section.
Getting ready
We will need to follow some instructions and install the prerequisites.
How to do it…
Let’s get started. In the following steps, we will explain how to compute a polynomial interpolation and the things we need to know:
They require the following parameters:
- points: An ndarray of floats, shape (n, D) data point coordinates. It can be either an array of shape (n, D) or a tuple of ndim arrays.
- values: An ndarray of float or complex shape (n,) data values.
- xi: A 2D ndarray of float or tuple of 1D array, shape (M, D). Points at which to interpolate data.
- method: A {‘linear’, ‘nearest’, ‘cubic’}—This is an optional method of interpolation.
One of the nearest return value is at the data point closest to the point of interpolation. See NearestNDInterpolator for more details.
- linear tessellates the input point set to n-dimensional simplices, and interpolates linearly on each simplex. See LinearNDInterpolator for more details.
- cubic (1D): Returns the value determined from a cubic spline.
- cubic (2D): Returns the value determined from a piecewise cubic, continuously differentiable (C1), and approximately curvature-minimizing polynomial surface. See CloughTocher2DInterpolator for more details.
- fill_value: float; optional. It is the value used to fill in for requested points outside of the convex hull of the input points. If it is not provided, then the default is nan. This option has no effect on the nearest method.
- rescale: bool; optional.
Rescale points to the unit cube before performing interpolation. This is useful if some of the input dimensions have non-commensurable units and differ by many orders of magnitude.
How it works…
One can see that the exact result is reproduced by all of the methods to some degree, but for this smooth function, the piecewise cubic interpolant gives the best results:
import matplotlib.pyplot as plt
import numpy as np
methods = [None, 'none', 'nearest', 'bilinear', 'bicubic', 'spline16',
'spline36', 'hanning', 'hamming', 'hermite', 'kaiser',
'quadric',
'catrom', 'gaussian', 'bessel', 'mitchell', 'sinc', 'lanczos']
# Fixing random state for reproducibility
np.random.seed(19680801)
grid = np.random.rand(4, 4)
fig, axes = plt.subplots(3, 6, figsize=(12, 6),
subplot_kw={'xticks': [], 'yticks': []})
fig.subplots_adjust(hspace=0.3, wspace=0.05)
for ax, interp_method in zip(axes.flat, methods):
ax.imshow(grid, interpolation=interp_method, cmap='viridis')
ax.set_title(interp_method)
plt.show()This is the result of the execution:
Univariate interpolation
In the next section, we will look at how to solve univariate interpolation.
Getting ready
We will need to follow some instructions and install the prerequisites.
How to do it…
The following table summarizes the different univariate interpolation modes coded in SciPy, together with the processes that we may use to resolve them:
Finding a cubic spline that interpolates a set of data
In this recipe, we will look at how to find a cubic spline that interpolates with the main method of spline.
Getting ready
We will need to follow some instructions and install the prerequisites.
How to do it…
We can use the following functions to solve the problems with this parameter:
- x: array_like, shape (n,). A 1D array containing values of the independent variable. The values must be real, finite, and in strictly increasing order.
- y: array_like. An array containing values of the dependent variable. It can have an arbitrary number of dimensions, but the length along axis must match the length of x. The values must be finite.
- axis: int; optional. The axis along which y is assumed to be varying, meaning for x[i], the corresponding values are np.take(y, i, axis=axis). The default is 0.
- bc_type: String or two-tuple; optional. Boundary condition type. Two additional equations, given by the boundary conditions, are required to determine all coefficients of polynomials on each segment. Refer to: https://docs.scipy.org/doc/scipy-0.19.1/reference/generated/scipy.interpolate.CubicSpline.html#r59.
If bc_type is a string, then the specified condition will be applied at both ends of a spline. The available conditions are:
- not-a-knot (default): The first and second segment at a curve end are the same polynomial. This is a good default when there is no information about boundary conditions.
- periodic: The interpolated function is assumed to be periodic in the period x[-1] – x[0]. The first and last value of y must be identical: y[0] == y[-1]. This boundary condition will result in y'[0] == y'[-1] and y”[0] == y”[-1].
- clamped: The first derivatives at the curve ends are zero. Assuming there is a 1D y, bc_type=((1, 0.0), (1, 0.0)) is the same condition.
- natural: The second derivatives at the curve ends are zero. Assuming there is a 1D y, bc_type=((2, 0.0), (2, 0.0)) is the same condition.
If bc_type is two-tuple, the first and the second value will be applied at the curve’s start and end respectively. The tuple value can be one of the previously mentioned strings (except periodic) or a tuple (order, deriv_values), allowing us to specify arbitrary derivatives at curve ends:
- order: The derivative order; it is 1 or 2.
- deriv_value: An array_like containing derivative values. The shape must be the same as y, excluding the axis dimension. For example, if y is 1D, then deriv_value must be a scalar. If y is 3D with shape (n0, n1, n2) and axis=2, then deriv_value must be 2D and have the shape (n0, n1).
- extrapolate: {bool, ‘periodic’, None}; optional. bool, determines whether or not to extrapolate to out-of-bounds points based on first and last intervals, or to return NaNs. periodic, periodic extrapolation is used. If none (default), extrapolate is set to periodic for bc_type=’periodic’ and to True otherwise.
How it works…
We have the following example:
%pylab inline
from scipy.interpolate import CubicSpline
import matplotlib.pyplot as plt
x = np.arange(10)
y = np.sin(x)
cs = CubicSpline(x, y)
xs = np.arange(-0.5, 9.6, 0.1)
plt.figure(figsize=(6.5, 4))
plt.plot(x, y, 'o', label='data')
plt.plot(xs, np.sin(xs), label='true')
plt.plot(xs, cs(xs), label="S")
plt.plot(xs, cs(xs, 1), label="S'")
plt.plot(xs, cs(xs, 2), label="S''")
plt.plot(xs, cs(xs, 3), label="S'''")
plt.xlim(-0.5, 9.5)
plt.legend(loc='lower left', ncol=2)
plt.show()We can see the result here:
We see the next example:
theta = 2 * np.pi * np.linspace(0, 1, 5)
y = np.c_[np.cos(theta), np.sin(theta)]
cs = CubicSpline(theta, y, bc_type='periodic')
print("ds/dx={:.1f} ds/dy={:.1f}".format(cs(0, 1)[0], cs(0, 1)[1]))
x=0.0 ds/dy=1.0
xs = 2 * np.pi * np.linspace(0, 1, 100)
plt.figure(figsize=(6.5, 4))
plt.plot(y[:, 0], y[:, 1], 'o', label='data')
plt.plot(np.cos(xs), np.sin(xs), label='true')
plt.plot(cs(xs)[:, 0], cs(xs)[:, 1], label='spline')
plt.axes().set_aspect('equal')
plt.legend(loc='center')
plt.show()In the following screenshot, we can see the final result:
Defining a B-spline for a given set of control points
In the next section, we will look at how to solve B-splines given some controlled data.
Getting ready
We need to follow some instructions and install the prerequisites.
How to do it…
- Univariate the spline in the B-spline basis
- Execute the following: S(x)=∑j=0n-1cjBj,k;t(x)S(x)=∑j=0n-1cjBj,k;t(x)
- Where it’s Bj,k;tBj,k;t are B-spline basis functions of degree k and knots t
- We can use the following parameters:
How it works …
Here, we construct a quadratic spline function on the base interval 2 <= x <= 4 and compare it with the naive way of evaluating the spline:
from scipy import interpolate
import numpy as np
import matplotlib.pyplot as plt
# sampling
x = np.linspace(0, 10, 10)
y = np.sin(x)
# spline trough all the sampled points
tck = interpolate.splrep(x, y)
x2 = np.linspace(0, 10, 200)
y2 = interpolate.splev(x2, tck)
# spline with all the middle points as knots (not working yet)
# knots = x[1:-1] # it should be something like this
knots = np.array([x[1]]) # not working with above line and just seeing what this line does
weights = np.concatenate(([1],np.ones(x.shape[0]-2)*.01,[1]))
tck = interpolate.splrep(x, y, t=knots, w=weights)
x3 = np.linspace(0, 10, 200)
y3 = interpolate.splev(x2, tck)
# plot
plt.plot(x, y, 'go', x2, y2, 'b', x3, y3,'r')
plt.show()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.
This is the result of solving the problem:
We successfully compute numerical computation and find interpolation function using polynomial and univariate interpolation coded in SciPy.
If you found this tutorial useful, do check out the book SciPy Recipes to get quick recipes for performing other mathematical operations like differential equation, K-means and Discrete Fourier Transform.
