Merge pull request #3320 from ev-br/poly_

ENH: interpolate: improve BPoly API and functionality

- remove unneeded attributes (axis and direction).
- implement in-place evaluation of derivatives.
- add tests for complex coefficients.
- deprecate PiecewisePolynomial

doc/release/0.14.0-notes.rst

@@ -124,6 +124,12 @@ Probability calculation aliases ``zprob``, ``fprob`` and ``ksprob`` are
 deprecated. Use instead the ``sf`` methods of the corresponding distributions
 or the ``special`` functions directly.
 
+``scipy.interpolate``
+---------------------
+
+``PiecewisePolynomial`` class is deprecated.
+
+
 Backwards incompatible changes
 ==============================
 

scipy/interpolate/_ppoly.pyx

@@ -6,6 +6,7 @@ local power basis.
 
 from .polyint import _Interpolator1D
 import numpy as np
+cimport numpy as cnp
 
 cimport cython
 
@@ -26,7 +27,6 @@ cdef extern:
                 double *vr, int *ldvr, double *work, int *lwork,
                 int *info)
 
-
 #------------------------------------------------------------------------------
 # Piecewise power basis polynomials
 #------------------------------------------------------------------------------
@@ -725,15 +725,17 @@ def _croots_poly1(double[:,:,::1] c, double_complex[:,:,::1] w):
 @cython.wraparound(False)
 @cython.boundscheck(False)
 @cython.cdivision(True)
-cdef double_or_complex evaluate_bpoly1(double_or_complex s, double_or_complex[:,:,::1] c, int ci, int cj) nogil:
+cdef double_or_complex evaluate_bpoly1(double_or_complex s,
+                                       double_or_complex[:,:,::1] c,
+                                       int ci, int cj) nogil:
     """
     Evaluate polynomial in the Berstein basis in a single interval.
 
-    A Berstein polynomial is defined as ..math::
+    A Berstein polynomial is defined as
 
-        b_{j, k} = comb(k, j) x^{j} (1-x)^{k-j}
+        .. math:: b_{j, k} = comb(k, j) x^{j} (1-x)^{k-j}
 
-    with ``0 <= x <= 1``
+    with ``0 <= x <= 1``.
 
     Parameters
     ----------
@@ -770,6 +772,65 @@ cdef double_or_complex evaluate_bpoly1(double_or_complex s, double_or_complex[:,
 
     return res
 
+
+@cython.wraparound(False)
+@cython.boundscheck(False)
+@cython.cdivision(True)
+cdef double_or_complex evaluate_bpoly1_deriv(double_or_complex s,
+                                             double_or_complex[:,:,::1] c,
+                                             int ci, int cj,
+                                             int nu,
+                                             double_or_complex[:,:,::1] wrk) nogil:
+    """
+    Evaluate the derivative of a polynomial in the Berstein basis 
+    in a single interval.
+
+    A Berstein polynomial is defined as
+
+        .. math:: b_{j, k} = comb(k, j) x^{j} (1-x)^{k-j}
+
+    with ``0 <= x <= 1``.
+
+    The algorithm is detailed in BPoly._construct_from_derivatives.
+
+    Parameters
+    ----------
+    s : double
+        Polynomial x-value
+    c : double[:,:,:]
+        Polynomial coefficients. c[:,ci,cj] will be used
+    ci, cj : int
+        Which of the coefs to use
+    nu : int
+        Order of the derivative to evaluate. Assumed strictly positive
+        (no checks are made).
+    wrk : double[:,:,::1]
+        A work array, shape (c.shape[0]-nu, 1, 1).
+
+    """
+    cdef int k, j, a
+    cdef double_or_complex res, term
+    cdef double comb, poch
+
+    k = c.shape[0] - 1  # polynomial order
+
+    if nu == 0:
+        res = evaluate_bpoly1(s, c, ci, cj)
+    else:
+        poch = 1.
+        for a in range(nu):
+            poch *= k - a
+
+        term = 0.
+        for a in range(k - nu + 1):
+            term, comb = 0., 1.
+            for j in range(nu+1):
+                term += c[j+a, ci, cj] * (-1)**(j+nu) * comb
+                comb *= 1. * (nu-j) / (j+1)
+            wrk[a, 0, 0] = term * poch
+        res = evaluate_bpoly1(s, wrk, 0, 0)
+    return res
+
 #
 # Evaluation; only differs from _ppoly by evaluate_poly1 -> evaluate_bpoly1
 #
@@ -779,9 +840,10 @@ cdef double_or_complex evaluate_bpoly1(double_or_complex s, double_or_complex[:,
 def evaluate_bernstein(double_or_complex[:,:,::1] c,
              double[::1] x,
              double[::1] xp,
-             int dx,
+             int nu,
              int extrapolate,
-             double_or_complex[:,::1] out):
+             double_or_complex[:,::1] out,
+             cnp.dtype dt):
     """
     Evaluate a piecewise polynomial in the Bernstein basis.
 
