Merge pull request scipy#5351 from ev-br/pchip_slopes

ENH: interpolate: refactor the PCHIP interpolant's algorithm

scipy/interpolate/_monotone.py

@@ -12,10 +12,10 @@
 
 
 class PchipInterpolator(BPoly):
-    """PCHIP 1-d monotonic cubic interpolation
+    r"""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
+    `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).
 
@@ -25,8 +25,8 @@ class PchipInterpolator(BPoly):
         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
+        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.
@@ -39,30 +39,46 @@ class PchipInterpolator(BPoly):
     __call__
     derivative
     antiderivative
+    roots
 
     See Also
     --------
     Akima1DInterpolator
 
     Notes
     -----
+    The interpolator preserves monotonicity in the interpolation data and does
+    not overshoot if the data is not smooth.
+
     The first derivatives are guaranteed to be continuous, but the second
-    derivatives may jump at x_k.
+    derivatives may jump at :math:`x_k`.
+
+    Determines the derivatives at the points :math:`x_k`, :math:`f'_k`,
+    by using PCHIP algorithm [1]_.
+
+    Let :math:`h_k = x_{k+1} - x_k`, and  :math:`d_k = (y_{k+1} - y_k) / h_k`
+    are the slopes at internal points :math:`x_k`.
+    If the signs of :math:`d_k` and :math:`d_{k-1}` are different or either of
+    them equals zero, then :math:`f'_k = 0`. Otherwise, it is given by the
+    weighted harmonic mean
 
-    Preserves monotonicity in the interpolation data and does not overshoot
-    if the data is not smooth.
+    .. math::
 
-    Determines the derivatives at the points x_k, d_k, by using PCHIP
-    algorithm:
+        \frac{w_1 + w_2}{f'_k} = \frac{w_1}{d_{k-1}} + \frac{w_2}{d_k}
 
-    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:
+    where :math:`w_1 = 2 h_k + h_{k-1}` and :math:`w_2 = h_k + 2 h_{k-1}`.
 
-       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})
+    The end slopes are set using a one-sided scheme [2]_.
+
+
+    References
+    ----------
+    .. [1] F. N. Fritsch and R. E. Carlson, Monotone Piecewise Cubic Interpolation,
+           SIAM J. Numer. Anal., 17(2), 238 (1980).
+           DOI:10.1137/0717021
+    .. [2] see, e.g., C. Moler, Numerical Computing with Matlab, 2004.
+           DOI: http://dx.doi.org/10.1137/1.9780898717952
 
-    where h_k is the spacing between x_k and x_{k+1}.
 
     """
     def __init__(self, x, y, axis=0, extrapolate=None):
@@ -89,11 +105,19 @@ def roots(self):
         return (PPoly.from_bernstein_basis(self._bpoly)).roots()
 
     @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])
+    def _edge_case(h0, h1, m0, m1):
+        # one-sided three-point estimate for the derivative
+        d = ((2*h0 + h1)*m0 - h0*m1) / (h0 + h1)
+
+        # try to preserve shape
+        mask = np.sign(d) != np.sign(m0)
+        mask2 = (np.sign(m0) != np.sign(m1)) & (np.abs(d) > 3.*np.abs(m0))
+        mmm = (~mask) & mask2
+
+        d[mask] = 0.
+        d[mmm] = 3.*m0[mmm]
+
+        return d
 
     @staticmethod
     def _find_derivatives(x, y):
@@ -115,22 +139,24 @@ def _find_derivatives(x, y):
         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))
+        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])
+            whmean = (w1/mk[:-1] + w2/mk[1:]) / (w1 + w2)
+
         dk = np.zeros_like(y)
         dk[1:-1][condition] = 0.0
-        dk[1:-1][~condition] = 1.0/whmean[~condition]
+        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])
+        # special case endpoints, as suggested in 
+        # Cleve Moler, Numerical Computing with MATLAB, Chap 3.4
+        dk[0] = PchipInterpolator._edge_case(hk[0], hk[1], mk[0], mk[1])
+        dk[-1] = PchipInterpolator._edge_case(hk[-1], hk[-2], mk[-1], mk[-2])
 
         return dk.reshape(y_shape)
 

scipy/interpolate/tests/test_polyint.py

@@ -492,5 +492,64 @@ def test_cast(self):
 
         assert_allclose(curve, curve1, atol=1e-14, rtol=1e-14)
 
+    def test_nag(self):
+        # Example from NAG C implementation,
+        # http://nag.com/numeric/cl/nagdoc_cl25/html/e01/e01bec.html
+        # suggested in gh-5326 as a smoke test for the way the derivatives
+        # are computed (see also gh-3453)
+        from scipy._lib.six import StringIO
+        dataStr = '''
+          7.99   0.00000E+0
+          8.09   0.27643E-4
+          8.19   0.43750E-1
+          8.70   0.16918E+0
+          9.20   0.46943E+0
+         10.00   0.94374E+0
+         12.00   0.99864E+0
+         15.00   0.99992E+0
+         20.00   0.99999E+0
+        '''
+        data = np.loadtxt(StringIO(dataStr))
+        pch = pchip(data[:,0], data[:,1])
+
+        resultStr = '''
+           7.9900       0.0000
+           9.1910       0.4640
+          10.3920       0.9645
+          11.5930       0.9965
+          12.7940       0.9992
+          13.9950       0.9998
+          15.1960       0.9999
+          16.3970       1.0000
+          17.5980       1.0000
+          18.7990       1.0000
+          20.0000       1.0000
+        '''
+        result = np.loadtxt(StringIO(resultStr))
+        assert_allclose(result[:,1], pch(result[:,0]), rtol=0., atol=5e-5)
+
+    def test_endslopes(self):
+        # this is a smoke test for gh-3453: PCHIP interpolator should not
+        # set edge slopes to zero if the data do not suggest zero edge derivatives
+        x = np.array([0.0, 0.1, 0.25, 0.35])
+        y1 = np.array([279.35, 0.5e3, 1.0e3, 2.5e3])
+        y2 = np.array([279.35, 2.5e3, 1.50e3, 1.0e3])
+        for pp in (pchip(x, y1), pchip(x, y2)):
+            for t in (x[0], x[-1]):
+                assert_(pp(t, 1) != 0)
+
+    def test_all_zeros(self):
+        x = np.arange(10)
+        y = np.zeros_like(x)
+
+        # this should work and not generate any warnings
+        with warnings.catch_warnings():
+            warnings.filterwarnings('error')
+            pch = pchip(x, y)
+
+        xx = np.linspace(0, 9, 101)
+        assert_equal(pch(xx), 0.)
+
+
 if __name__ == '__main__':
     run_module_suite()