.. currentmodule:: scipy.interpolate .. _tutorial-extrapolation: ============================= Extrapolation tips and tricks ============================= Handling of extrapolation---evaluation of the interpolators on query points outside of the domain of interpolated data---is not fully consistent among different routines in `scipy.interpolate`. Different interpolators use different sets of keyword arguments to control the behavior outside of the data domain: some use ``extrapolate=True/False/None``, some allow the ``fill_value`` keyword. Refer to the API documentation for details for each specific interpolation routine. Depending on a particular problem, the available keywords may or may not be sufficient. Special attention needs to be paid to extrapolation of non-linear interpolants. Very often the extrapolated results make less sense with increasing distance from the data domain. This is of course to be expected: an interpolant only knows the data within the data domain. When the default extrapolated results are not adequate, users need to implement the desired extrapolation mode themselves. In this tutorial, we consider several worked examples where we demonstrate both the use of available keywords and manual implementation of desired extrapolation modes. These examples may or may not be applicable to your particular problem; they are not necessarily best practices; and they are deliberately pared down to bare essentials needed to demonstrate the main ideas, in a hope that they serve as an inspiration for your handling of your particular problem. .. _tutorial-extrapolation-left-right: `interp1d` : replicate `numpy.interp` left and right fill values ==================================================================== TL;DR: Use ``fill_value=(left, right)`` `numpy.interp` uses constant extrapolation, and defaults to extending the first and last values of the ``y`` array in the interpolation interval: the output of ``np.interp(xnew, x, y)`` is ``y[0]`` for ``xnew < x[0]`` and ``y[-1]`` for ``xnew > x[-1]``. By default, `interp1d` refuses to extrapolate, and raises a ``ValueError`` when evaluated on a data point outside of the interpolation range. This can be switched off by the ``bounds_error=False`` argument: then ``interp1d`` sets the out-of-range values with the ``fill_value``, which is ``nan`` by default. To mimic the behavior of `numpy.interp` with `interp1d`, you can use the fact that it supports a 2-tuple as the ``fill_value``. The tuple elements are then used to fill for ``xnew < min(x)`` and ``x > max(x)``, respectively. For multidimensional ``y``, these elements must have the same shape as ``y`` or be broadcastable to it. To illustrate: .. plot:: import numpy as np import matplotlib.pyplot as plt from scipy.interpolate import interp1d x = np.linspace(0, 1.5*np.pi, 11) y = np.column_stack((np.cos(x), np.sin(x))) # y.shape is (11, 2) func = interp1d(x, y, axis=0, # interpolate along columns bounds_error=False, kind='linear', fill_value=(y[0], y[-1])) xnew = np.linspace(-np.pi, 2.5*np.pi, 51) ynew = func(xnew) fix, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4)) ax1.plot(xnew, ynew[:, 0]) ax1.plot(x, y[:, 0], 'o') ax2.plot(xnew, ynew[:, 1]) ax2.plot(x, y[:, 1], 'o') plt.tight_layout() .. _tutorial-extrapolation-cubicspline-extend: CubicSpline extend the boundary conditions ========================================== `CubicSpline` needs two extra boundary conditions, which are controlled by the ``bc_type`` parameter. This parameter can either list explicit values of derivatives at the edges, or use helpful aliases. For instance, ``bc_type="clamped"`` sets the first derivatives to zero, ``bc_type="natural"`` sets the second derivatives to zero (two other recognized string values are "periodic" and "not-a-knot") While the extrapolation is controlled by the boundary condition, the relation is not very intuitive. For instance, one can expect that for ``bc_type="natural"``, the extrapolation is linear. This expectation is too strong: each boundary condition sets the derivatives at a single point, *at the boundary* only. Extrapolation is done from the first and last polynomial pieces, which — for a natural spline — is a cubic with a zero second derivative at a given point. One other way of seeing why this expectation is too strong is to consider a dataset with only three data points, where the spline has two polynomial pieces. To extrapolate linearly, this expectation implies that both of these pieces are linear. But then, two linear pieces cannot match at a middle point with a continuous 2nd derivative! (Unless of course, if all three data points actually lie on a single straight line). To illustrate the behavior we consider a synthetic dataset and compare several boundary conditions: .. plot:: import numpy as np import matplotlib.pyplot as plt from scipy.interpolate import CubicSpline xs = [1, 2, 3, 4, 5, 6, 7, 8] ys = [4.5, 3.6, 1.6, 0.0, -3.3, -3.1, -1.8, -1.7] notaknot = CubicSpline(xs, ys, bc_type='not-a-knot') natural = CubicSpline(xs, ys, bc_type='natural') clamped = CubicSpline(xs, ys, bc_type='clamped') xnew = np.linspace(min(xs) - 4, max(xs) + 4, 101) splines = [notaknot, natural, clamped] titles = ['not-a-knot', 'natural', 'clamped'] fig, axs = plt.subplots(3, 3, figsize=(12, 12)) for i in [0, 1, 2]: for j, spline, title in zip(range(3), splines, titles): axs[i, j].plot(xs, spline(xs, nu=i),'o') axs[i, j].plot(xnew, spline(xnew, nu=i),'-') axs[i, j].set_title(f'{title}, deriv={i}') plt.tight_layout() plt.show() It is clearly seen that the natural spline does have the zero second derivative at the boundaries, but extrapolation is non-linear. ``bc_type="clamped"`` shows a similar behavior: first derivatives are only equal to zero exactly at the boundary. In all cases, extrapolation is done by extending the first and last polynomial pieces of the spline, whatever they happen to be. One possible way to force the extrapolation is to extend the interpolation domain to add first and last polynomial pieces which have desired properties. Here we use ``extend`` method of the `CubicSpline` superclass, `PPoly`, to add two extra breakpoints and to make sure that the additional polynomial pieces maintain the values of the derivatives. Then the extrapolation proceeds using these two additional intervals. .. plot:: import numpy as np import matplotlib.pyplot as plt from scipy.interpolate import CubicSpline def add_boundary_knots(spline): """ Add knots infinitesimally to the left and right. Additional intervals are added to have zero 2nd and 3rd derivatives, and to maintain the first derivative from whatever boundary condition was selected. The spline is modified in place. """ # determine the slope at the left edge leftx = spline.x[0] lefty = spline(leftx) leftslope = spline(leftx, nu=1) # add a new breakpoint just to the left and use the # known slope to construct the PPoly coefficients. leftxnext = np.nextafter(leftx, leftx - 1) leftynext = lefty + leftslope*(leftxnext - leftx) leftcoeffs = np.array([0, 0, leftslope, leftynext]) spline.extend(leftcoeffs[..., None], np.r_[leftxnext]) # repeat with additional knots to the right rightx = spline.x[-1] righty = spline(rightx) rightslope = spline(rightx,nu=1) rightxnext = np.nextafter(rightx, rightx + 1) rightynext = righty + rightslope * (rightxnext - rightx) rightcoeffs = np.array([0, 0, rightslope, rightynext]) spline.extend(rightcoeffs[..., None], np.r_[rightxnext]) xs = [1, 2, 3, 4, 5, 6, 7, 8] ys = [4.5, 3.6, 1.6, 0.0, -3.3, -3.1, -1.8, -1.7] notaknot = CubicSpline(xs,ys, bc_type='not-a-knot') # not-a-knot does not require additional intervals natural = CubicSpline(xs,ys, bc_type='natural') # extend the natural natural spline with linear extrapolating knots add_boundary_knots(natural) clamped = CubicSpline(xs,ys, bc_type='clamped') # extend the clamped spline with constant extrapolating knots add_boundary_knots(clamped) xnew = np.linspace(min(xs) - 5, max(xs) + 5, 201) fig, axs = plt.subplots(3, 3,figsize=(12,12)) splines = [notaknot, natural, clamped] titles = ['not-a-knot', 'natural', 'clamped'] for i in [0, 1, 2]: for j, spline, title in zip(range(3), splines, titles): axs[i, j].plot(xs, spline(xs, nu=i),'o') axs[i, j].plot(xnew, spline(xnew, nu=i),'-') axs[i, j].set_title(f'{title}, deriv={i}') plt.tight_layout() plt.show() .. _tutorial-extrapolation-asymptotics: Manually implement the asymptotics ================================== The previous trick of extending the interpolation domain relies on the `CubicSpline.extend` method. A somewhat more general alternative is to implement a wrapper which handles the out-of-bounds behavior explicitly. Let us consider a worked example. The setup ~~~~~~~~~ Suppose we want to solve at a given value of :math:`a` the equation .. math:: a x = 1/\tan{x}\;. (One application where these kinds of equations appear is solving for energy levels of a quantum particle). For simplicity, let’s only consider :math:`x\in (0, \pi/2)`. Solving this equation *once* is straightforward: .. plot:: import numpy as np import matplotlib.pyplot as plt from scipy.optimize import brentq def f(x, a): return a*x - 1/np.tan(x) a = 3 x0 = brentq(f, 1e-16, np.pi/2, args=(a,)) # here we shift the left edge # by a machine epsilon to avoid # a division by zero at x=0 xx = np.linspace(0.2, np.pi/2, 101) plt.plot(xx, a*xx, '--') plt.plot(xx, 1/np.tan(xx), '--') plt.plot(x0, a*x0, 'o', ms=12) plt.text(0.1, 0.9, fr'$x_0 = {x0:.3f}$', transform=plt.gca().transAxes, fontsize=16) plt.show() However, if we need to solve it multiple times (e.g. to find *a series* of roots due to periodicity of the ``tan`` function), repeated calls to `scipy.optimize.brentq` become prohibitively expensive. To circumvent this difficulty, we tabulate :math:`y = ax - 1/\tan{x}` and interpolate it on the tabulated grid. In fact, we will use the *inverse* interpolation: we interpolate the values of :math:`x` versus :math:`у`. This way, solving the original equation becomes simply an evaluation of the interpolated function at zero :math:`y` argument. To improve the interpolation accuracy we will use the knowledge of the derivatives of the tabulated function. We will use `BPoly.from_derivatives` to construct a cubic interpolant (equivalently, we could have used `CubicHermiteSpline`) .. plot:: import numpy as np import matplotlib.pyplot as plt from scipy.interpolate import BPoly def f(x, a): return a*x - 1/np.tan(x) xleft, xright = 0.2, np.pi/2 x = np.linspace(xleft, xright, 11) fig, ax = plt.subplots(1, 2, figsize=(12, 4)) for j, a in enumerate([3, 93]): y = f(x, a) dydx = a + 1./np.sin(x)**2 # d(ax - 1/tan(x)) / dx dxdy = 1 / dydx # dx/dy = 1 / (dy/dx) xdx = np.c_[x, dxdy] spl = BPoly.from_derivatives(y, xdx) # inverse interpolation yy = np.linspace(f(xleft, a), f(xright, a), 51) ax[j].plot(yy, spl(yy), '--') ax[j].plot(y, x, 'o') ax[j].set_xlabel(r'$y$') ax[j].set_ylabel(r'$x$') ax[j].set_title(rf'$a = {a}$') ax[j].plot(0, spl(0), 'o', ms=12) ax[j].text(0.1, 0.85, fr'$x_0 = {spl(0):.3f}$', transform=ax[j].transAxes, fontsize=18) ax[j].grid(True) plt.tight_layout() plt.show() Note that for :math:`a=3`, ``spl(0)`` agrees with the ``brentq`` call above, while for :math:`a = 93`, the difference is substantial. The reason the procedure starts failing for large :math:`a` is that the straight line :math:`y = ax` tends towards the vertical axis, and the root of the original equation tends towards :math:`x=0`. Since we tabulated the original function at a finite grid, ``spl(0)`` involves extrapolation for too-large values of :math:`a`. Relying on extrapolation is prone to losing accuracy and is best avoided. Use the known asymptotics ~~~~~~~~~~~~~~~~~~~~~~~~~ Looking at the original equation, we note that for :math:`x\to 0`, :math:`\tan(x) = x + O(x^3)`, and the original equation becomes .. math:: ax = 1/x \;, so that :math:`x_0 \approx 1/\sqrt{a}` for :math:`a \gg 1`. We will use this to cook up a class which switches from interpolation to using this known asymptotic behavior for out-of-range data. A bare-bones implementation may look like this .. code:: python class RootWithAsymptotics: def __init__(self, a): # construct the interpolant xleft, xright = 0.2, np.pi/2 x = np.linspace(xleft, xright, 11) y = f(x, a) dydx = a + 1./np.sin(x)**2 # d(ax - 1/tan(x)) / dx dxdy = 1 / dydx # dx/dy = 1 / (dy/dx) # inverse interpolation self.spl = BPoly.from_derivatives(y, np.c_[x, dxdy]) self.a = a def root(self): out = self.spl(0) asympt = 1./np.sqrt(self.a) return np.where(spl.x.min() < asympt, out, asympt) And then >>> r = RootWithAsymptotics(93) >>> r.root() array(0.10369517) which differs from the extrapolated result and agrees with the ``brentq`` call. Note that this implementation is intentionally pared down. From the API perspective, you may want to instead implement the ``__call__`` method so that the full dependence of ``x`` on ``y`` is available. From the numerical perspective, more work is needed to make sure that the switch between interpolation and asymptotics occurs deep enough into the asymptotic regime, so that the resulting function is smooth enough at the switch-over point. Also in this example we artificially limited the problem to only consider a single periodicity interval of the ``tan`` function, and only dealt with :math:`a > 0`. For negative values of :math:`a`, we would need to implement the other asymptotics, for :math:`x\to \pi`. However the basic idea is the same. .. _tutorial-extrapolation-CT_NN: Extrapolation in ``D > 1`` ========================== The basic idea of implementing extrapolation manually in a wrapper class or function can be easily generalized to higher dimensions. As an example, we consider a C1-smooth interpolation problem of 2D data using `CloughTocher2DInterpolator`. By default, it fills the out of bounds values with ``nan``\ s, and we want to instead use for each query point the value of its nearest neighbor. Since `CloughTocher2DInterpolator` accepts either 2D data or a Delaunay triangulation of the data points, the efficient way of finding nearest neighbors of query points would be to construct the triangulation (using `scipy.spatial` tools) and use it to find nearest neighbors on the convex hull of the data. We will instead use a simpler, naive method and rely on looping over the whole dataset using NumPy broadcasting. .. plot:: import numpy as np import matplotlib.pyplot as plt from scipy.interpolate import CloughTocher2DInterpolator as CT def my_CT(xy, z): """CT interpolator + nearest-neighbor extrapolation. Parameters ---------- xy : ndarray, shape (npoints, ndim) Coordinates of data points z : ndarray, shape (npoints) Values at data points Returns ------- func : callable A callable object which mirrors the CT behavior, with an additional neareast-neighbor extrapolation outside of the data range. """ x = xy[:, 0] y = xy[:, 1] f = CT(xy, z) # this inner function will be returned to a user def new_f(xx, yy): # evaluate the CT interpolator. Out-of-bounds values are nan. zz = f(xx, yy) nans = np.isnan(zz) if nans.any(): # for each nan point, find its nearest neighbor inds = np.argmin( (x[:, None] - xx[nans])**2 + (y[:, None] - yy[nans])**2 , axis=0) # ... and use its value zz[nans] = z[inds] return zz return new_f # Now illustrate the difference between the original ``CT`` interpolant # and ``my_CT`` on a small example: x = np.array([1, 1, 1, 2, 2, 2, 4, 4, 4]) y = np.array([1, 2, 3, 1, 2, 3, 1, 2, 3]) z = np.array([0, 7, 8, 3, 4, 7, 1, 3, 4]) xy = np.c_[x, y] lut = CT(xy, z) lut2 = my_CT(xy, z) X = np.linspace(min(x) - 0.5, max(x) + 0.5, 71) Y = np.linspace(min(y) - 0.5, max(y) + 0.5, 71) X, Y = np.meshgrid(X, Y) fig = plt.figure() ax = fig.add_subplot(projection='3d') ax.plot_wireframe(X, Y, lut(X, Y), label='CT') ax.plot_wireframe(X, Y, lut2(X, Y), color='m', cstride=10, rstride=10, alpha=0.7, label='CT + n.n.') ax.scatter(x, y, z, 'o', color='k', s=48, label='data') ax.legend() plt.tight_layout()