@@ -795,7 +857,7 @@ def evaluate_bernstein(double_or_complex[:,:,::1] c,
         Breakpoints of polynomials
     xp : ndarray, shape (r,)
         Points to evaluate the piecewise polynomial at.
-    dx : int
+    nu : int
         Order of derivative to evaluate.  The derivative is evaluated
         piecewise and may have discontinuities.
     extrapolate : int, optional
@@ -810,11 +872,12 @@ def evaluate_bernstein(double_or_complex[:,:,::1] c,
     cdef int ip, jp
     cdef int interval
     cdef double xval
-    cdef double_or_complex s
+    cdef double_or_complex s, ds, ds_nu
+    cdef double_or_complex[:,:,::1] wrk
 
     # check derivative order
-    if dx != 0:
-        raise NotImplementedError("Cannot do derivatives in the B-basis yet.")
+    if nu < 0:
+        raise NotImplementedError("Cannot do antiderivatives in the B-basis yet.")
 
     # shape checks
     if out.shape[0] != xp.shape[0]:
@@ -824,6 +887,9 @@ def evaluate_bernstein(double_or_complex[:,:,::1] c,
     if c.shape[1] != x.shape[0] - 1:
         raise ValueError("x and c have incompatible shapes")
 
+    if nu > 0:
+        wrk = np.empty((c.shape[0]-nu, 1, 1), dtype=dt)
+
     # evaluate
     interval = 0
 
@@ -841,7 +907,12 @@ def evaluate_bernstein(double_or_complex[:,:,::1] c,
             interval = i
 
         # Evaluate the local polynomial(s)
+        ds = x[interval+1] - x[interval]
+        ds_nu = ds**nu
         for jp in range(c.shape[2]):
-            s = (xval - x[interval]) / (x[interval+1] - x[interval]) 
-            out[ip, jp] = evaluate_bpoly1(s, c, interval, jp)
-
+            s = (xval - x[interval]) / ds
+            if nu == 0:
+                out[ip, jp] = evaluate_bpoly1(s, c, interval, jp)
+            else:
+                out[ip, jp] = evaluate_bpoly1_deriv(s, c, interval, jp,
+                        nu, wrk) / ds_nu

scipy/interpolate/interpolate.py

@@ -1054,7 +1054,7 @@ class BPoly(_PPolyBase):
     def _evaluate(self, x, nu, extrapolate, out):
         _ppoly.evaluate_bernstein(
             self.c.reshape(self.c.shape[0], self.c.shape[1], -1),
-            self.x, x, nu, bool(extrapolate), out)
+            self.x, x, nu, bool(extrapolate), out, self.c.dtype)
 
     def derivative(self, nu=1):
         """
@@ -1148,8 +1148,7 @@ def from_power_basis(cls, pp, extrapolate=None):
         return cls.construct_fast(c, pp.x, extrapolate)
 
     @classmethod
-    def from_derivatives(cls, xi, yi, orders=None, direction=None,
-                         axis=0, extrapolate=None):
+    def from_derivatives(cls, xi, yi, orders=None, extrapolate=None):
         """Construct a piecewise polynomial in the Bernstein basis,
         compatible with the specified values and derivatives at breakpoints.
 
@@ -1162,8 +1161,6 @@ def from_derivatives(cls, xi, yi, orders=None, direction=None,
         orders : None or int or array_like of ints. Default: None.
             Specifies the degree of local polynomials. If not None, some
             derivatives are ignored.
-        axis : int, optional
-            Interpolation axis, default is 0.
         extrapolate : bool, optional
             Whether to extrapolate to ouf-of-bounds points based on first
             and last intervals, or to return NaNs. Default: True.
@@ -1190,12 +1187,12 @@ def from_derivatives(cls, xi, yi, orders=None, direction=None,
         Examples
         --------
 
-        >>> BPoly.from_derivatives([[1, 2], [3, 4]], [0, 1])
+        >>> BPoly.from_derivatives([0, 1], [[1, 2], [3, 4]])
 
         Creates a polynomial `f(x)` of degree 3, defined on `[0, 1]`
         such that `f(0) = 1, df/dx(0) = 2, f(1) = 3, df/dx(1) = 4`
 
-        >>> BPoly.from_derivatives([[0, 1], [0], [2], [0, 1, 2])
+        >>> BPoly.from_derivatives([0, 1, 2], [[0, 1], [0], [2]])
 
         Creates a piecewise polynomial `f(x)`, such that
         `f(0) = f(1) = 0`, `f(2) = 2`, and `df/dx(0) = 1`.
@@ -1212,11 +1209,6 @@ def from_derivatives(cls, xi, yi, orders=None, direction=None,
         So that f'(1-0) = -1 and f'(1+0) = 2
 
         """
-        if axis != 0:
-            raise NotImplementedError
-        if direction is not None:
-            raise NotImplementedError
-
         xi = np.asarray(xi)
         if len(xi) != len(yi):
             raise ValueError("xi and yi need to have the same length")
@@ -1227,7 +1219,11 @@ def from_derivatives(cls, xi, yi, orders=None, direction=None,
         m = len(xi) - 1
 
         # global poly order is k-1, local orders are <=k and can vary
-        k = max(len(yi[i]) + len(yi[i+1]) for i in range(m))
+        try:
+            k = max(len(yi[i]) + len(yi[i+1]) for i in range(m))
+        except TypeError:
+            raise ValueError("Using a 1D array for y? Please .reshape(-1, 1).")
+
         if orders is None:
             orders = [None] * m
         else:
@@ -1256,7 +1252,7 @@ def from_derivatives(cls, xi, yi, orders=None, direction=None,
                     raise ValueError("`order` input incompatible with"
                             " length y1 or y2.")
 
-            b = BPoly._construct_from_derivatives(xi[i], xi[i+1],  y1[:n1], y2[:n2])
+            b = BPoly._construct_from_derivatives(xi[i], xi[i+1], y1[:n1], y2[:n2])
             if len(b) < k:
                 b = BPoly._raise_degree(b, k - len(b))
             c.append(b)
@@ -1271,10 +1267,10 @@ def _construct_from_derivatives(xa, xb, ya, yb):
 
         Return the coefficients of a polynomial in the Bernstein basis
         defined on `[xa, xb]` and having the values and derivatives at the
-        endpoints `xa` and `xb` as specified by `ya` and `yb`, respectively.
+        endpoints ``xa`` and ``xb`` as specified by ``ya`` and ``yb``.
         The polynomial constructed is of the minimal possible degree, i.e.,
-        if the lengths of `ya` and `yb` are `na` and `nb`, the degree of the
-        polynomial is `na + nb - 1`.
+        if the lengths of ``ya`` and ``yb`` are ``na`` and ``nb``, the degree
+        of the polynomial is ``na + nb - 1``.
 
         Parameters
         ----------
@@ -1283,10 +1279,10 @@ def _construct_from_derivatives(xa, xb, ya, yb):
         xb : float
             Right-hand end point of the interval
         ya : array_like
-            Derivatives at `xa`. `ya[0]` is the value of the function, and
-            `ya[i]` for `i > 0` is the value of the `i`-th derivative.
+            Derivatives at ``xa``. ``ya[0]`` is the value of the function, and
+            ``ya[i]`` for ``i > 0`` is the value of the ``i``-th derivative.
         yb : array_like
-            Derivatives at `xb`.
+            Derivatives at ``xb``.
 
         Returns
         -------

scipy/interpolate/polyint.py

@@ -1,5 +1,7 @@
 from __future__ import division, print_function, absolute_import
 
+import warnings
+
 import numpy as np
 from scipy.misc import factorial
 
@@ -703,6 +705,10 @@ class PiecewisePolynomial(_Interpolator1DWithDerivatives):
     def __init__(self, xi, yi, orders=None, direction=None, axis=0):
         _Interpolator1DWithDerivatives.__init__(self, axis=axis)
 
+        warnings.warn('PiecewisePolynomial is deprecated in scipy 0.14. '
+                      'Use BPoly.from_derivatives instead.',
+                category=DeprecationWarning)
+
         if axis != 0:
             try:
                 yi = np.asarray(yi)

scipy/interpolate/tests/test_interpolate.py

@@ -425,6 +425,18 @@ def test_shape(self):
                 # constructing the object array here for old numpy)
                 assert_raises(ValueError, p, np.array([[0.1, 0.2], [0.4]]))
 
+    def test_complex_coef(self):
+        np.random.seed(12345)
+        x = np.sort(np.random.random(13))
+        c = np.random.random((8, 12)) * (1. + 0.3j)
+        c_re, c_im = c.real, c.imag
+        xp = np.random.random(5)
+        for cls in (PPoly, BPoly):
+            p, p_re, p_im = cls(c, x), cls(c_re, x), cls(c_im, x)
+            for nu in [0, 1, 2]:
+                assert_allclose(p(xp, nu).real, p_re(xp, nu))
+                assert_allclose(p(xp, nu).imag, p_im(xp, nu))
+
 
 class TestPolySubclassing(TestCase):
     class P(PPoly):
@@ -867,6 +879,12 @@ def test_derivative(self):
         assert_allclose(bp_der(0.4), -6*(0.6))
         assert_allclose(bp_der(1.7), 0.7)
 
+        # derivatives in-place
+        assert_allclose([bp(0.4, nu=1), bp(0.4, nu=2), bp(0.4, nu=3)],
+                        [-6*(1-0.4), 6., 0.])
+        assert_allclose([bp(1.7, nu=1), bp(1.7, nu=2), bp(1.7, nu=3)],
+                        [0.7, 1., 0])
+
     def test_derivative_ppoly(self):
         # make sure it's consistent w/ power basis
         np.random.seed(1234)
@@ -882,6 +900,17 @@ def test_derivative_ppoly(self):
             xp = np.linspace(x[0], x[-1], 21)
             assert_allclose(bp(xp), pp(xp))
 
+    def test_deriv_inplace(self):
+        np.random.seed(1234)
+        m, k = 5, 8   # number of intervals, order
+        x = np.sort(np.random.random(m))
+        c = np.random.random((k, m-1))
+        bp = BPoly(c, x)
+
+        xp = np.linspace(x[0], x[-1], 21)
+        for i in range(k):
+            assert_allclose(bp(xp, i), bp.derivative(i)(xp))
+
 
 class TestPolyConversions(TestCase):
     def test_bp_from_pp(self):