Merge pull request #3267 from ev-br/pchip_

ENH: interpolate: reimplement pchip interpolator

Convert pchip interpolator into a wrapper around BPoly. The interpolation
algorithm is the same as previously, but evaluation is done via the new
piecewise polynomial functionality.

scipy/interpolate/__init__.py

@@ -163,6 +163,8 @@
 
 from .polyint import *
 
+from ._monotone import *
+
 from .ndgriddata import *
 
 __all__ = [s for s in dir() if not s.startswith('_')]

scipy/interpolate/_monotone.py

@@ -0,0 +1,224 @@
+from __future__ import division, print_function, absolute_import
+
+import numpy as np
+
+from . import BPoly
+
+__all__ = ["PchipInterpolator", "pchip_interpolate", "pchip"]
+
+class PchipInterpolator(object):
+    """PCHIP 1-d monotonic cubic interpolation
+
+    x and y are arrays of values used to approximate some function f,
+    with ``y = f(x)``.  The interpolant uses monotonic cubic splines
+    to find the value of new points. (PCHIP stands for Piecewise Cubic
+    Hermite Interpolating Polynomial).
+
+    Parameters
+    ----------
+    x : ndarray
+        A 1-D array of monotonically increasing real values.  `x` cannot
+        include duplicate values (otherwise f is overspecified)
+    y : ndarray
+        A 1-D array of real values.  `y`'s length along the interpolation
+        axis must be equal to the length of `x`. If N-D array, use axis 
+        parameter to select correct axis.
+    axis : int, optional
+        Axis in the y array corresponding to the x-coordinate values.
+    extrapolate : bool, optional
+        Whether to extrapolate to ouf-of-bounds points based on first
+        and last intervals, or to return NaNs.
+
+    Methods
+    -------
+    __call__
+    derivative
+
+    Notes
+    -----
+    The first derivatives are guaranteed to be continuous, but the second
+    derivatives may jump at x_k. 
+
+    Preserves monotonicity in the interpolation data and does not overshoot
+    if the data is not smooth.
+
+    Determines the derivatives at the points x_k, d_k, by using PCHIP algorithm:
+      
+    Let m_k be the slope of the kth segment (between k and k+1)
+    If m_k=0 or m_{k-1}=0 or sgn(m_k) != sgn(m_{k-1}) then d_k == 0
+    else use weighted harmonic mean:
+
+       w_1 = 2h_k + h_{k-1}, w_2 = h_k + 2h_{k-1}
+       1/d_k = 1/(w_1 + w_2)*(w_1 / m_k + w_2 / m_{k-1})
+
+    where h_k is the spacing between x_k and x_{k+1}.
+
+    """
+    def __init__(self, x, y, axis=0, extrapolate=None):
+        x = np.asarray(x)
+        y = np.asarray(y)
+
+        axis = axis % y.ndim
+
+        xp = x.reshape((x.shape[0],) + (1,)*(y.ndim-1))
+        yp = np.rollaxis(y, axis)
+
+        dk = self._find_derivatives(xp, yp)
+        data = np.hstack((yp[:, None, ...], dk[:, None, ...]))
+
+        self._bpoly = BPoly.from_derivatives(x, data, orders=None,
+                extrapolate=extrapolate)
+        self.axis = axis
+
+    def __call__(self, x, nu=0, extrapolate=None):
+        """
+        Evaluate the PCHIP interpolant or its derivative.
+
+        Parameters
+        ----------
+        x : array-like
+            Points to evaluate the interpolant at.
+        nu : int, optional
+            Order of derivative to evaluate. Must be non-negative.
+        extrapolate : bool, optional
+            Whether to extrapolate to ouf-of-bounds points based on first
+            and last intervals, or to return NaNs.
+
+        Returns
+        -------
+        y : array-like
+            Interpolated values. Shape is determined by replacing
+            the interpolation axis in the original array with the shape of x.
+
+        """
+        out = self._bpoly(x, nu, extrapolate)
+        return self._reshaper(x, out)
+
+    def derivative(self, der=1):
+        """
+        Construct a piecewise polynomial representing the derivative.
+
+        Parameters
+        ----------
+        der : int, optional
+            Order of derivative to evaluate. (Default: 1)
+            If negative, the antiderivative is returned.
+
+        Returns
+        ------- 
+        Piecewise polynomial of order k2 = k - der representing the derivative
+        of this polynomial.
+
+        """
+        t = object.__new__(self.__class__)
+        t.axis = self.axis
+        t._bpoly = self._bpoly.derivative(der)
+        return t
+
+    def roots(self):
+        """
+        Return the roots of the interpolated function.
+        """
+        return (PPoly.from_bernstein_basis(self._bpoly)).roots()
+
+    def _reshaper(self, x, out):
+        x = np.asarray(x)
+        l = x.ndim
+        transp = (tuple(range(l, l+self.axis)) + tuple(range(l)) +
+                tuple(range(l+self.axis, out.ndim)))
+        return out.transpose(transp)
+
+    @staticmethod
+    def _edge_case(m0, d1, out):
+        m0 = np.atleast_1d(m0)
+        d1 = np.atleast_1d(d1)
+        mask = (d1 != 0) & (m0 != 0)
+        out[mask] = 1.0/(1.0/m0[mask]+1.0/d1[mask])
+
+    @staticmethod
+    def _find_derivatives(x, y):
+        # Determine the derivatives at the points y_k, d_k, by using
+        #  PCHIP algorithm is:
+        # We choose the derivatives at the point x_k by
+        # Let m_k be the slope of the kth segment (between k and k+1)
+        # If m_k=0 or m_{k-1}=0 or sgn(m_k) != sgn(m_{k-1}) then d_k == 0
+        # else use weighted harmonic mean:
+        #   w_1 = 2h_k + h_{k-1}, w_2 = h_k + 2h_{k-1}
+        #   1/d_k = 1/(w_1 + w_2)*(w_1 / m_k + w_2 / m_{k-1})
+        #   where h_k is the spacing between x_k and x_{k+1}
+        y_shape = y.shape
+        if y.ndim == 1:
+            # So that _edge_case doesn't end up assigning to scalars
+            x = x[:,None]
+            y = y[:,None]
+
+        hk = x[1:] - x[:-1]
+        mk = (y[1:] - y[:-1]) / hk
+        smk = np.sign(mk)
+        condition = ((smk[1:] != smk[:-1]) | (mk[1:] == 0) | (mk[:-1] == 0))
+
+        w1 = 2*hk[1:] + hk[:-1]
+        w2 = hk[1:] + 2*hk[:-1]
+        # values where division by zero occurs will be excluded
+        # by 'condition' afterwards
+        with np.errstate(divide='ignore'):
+            whmean = 1.0/(w1+w2)*(w1/mk[1:] + w2/mk[:-1])
+        dk = np.zeros_like(y)
+        dk[1:-1][condition] = 0.0
+        dk[1:-1][~condition] = 1.0/whmean[~condition]
+
+        # For end-points choose d_0 so that 1/d_0 = 1/m_0 + 1/d_1 unless
+        #  one of d_1 or m_0 is 0, then choose d_0 = 0
+        PchipInterpolator._edge_case(mk[0],dk[1], dk[0])
+        PchipInterpolator._edge_case(mk[-1],dk[-2], dk[-1])
+
+        return dk.reshape(y_shape)
+
+
+def pchip_interpolate(xi, yi, x, der=0, axis=0):
+    """
+    Convenience function for pchip interpolation. 
+    xi and yi are arrays of values used to approximate some function f,
+    with ``yi = f(xi)``.  The interpolant uses monotonic cubic splines
+    to find the value of new points x and the derivatives there.
+
+    See `PchipInterpolator` for details.
+
+    Parameters
+    ----------
+    xi : array_like
+        A sorted list of x-coordinates, of length N.
+    yi :  array_like
+        A 1-D array of real values.  `yi`'s length along the interpolation
+        axis must be equal to the length of `xi`. If N-D array, use axis 
+        parameter to select correct axis.
+    x : scalar or array_like
+        Of length M.
+    der : integer or list
+        How many derivatives to extract; None for all potentially
+        nonzero derivatives (that is a number equal to the number
+        of points), or a list of derivatives to extract. This number
+        includes the function value as 0th derivative.
+    axis : int, optional
+        Axis in the yi array corresponding to the x-coordinate values.
+
+    See Also
+    --------
+    PchipInterpolator
+
+    Returns
+    -------
+    y : scalar or array_like
+        The result, of length R or length M or M by R,
+
+    """
+    P = PchipInterpolator(xi, yi, axis=axis)
+    if der == 0:
+        return P(x)
+    elif _isscalar(der):
+        return P(x, der=der)
+    else:
+        return [P(x, nu) for nu in der]
+
+# Backwards compatibility
+pchip = PchipInterpolator

