diff --git a/doc/release/0.14.0-notes.rst b/doc/release/0.14.0-notes.rst index 7f2ff726abb6..f47c6989d41e 100644 --- a/doc/release/0.14.0-notes.rst +++ b/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 ============================== diff --git a/scipy/interpolate/_ppoly.pyx b/scipy/interpolate/_ppoly.pyx index 85030b7648b2..84e145a3a32f 100644 --- a/scipy/interpolate/_ppoly.pyx +++ b/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 diff --git a/scipy/interpolate/interpolate.py b/scipy/interpolate/interpolate.py index a6f972dfd42e..760d8140fc89 100644 --- a/scipy/interpolate/interpolate.py +++ b/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 ------- diff --git a/scipy/interpolate/polyint.py b/scipy/interpolate/polyint.py index 21ba8db51876..06203c68be49 100644 --- a/scipy/interpolate/polyint.py +++ b/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) diff --git a/scipy/interpolate/tests/test_interpolate.py b/scipy/interpolate/tests/test_interpolate.py index de8913a7250a..d75b156beb53 100644 --- a/scipy/interpolate/tests/test_interpolate.py +++ b/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):