scipy/interpolate/polyint.py

@@ -8,8 +8,7 @@
 
 __all__ = ["KroghInterpolator", "krogh_interpolate", "BarycentricInterpolator",
            "barycentric_interpolate", "PiecewisePolynomial",
-           "piecewise_polynomial_interpolate", "approximate_taylor_polynomial",
-           "PchipInterpolator", "pchip_interpolate", "pchip"]
+           "piecewise_polynomial_interpolate", "approximate_taylor_polynomial"]
 
 
 def _isscalar(x):
@@ -923,145 +922,3 @@ def piecewise_polynomial_interpolate(xi,yi,x,orders=None,der=0,axis=0):
         return P.derivative(x,der=der)
     else:
         return P.derivatives(x,der=np.amax(der)+1)[der]
-
-
-class PchipInterpolator(PiecewisePolynomial):
-    """PCHIP 1-d monotonic cubic interpolation
-
-    x and y are arrays of values used to approximate some function f,
-    with ``y = f(x)``.  The interpolant uses monotonic cubic splines
-    to find the value of new points.
-
-    Parameters
-    ----------
-    x : ndarray
-        A 1-D array of monotonically increasing real values.  `x` cannot
-        include duplicate values (otherwise f is overspecified)
-    y : ndarray
-        A 1-D array of real values.  `y`'s length along the interpolation
-        axis must be equal to the length of `x`. If N-D array, use axis 
-        parameter to select correct axis.
-    axis : int, optional
-        Axis in the y array corresponding to the x-coordinate values.
-
-    Notes
-    -----
-    Assumes x is sorted in monotonic order (e.g. ``x[1] > x[0]``).
-
-    """
-    def __init__(self, x, y, axis=0):
-        x = np.asarray(x)
-        y = np.asarray(y)
-
-        axis = axis % y.ndim
-
-        xp = x.reshape((x.shape[0],) + (1,)*(y.ndim-1))
-        yp = np.rollaxis(y, axis)
-
-        data = np.empty((yp.shape[0], 2) + yp.shape[1:], y.dtype)
-        data[:,0] = yp
-        data[:,1] = PchipInterpolator._find_derivatives(xp, yp)
-
-        s = list(range(2, y.ndim + 1))
-        s.insert(axis, 1)
-        s.insert(axis, 0)
-        data = data.transpose(s)
-
-        PiecewisePolynomial.__init__(self, x, data, orders=3, direction=None,
-                                     axis=axis)
-
-    @staticmethod
-    def _edge_case(m0, d1, out):
-        m0 = np.atleast_1d(m0)
-        d1 = np.atleast_1d(d1)
-        mask = (d1 != 0) & (m0 != 0)
-        out[mask] = 1.0/(1.0/m0[mask]+1.0/d1[mask])
-
-    @staticmethod
-    def _find_derivatives(x, y):
-        # Determine the derivatives at the points y_k, d_k, by using
-        #  PCHIP algorithm is:
-        # We choose the derivatives at the point x_k by
-        # Let m_k be the slope of the kth segment (between k and k+1)
-        # If m_k=0 or m_{k-1}=0 or sgn(m_k) != sgn(m_{k-1}) then d_k == 0
-        # else use weighted harmonic mean:
-        #   w_1 = 2h_k + h_{k-1}, w_2 = h_k + 2h_{k-1}
-        #   1/d_k = 1/(w_1 + w_2)*(w_1 / m_k + w_2 / m_{k-1})
-        #   where h_k is the spacing between x_k and x_{k+1}
-
-        y_shape = y.shape
-        if y.ndim == 1:
-            # So that _edge_case doesn't end up assigning to scalars
-            x = x[:,None]
-            y = y[:,None]
-
-        hk = x[1:] - x[:-1]
-        mk = (y[1:] - y[:-1]) / hk
-        smk = np.sign(mk)
-        condition = ((smk[1:] != smk[:-1]) | (mk[1:] == 0) | (mk[:-1] == 0))
-
-        w1 = 2*hk[1:] + hk[:-1]
-        w2 = hk[1:] + 2*hk[:-1]
-        # values where division by zero occurs will be excluded
-        # by 'condition' afterwards
-        with np.errstate(divide='ignore'):
-            whmean = 1.0/(w1+w2)*(w1/mk[1:] + w2/mk[:-1])
-        dk = np.zeros_like(y)
-        dk[1:-1][condition] = 0.0
-        dk[1:-1][~condition] = 1.0/whmean[~condition]
-
-        # For end-points choose d_0 so that 1/d_0 = 1/m_0 + 1/d_1 unless
-        #  one of d_1 or m_0 is 0, then choose d_0 = 0
-        PchipInterpolator._edge_case(mk[0],dk[1], dk[0])
-        PchipInterpolator._edge_case(mk[-1],dk[-2], dk[-1])
-
-        return dk.reshape(y_shape)
-
-
-def pchip_interpolate(xi, yi, x, der=0, axis=0):
-    """
-    Convenience function for pchip interpolation. 
-    xi and yi are arrays of values used to approximate some function f,
-    with ``yi = f(xi)``.  The interpolant uses monotonic cubic splines
-    to find the value of new points x and the derivatives there.
-
-    See `PchipInterpolator` for details.
-
-    Parameters
-    ----------
-    xi : array_like
-        A sorted list of x-coordinates, of length N.
-    yi :  array_like
-        A 1-D array of real values.  `yi`'s length along the interpolation
-        axis must be equal to the length of `xi`. If N-D array, use axis 
-        parameter to select correct axis.
-    x : scalar or array_like
-        Of length M.
-    der : integer or list
-        How many derivatives to extract; None for all potentially
-        nonzero derivatives (that is a number equal to the number
-        of points), or a list of derivatives to extract. This number
-        includes the function value as 0th derivative.
-    axis : int, optional
-        Axis in the yi array corresponding to the x-coordinate values.
-
-    See Also
-    --------
-    PchipInterpolator
-
-    Returns
-    -------
-    y : scalar or array_like
-        The result, of length R or length M or M by R,
-
-    """
-    P = PchipInterpolator(xi, yi, axis=axis)
-    if der == 0:
-        return P(x)
-    elif _isscalar(der):
-        return P.derivative(x,der=der)
-    else:
-        return P.derivatives(x,der=np.amax(der)+1)[der]
-
-# Backwards compatibility
-pchip = PchipInterpolator