from numpy.testing import (assert_, assert_equal, assert_almost_equal, assert_array_almost_equal, assert_array_equal, assert_allclose, suppress_warnings) from pytest import raises as assert_raises import pytest from numpy import mgrid, pi, sin, ogrid, poly1d, linspace import numpy as np from scipy.interpolate import (interp1d, interp2d, lagrange, PPoly, BPoly, splrep, splev, splantider, splint, sproot, Akima1DInterpolator, NdPPoly, BSpline) from scipy.special import poch, gamma from scipy.interpolate import _ppoly from scipy._lib._gcutils import assert_deallocated, IS_PYPY from scipy.integrate import nquad from scipy.special import binom class TestInterp2D: def test_interp2d(self): y, x = mgrid[0:2:20j, 0:pi:21j] z = sin(x+0.5*y) with suppress_warnings() as sup: sup.filter(DeprecationWarning) II = interp2d(x, y, z) assert_almost_equal(II(1.0, 2.0), sin(2.0), decimal=2) v, u = ogrid[0:2:24j, 0:pi:25j] assert_almost_equal(II(u.ravel(), v.ravel()), sin(u+0.5*v), decimal=2) def test_interp2d_meshgrid_input(self): # Ticket #703 x = linspace(0, 2, 16) y = linspace(0, pi, 21) z = sin(x[None, :] + y[:, None]/2.) with suppress_warnings() as sup: sup.filter(DeprecationWarning) II = interp2d(x, y, z) assert_almost_equal(II(1.0, 2.0), sin(2.0), decimal=2) def test_interp2d_meshgrid_input_unsorted(self): np.random.seed(1234) x = linspace(0, 2, 16) y = linspace(0, pi, 21) z = sin(x[None, :] + y[:, None] / 2.) with suppress_warnings() as sup: sup.filter(DeprecationWarning) ip1 = interp2d(x.copy(), y.copy(), z, kind='cubic') np.random.shuffle(x) z = sin(x[None, :] + y[:, None]/2.) ip2 = interp2d(x.copy(), y.copy(), z, kind='cubic') np.random.shuffle(x) np.random.shuffle(y) z = sin(x[None, :] + y[:, None] / 2.) ip3 = interp2d(x, y, z, kind='cubic') x = linspace(0, 2, 31) y = linspace(0, pi, 30) assert_equal(ip1(x, y), ip2(x, y)) assert_equal(ip1(x, y), ip3(x, y)) def test_interp2d_eval_unsorted(self): y, x = mgrid[0:2:20j, 0:pi:21j] z = sin(x + 0.5*y) with suppress_warnings() as sup: sup.filter(DeprecationWarning) func = interp2d(x, y, z) xe = np.array([3, 4, 5]) ye = np.array([5.3, 7.1]) assert_allclose(func(xe, ye), func(xe, ye[::-1])) assert_raises(ValueError, func, xe, ye[::-1], 0, 0, True) def test_interp2d_linear(self): # Ticket #898 a = np.zeros([5, 5]) a[2, 2] = 1.0 x = y = np.arange(5) with suppress_warnings() as sup: sup.filter(DeprecationWarning) b = interp2d(x, y, a, 'linear') assert_almost_equal(b(2.0, 1.5), np.array([0.5]), decimal=2) assert_almost_equal(b(2.0, 2.5), np.array([0.5]), decimal=2) def test_interp2d_bounds(self): x = np.linspace(0, 1, 5) y = np.linspace(0, 2, 7) z = x[None, :]**2 + y[:, None] ix = np.linspace(-1, 3, 31) iy = np.linspace(-1, 3, 33) with suppress_warnings() as sup: sup.filter(DeprecationWarning) b = interp2d(x, y, z, bounds_error=True) assert_raises(ValueError, b, ix, iy) b = interp2d(x, y, z, fill_value=np.nan) iz = b(ix, iy) mx = (ix < 0) | (ix > 1) my = (iy < 0) | (iy > 2) assert_(np.isnan(iz[my, :]).all()) assert_(np.isnan(iz[:, mx]).all()) assert_(np.isfinite(iz[~my, :][:, ~mx]).all()) class TestInterp1D: def setup_method(self): self.x5 = np.arange(5.) self.x10 = np.arange(10.) self.y10 = np.arange(10.) self.x25 = self.x10.reshape((2,5)) self.x2 = np.arange(2.) self.y2 = np.arange(2.) self.x1 = np.array([0.]) self.y1 = np.array([0.]) self.y210 = np.arange(20.).reshape((2, 10)) self.y102 = np.arange(20.).reshape((10, 2)) self.y225 = np.arange(20.).reshape((2, 2, 5)) self.y25 = np.arange(10.).reshape((2, 5)) self.y235 = np.arange(30.).reshape((2, 3, 5)) self.y325 = np.arange(30.).reshape((3, 2, 5)) # Edge updated test matrix 1 # array([[ 30, 1, 2, 3, 4, 5, 6, 7, 8, -30], # [ 30, 11, 12, 13, 14, 15, 16, 17, 18, -30]]) self.y210_edge_updated = np.arange(20.).reshape((2, 10)) self.y210_edge_updated[:, 0] = 30 self.y210_edge_updated[:, -1] = -30 # Edge updated test matrix 2 # array([[ 30, 30], # [ 2, 3], # [ 4, 5], # [ 6, 7], # [ 8, 9], # [ 10, 11], # [ 12, 13], # [ 14, 15], # [ 16, 17], # [-30, -30]]) self.y102_edge_updated = np.arange(20.).reshape((10, 2)) self.y102_edge_updated[0, :] = 30 self.y102_edge_updated[-1, :] = -30 self.fill_value = -100.0 def test_validation(self): # Make sure that appropriate exceptions are raised when invalid values # are given to the constructor. # These should all work. for kind in ('nearest', 'nearest-up', 'zero', 'linear', 'slinear', 'quadratic', 'cubic', 'previous', 'next'): interp1d(self.x10, self.y10, kind=kind) interp1d(self.x10, self.y10, kind=kind, fill_value="extrapolate") interp1d(self.x10, self.y10, kind='linear', fill_value=(-1, 1)) interp1d(self.x10, self.y10, kind='linear', fill_value=np.array([-1])) interp1d(self.x10, self.y10, kind='linear', fill_value=(-1,)) interp1d(self.x10, self.y10, kind='linear', fill_value=-1) interp1d(self.x10, self.y10, kind='linear', fill_value=(-1, -1)) interp1d(self.x10, self.y10, kind=0) interp1d(self.x10, self.y10, kind=1) interp1d(self.x10, self.y10, kind=2) interp1d(self.x10, self.y10, kind=3) interp1d(self.x10, self.y210, kind='linear', axis=-1, fill_value=(-1, -1)) interp1d(self.x2, self.y210, kind='linear', axis=0, fill_value=np.ones(10)) interp1d(self.x2, self.y210, kind='linear', axis=0, fill_value=(np.ones(10), np.ones(10))) interp1d(self.x2, self.y210, kind='linear', axis=0, fill_value=(np.ones(10), -1)) # x array must be 1D. assert_raises(ValueError, interp1d, self.x25, self.y10) # y array cannot be a scalar. assert_raises(ValueError, interp1d, self.x10, np.array(0)) # Check for x and y arrays having the same length. assert_raises(ValueError, interp1d, self.x10, self.y2) assert_raises(ValueError, interp1d, self.x2, self.y10) assert_raises(ValueError, interp1d, self.x10, self.y102) interp1d(self.x10, self.y210) interp1d(self.x10, self.y102, axis=0) # Check for x and y having at least 1 element. assert_raises(ValueError, interp1d, self.x1, self.y10) assert_raises(ValueError, interp1d, self.x10, self.y1) # Bad fill values assert_raises(ValueError, interp1d, self.x10, self.y10, kind='linear', fill_value=(-1, -1, -1)) # doesn't broadcast assert_raises(ValueError, interp1d, self.x10, self.y10, kind='linear', fill_value=[-1, -1, -1]) # doesn't broadcast assert_raises(ValueError, interp1d, self.x10, self.y10, kind='linear', fill_value=np.array((-1, -1, -1))) # doesn't broadcast assert_raises(ValueError, interp1d, self.x10, self.y10, kind='linear', fill_value=[[-1]]) # doesn't broadcast assert_raises(ValueError, interp1d, self.x10, self.y10, kind='linear', fill_value=[-1, -1]) # doesn't broadcast assert_raises(ValueError, interp1d, self.x10, self.y10, kind='linear', fill_value=np.array([])) # doesn't broadcast assert_raises(ValueError, interp1d, self.x10, self.y10, kind='linear', fill_value=()) # doesn't broadcast assert_raises(ValueError, interp1d, self.x2, self.y210, kind='linear', axis=0, fill_value=[-1, -1]) # doesn't broadcast assert_raises(ValueError, interp1d, self.x2, self.y210, kind='linear', axis=0, fill_value=(0., [-1, -1])) # above doesn't bc def test_init(self): # Check that the attributes are initialized appropriately by the # constructor. assert_(interp1d(self.x10, self.y10).copy) assert_(not interp1d(self.x10, self.y10, copy=False).copy) assert_(interp1d(self.x10, self.y10).bounds_error) assert_(not interp1d(self.x10, self.y10, bounds_error=False).bounds_error) assert_(np.isnan(interp1d(self.x10, self.y10).fill_value)) assert_equal(interp1d(self.x10, self.y10, fill_value=3.0).fill_value, 3.0) assert_equal(interp1d(self.x10, self.y10, fill_value=(1.0, 2.0)).fill_value, (1.0, 2.0)) assert_equal(interp1d(self.x10, self.y10).axis, 0) assert_equal(interp1d(self.x10, self.y210).axis, 1) assert_equal(interp1d(self.x10, self.y102, axis=0).axis, 0) assert_array_equal(interp1d(self.x10, self.y10).x, self.x10) assert_array_equal(interp1d(self.x10, self.y10).y, self.y10) assert_array_equal(interp1d(self.x10, self.y210).y, self.y210) def test_assume_sorted(self): # Check for unsorted arrays interp10 = interp1d(self.x10, self.y10) interp10_unsorted = interp1d(self.x10[::-1], self.y10[::-1]) assert_array_almost_equal(interp10_unsorted(self.x10), self.y10) assert_array_almost_equal(interp10_unsorted(1.2), np.array([1.2])) assert_array_almost_equal(interp10_unsorted([2.4, 5.6, 6.0]), interp10([2.4, 5.6, 6.0])) # Check assume_sorted keyword (defaults to False) interp10_assume_kw = interp1d(self.x10[::-1], self.y10[::-1], assume_sorted=False) assert_array_almost_equal(interp10_assume_kw(self.x10), self.y10) interp10_assume_kw2 = interp1d(self.x10[::-1], self.y10[::-1], assume_sorted=True) # Should raise an error for unsorted input if assume_sorted=True assert_raises(ValueError, interp10_assume_kw2, self.x10) # Check that if y is a 2-D array, things are still consistent interp10_y_2d = interp1d(self.x10, self.y210) interp10_y_2d_unsorted = interp1d(self.x10[::-1], self.y210[:, ::-1]) assert_array_almost_equal(interp10_y_2d(self.x10), interp10_y_2d_unsorted(self.x10)) def test_linear(self): for kind in ['linear', 'slinear']: self._check_linear(kind) def _check_linear(self, kind): # Check the actual implementation of linear interpolation. interp10 = interp1d(self.x10, self.y10, kind=kind) assert_array_almost_equal(interp10(self.x10), self.y10) assert_array_almost_equal(interp10(1.2), np.array([1.2])) assert_array_almost_equal(interp10([2.4, 5.6, 6.0]), np.array([2.4, 5.6, 6.0])) # test fill_value="extrapolate" extrapolator = interp1d(self.x10, self.y10, kind=kind, fill_value='extrapolate') assert_allclose(extrapolator([-1., 0, 9, 11]), [-1, 0, 9, 11], rtol=1e-14) opts = dict(kind=kind, fill_value='extrapolate', bounds_error=True) assert_raises(ValueError, interp1d, self.x10, self.y10, **opts) def test_linear_dtypes(self): # regression test for gh-5898, where 1D linear interpolation has been # delegated to numpy.interp for all float dtypes, and the latter was # not handling e.g. np.float128. for dtyp in [np.float16, np.float32, np.float64, np.longdouble]: x = np.arange(8, dtype=dtyp) y = x yp = interp1d(x, y, kind='linear')(x) assert_equal(yp.dtype, dtyp) assert_allclose(yp, y, atol=1e-15) # regression test for gh-14531, where 1D linear interpolation has been # has been extended to delegate to numpy.interp for integer dtypes x = [0, 1, 2] y = [np.nan, 0, 1] yp = interp1d(x, y)(x) assert_allclose(yp, y, atol=1e-15) def test_slinear_dtypes(self): # regression test for gh-7273: 1D slinear interpolation fails with # float32 inputs dt_r = [np.float16, np.float32, np.float64] dt_rc = dt_r + [np.complex64, np.complex128] spline_kinds = ['slinear', 'zero', 'quadratic', 'cubic'] for dtx in dt_r: x = np.arange(0, 10, dtype=dtx) for dty in dt_rc: y = np.exp(-x/3.0).astype(dty) for dtn in dt_r: xnew = x.astype(dtn) for kind in spline_kinds: f = interp1d(x, y, kind=kind, bounds_error=False) assert_allclose(f(xnew), y, atol=1e-7, err_msg=f"{dtx}, {dty} {dtn}") def test_cubic(self): # Check the actual implementation of spline interpolation. interp10 = interp1d(self.x10, self.y10, kind='cubic') assert_array_almost_equal(interp10(self.x10), self.y10) assert_array_almost_equal(interp10(1.2), np.array([1.2])) assert_array_almost_equal(interp10(1.5), np.array([1.5])) assert_array_almost_equal(interp10([2.4, 5.6, 6.0]), np.array([2.4, 5.6, 6.0]),) def test_nearest(self): # Check the actual implementation of nearest-neighbour interpolation. # Nearest asserts that half-integer case (1.5) rounds down to 1 interp10 = interp1d(self.x10, self.y10, kind='nearest') assert_array_almost_equal(interp10(self.x10), self.y10) assert_array_almost_equal(interp10(1.2), np.array(1.)) assert_array_almost_equal(interp10(1.5), np.array(1.)) assert_array_almost_equal(interp10([2.4, 5.6, 6.0]), np.array([2., 6., 6.]),) # test fill_value="extrapolate" extrapolator = interp1d(self.x10, self.y10, kind='nearest', fill_value='extrapolate') assert_allclose(extrapolator([-1., 0, 9, 11]), [0, 0, 9, 9], rtol=1e-14) opts = dict(kind='nearest', fill_value='extrapolate', bounds_error=True) assert_raises(ValueError, interp1d, self.x10, self.y10, **opts) def test_nearest_up(self): # Check the actual implementation of nearest-neighbour interpolation. # Nearest-up asserts that half-integer case (1.5) rounds up to 2 interp10 = interp1d(self.x10, self.y10, kind='nearest-up') assert_array_almost_equal(interp10(self.x10), self.y10) assert_array_almost_equal(interp10(1.2), np.array(1.)) assert_array_almost_equal(interp10(1.5), np.array(2.)) assert_array_almost_equal(interp10([2.4, 5.6, 6.0]), np.array([2., 6., 6.]),) # test fill_value="extrapolate" extrapolator = interp1d(self.x10, self.y10, kind='nearest-up', fill_value='extrapolate') assert_allclose(extrapolator([-1., 0, 9, 11]), [0, 0, 9, 9], rtol=1e-14) opts = dict(kind='nearest-up', fill_value='extrapolate', bounds_error=True) assert_raises(ValueError, interp1d, self.x10, self.y10, **opts) def test_previous(self): # Check the actual implementation of previous interpolation. interp10 = interp1d(self.x10, self.y10, kind='previous') assert_array_almost_equal(interp10(self.x10), self.y10) assert_array_almost_equal(interp10(1.2), np.array(1.)) assert_array_almost_equal(interp10(1.5), np.array(1.)) assert_array_almost_equal(interp10([2.4, 5.6, 6.0]), np.array([2., 5., 6.]),) # test fill_value="extrapolate" extrapolator = interp1d(self.x10, self.y10, kind='previous', fill_value='extrapolate') assert_allclose(extrapolator([-1., 0, 9, 11]), [np.nan, 0, 9, 9], rtol=1e-14) # Tests for gh-9591 interpolator1D = interp1d(self.x10, self.y10, kind="previous", fill_value='extrapolate') assert_allclose(interpolator1D([-1, -2, 5, 8, 12, 25]), [np.nan, np.nan, 5, 8, 9, 9]) interpolator2D = interp1d(self.x10, self.y210, kind="previous", fill_value='extrapolate') assert_allclose(interpolator2D([-1, -2, 5, 8, 12, 25]), [[np.nan, np.nan, 5, 8, 9, 9], [np.nan, np.nan, 15, 18, 19, 19]]) interpolator2DAxis0 = interp1d(self.x10, self.y102, kind="previous", axis=0, fill_value='extrapolate') assert_allclose(interpolator2DAxis0([-2, 5, 12]), [[np.nan, np.nan], [10, 11], [18, 19]]) opts = dict(kind='previous', fill_value='extrapolate', bounds_error=True) assert_raises(ValueError, interp1d, self.x10, self.y10, **opts) # Tests for gh-16813 interpolator1D = interp1d([0, 1, 2], [0, 1, -1], kind="previous", fill_value='extrapolate', assume_sorted=True) assert_allclose(interpolator1D([-2, -1, 0, 1, 2, 3, 5]), [np.nan, np.nan, 0, 1, -1, -1, -1]) interpolator1D = interp1d([2, 0, 1], # x is not ascending [-1, 0, 1], kind="previous", fill_value='extrapolate', assume_sorted=False) assert_allclose(interpolator1D([-2, -1, 0, 1, 2, 3, 5]), [np.nan, np.nan, 0, 1, -1, -1, -1]) interpolator2D = interp1d(self.x10, self.y210_edge_updated, kind="previous", fill_value='extrapolate') assert_allclose(interpolator2D([-1, -2, 5, 8, 12, 25]), [[np.nan, np.nan, 5, 8, -30, -30], [np.nan, np.nan, 15, 18, -30, -30]]) interpolator2DAxis0 = interp1d(self.x10, self.y102_edge_updated, kind="previous", axis=0, fill_value='extrapolate') assert_allclose(interpolator2DAxis0([-2, 5, 12]), [[np.nan, np.nan], [10, 11], [-30, -30]]) def test_next(self): # Check the actual implementation of next interpolation. interp10 = interp1d(self.x10, self.y10, kind='next') assert_array_almost_equal(interp10(self.x10), self.y10) assert_array_almost_equal(interp10(1.2), np.array(2.)) assert_array_almost_equal(interp10(1.5), np.array(2.)) assert_array_almost_equal(interp10([2.4, 5.6, 6.0]), np.array([3., 6., 6.]),) # test fill_value="extrapolate" extrapolator = interp1d(self.x10, self.y10, kind='next', fill_value='extrapolate') assert_allclose(extrapolator([-1., 0, 9, 11]), [0, 0, 9, np.nan], rtol=1e-14) # Tests for gh-9591 interpolator1D = interp1d(self.x10, self.y10, kind="next", fill_value='extrapolate') assert_allclose(interpolator1D([-1, -2, 5, 8, 12, 25]), [0, 0, 5, 8, np.nan, np.nan]) interpolator2D = interp1d(self.x10, self.y210, kind="next", fill_value='extrapolate') assert_allclose(interpolator2D([-1, -2, 5, 8, 12, 25]), [[0, 0, 5, 8, np.nan, np.nan], [10, 10, 15, 18, np.nan, np.nan]]) interpolator2DAxis0 = interp1d(self.x10, self.y102, kind="next", axis=0, fill_value='extrapolate') assert_allclose(interpolator2DAxis0([-2, 5, 12]), [[0, 1], [10, 11], [np.nan, np.nan]]) opts = dict(kind='next', fill_value='extrapolate', bounds_error=True) assert_raises(ValueError, interp1d, self.x10, self.y10, **opts) # Tests for gh-16813 interpolator1D = interp1d([0, 1, 2], [0, 1, -1], kind="next", fill_value='extrapolate', assume_sorted=True) assert_allclose(interpolator1D([-2, -1, 0, 1, 2, 3, 5]), [0, 0, 0, 1, -1, np.nan, np.nan]) interpolator1D = interp1d([2, 0, 1], # x is not ascending [-1, 0, 1], kind="next", fill_value='extrapolate', assume_sorted=False) assert_allclose(interpolator1D([-2, -1, 0, 1, 2, 3, 5]), [0, 0, 0, 1, -1, np.nan, np.nan]) interpolator2D = interp1d(self.x10, self.y210_edge_updated, kind="next", fill_value='extrapolate') assert_allclose(interpolator2D([-1, -2, 5, 8, 12, 25]), [[30, 30, 5, 8, np.nan, np.nan], [30, 30, 15, 18, np.nan, np.nan]]) interpolator2DAxis0 = interp1d(self.x10, self.y102_edge_updated, kind="next", axis=0, fill_value='extrapolate') assert_allclose(interpolator2DAxis0([-2, 5, 12]), [[30, 30], [10, 11], [np.nan, np.nan]]) def test_zero(self): # Check the actual implementation of zero-order spline interpolation. interp10 = interp1d(self.x10, self.y10, kind='zero') assert_array_almost_equal(interp10(self.x10), self.y10) assert_array_almost_equal(interp10(1.2), np.array(1.)) assert_array_almost_equal(interp10(1.5), np.array(1.)) assert_array_almost_equal(interp10([2.4, 5.6, 6.0]), np.array([2., 5., 6.])) def bounds_check_helper(self, interpolant, test_array, fail_value): # Asserts that a ValueError is raised and that the error message # contains the value causing this exception. assert_raises(ValueError, interpolant, test_array) try: interpolant(test_array) except ValueError as err: assert (f"{fail_value}" in str(err)) def _bounds_check(self, kind='linear'): # Test that our handling of out-of-bounds input is correct. extrap10 = interp1d(self.x10, self.y10, fill_value=self.fill_value, bounds_error=False, kind=kind) assert_array_equal(extrap10(11.2), np.array(self.fill_value)) assert_array_equal(extrap10(-3.4), np.array(self.fill_value)) assert_array_equal(extrap10([[[11.2], [-3.4], [12.6], [19.3]]]), np.array(self.fill_value),) assert_array_equal(extrap10._check_bounds( np.array([-1.0, 0.0, 5.0, 9.0, 11.0])), np.array([[True, False, False, False, False], [False, False, False, False, True]])) raises_bounds_error = interp1d(self.x10, self.y10, bounds_error=True, kind=kind) self.bounds_check_helper(raises_bounds_error, -1.0, -1.0) self.bounds_check_helper(raises_bounds_error, 11.0, 11.0) self.bounds_check_helper(raises_bounds_error, [0.0, -1.0, 0.0], -1.0) self.bounds_check_helper(raises_bounds_error, [0.0, 1.0, 21.0], 21.0) raises_bounds_error([0.0, 5.0, 9.0]) def _bounds_check_int_nan_fill(self, kind='linear'): x = np.arange(10).astype(int) y = np.arange(10).astype(int) c = interp1d(x, y, kind=kind, fill_value=np.nan, bounds_error=False) yi = c(x - 1) assert_(np.isnan(yi[0])) assert_array_almost_equal(yi, np.r_[np.nan, y[:-1]]) def test_bounds(self): for kind in ('linear', 'cubic', 'nearest', 'previous', 'next', 'slinear', 'zero', 'quadratic'): self._bounds_check(kind) self._bounds_check_int_nan_fill(kind) def _check_fill_value(self, kind): interp = interp1d(self.x10, self.y10, kind=kind, fill_value=(-100, 100), bounds_error=False) assert_array_almost_equal(interp(10), 100) assert_array_almost_equal(interp(-10), -100) assert_array_almost_equal(interp([-10, 10]), [-100, 100]) # Proper broadcasting: # interp along axis of length 5 # other dim=(2, 3), (3, 2), (2, 2), or (2,) # one singleton fill_value (works for all) for y in (self.y235, self.y325, self.y225, self.y25): interp = interp1d(self.x5, y, kind=kind, axis=-1, fill_value=100, bounds_error=False) assert_array_almost_equal(interp(10), 100) assert_array_almost_equal(interp(-10), 100) assert_array_almost_equal(interp([-10, 10]), 100) # singleton lower, singleton upper interp = interp1d(self.x5, y, kind=kind, axis=-1, fill_value=(-100, 100), bounds_error=False) assert_array_almost_equal(interp(10), 100) assert_array_almost_equal(interp(-10), -100) if y.ndim == 3: result = [[[-100, 100]] * y.shape[1]] * y.shape[0] else: result = [[-100, 100]] * y.shape[0] assert_array_almost_equal(interp([-10, 10]), result) # one broadcastable (3,) fill_value fill_value = [100, 200, 300] for y in (self.y325, self.y225): assert_raises(ValueError, interp1d, self.x5, y, kind=kind, axis=-1, fill_value=fill_value, bounds_error=False) interp = interp1d(self.x5, self.y235, kind=kind, axis=-1, fill_value=fill_value, bounds_error=False) assert_array_almost_equal(interp(10), [[100, 200, 300]] * 2) assert_array_almost_equal(interp(-10), [[100, 200, 300]] * 2) assert_array_almost_equal(interp([-10, 10]), [[[100, 100], [200, 200], [300, 300]]] * 2) # one broadcastable (2,) fill_value fill_value = [100, 200] assert_raises(ValueError, interp1d, self.x5, self.y235, kind=kind, axis=-1, fill_value=fill_value, bounds_error=False) for y in (self.y225, self.y325, self.y25): interp = interp1d(self.x5, y, kind=kind, axis=-1, fill_value=fill_value, bounds_error=False) result = [100, 200] if y.ndim == 3: result = [result] * y.shape[0] assert_array_almost_equal(interp(10), result) assert_array_almost_equal(interp(-10), result) result = [[100, 100], [200, 200]] if y.ndim == 3: result = [result] * y.shape[0] assert_array_almost_equal(interp([-10, 10]), result) # broadcastable (3,) lower, singleton upper fill_value = (np.array([-100, -200, -300]), 100) for y in (self.y325, self.y225): assert_raises(ValueError, interp1d, self.x5, y, kind=kind, axis=-1, fill_value=fill_value, bounds_error=False) interp = interp1d(self.x5, self.y235, kind=kind, axis=-1, fill_value=fill_value, bounds_error=False) assert_array_almost_equal(interp(10), 100) assert_array_almost_equal(interp(-10), [[-100, -200, -300]] * 2) assert_array_almost_equal(interp([-10, 10]), [[[-100, 100], [-200, 100], [-300, 100]]] * 2) # broadcastable (2,) lower, singleton upper fill_value = (np.array([-100, -200]), 100) assert_raises(ValueError, interp1d, self.x5, self.y235, kind=kind, axis=-1, fill_value=fill_value, bounds_error=False) for y in (self.y225, self.y325, self.y25): interp = interp1d(self.x5, y, kind=kind, axis=-1, fill_value=fill_value, bounds_error=False) assert_array_almost_equal(interp(10), 100) result = [-100, -200] if y.ndim == 3: result = [result] * y.shape[0] assert_array_almost_equal(interp(-10), result) result = [[-100, 100], [-200, 100]] if y.ndim == 3: result = [result] * y.shape[0] assert_array_almost_equal(interp([-10, 10]), result) # broadcastable (3,) lower, broadcastable (3,) upper fill_value = ([-100, -200, -300], [100, 200, 300]) for y in (self.y325, self.y225): assert_raises(ValueError, interp1d, self.x5, y, kind=kind, axis=-1, fill_value=fill_value, bounds_error=False) for ii in range(2): # check ndarray as well as list here if ii == 1: fill_value = tuple(np.array(f) for f in fill_value) interp = interp1d(self.x5, self.y235, kind=kind, axis=-1, fill_value=fill_value, bounds_error=False) assert_array_almost_equal(interp(10), [[100, 200, 300]] * 2) assert_array_almost_equal(interp(-10), [[-100, -200, -300]] * 2) assert_array_almost_equal(interp([-10, 10]), [[[-100, 100], [-200, 200], [-300, 300]]] * 2) # broadcastable (2,) lower, broadcastable (2,) upper fill_value = ([-100, -200], [100, 200]) assert_raises(ValueError, interp1d, self.x5, self.y235, kind=kind, axis=-1, fill_value=fill_value, bounds_error=False) for y in (self.y325, self.y225, self.y25): interp = interp1d(self.x5, y, kind=kind, axis=-1, fill_value=fill_value, bounds_error=False) result = [100, 200] if y.ndim == 3: result = [result] * y.shape[0] assert_array_almost_equal(interp(10), result) result = [-100, -200] if y.ndim == 3: result = [result] * y.shape[0] assert_array_almost_equal(interp(-10), result) result = [[-100, 100], [-200, 200]] if y.ndim == 3: result = [result] * y.shape[0] assert_array_almost_equal(interp([-10, 10]), result) # one broadcastable (2, 2) array-like fill_value = [[100, 200], [1000, 2000]] for y in (self.y235, self.y325, self.y25): assert_raises(ValueError, interp1d, self.x5, y, kind=kind, axis=-1, fill_value=fill_value, bounds_error=False) for ii in range(2): if ii == 1: fill_value = np.array(fill_value) interp = interp1d(self.x5, self.y225, kind=kind, axis=-1, fill_value=fill_value, bounds_error=False) assert_array_almost_equal(interp(10), [[100, 200], [1000, 2000]]) assert_array_almost_equal(interp(-10), [[100, 200], [1000, 2000]]) assert_array_almost_equal(interp([-10, 10]), [[[100, 100], [200, 200]], [[1000, 1000], [2000, 2000]]]) # broadcastable (2, 2) lower, broadcastable (2, 2) upper fill_value = ([[-100, -200], [-1000, -2000]], [[100, 200], [1000, 2000]]) for y in (self.y235, self.y325, self.y25): assert_raises(ValueError, interp1d, self.x5, y, kind=kind, axis=-1, fill_value=fill_value, bounds_error=False) for ii in range(2): if ii == 1: fill_value = (np.array(fill_value[0]), np.array(fill_value[1])) interp = interp1d(self.x5, self.y225, kind=kind, axis=-1, fill_value=fill_value, bounds_error=False) assert_array_almost_equal(interp(10), [[100, 200], [1000, 2000]]) assert_array_almost_equal(interp(-10), [[-100, -200], [-1000, -2000]]) assert_array_almost_equal(interp([-10, 10]), [[[-100, 100], [-200, 200]], [[-1000, 1000], [-2000, 2000]]]) def test_fill_value(self): # test that two-element fill value works for kind in ('linear', 'nearest', 'cubic', 'slinear', 'quadratic', 'zero', 'previous', 'next'): self._check_fill_value(kind) def test_fill_value_writeable(self): # backwards compat: fill_value is a public writeable attribute interp = interp1d(self.x10, self.y10, fill_value=123.0) assert_equal(interp.fill_value, 123.0) interp.fill_value = 321.0 assert_equal(interp.fill_value, 321.0) def _nd_check_interp(self, kind='linear'): # Check the behavior when the inputs and outputs are multidimensional. # Multidimensional input. interp10 = interp1d(self.x10, self.y10, kind=kind) assert_array_almost_equal(interp10(np.array([[3., 5.], [2., 7.]])), np.array([[3., 5.], [2., 7.]])) # Scalar input -> 0-dim scalar array output assert_(isinstance(interp10(1.2), np.ndarray)) assert_equal(interp10(1.2).shape, ()) # Multidimensional outputs. interp210 = interp1d(self.x10, self.y210, kind=kind) assert_array_almost_equal(interp210(1.), np.array([1., 11.])) assert_array_almost_equal(interp210(np.array([1., 2.])), np.array([[1., 2.], [11., 12.]])) interp102 = interp1d(self.x10, self.y102, axis=0, kind=kind) assert_array_almost_equal(interp102(1.), np.array([2.0, 3.0])) assert_array_almost_equal(interp102(np.array([1., 3.])), np.array([[2., 3.], [6., 7.]])) # Both at the same time! x_new = np.array([[3., 5.], [2., 7.]]) assert_array_almost_equal(interp210(x_new), np.array([[[3., 5.], [2., 7.]], [[13., 15.], [12., 17.]]])) assert_array_almost_equal(interp102(x_new), np.array([[[6., 7.], [10., 11.]], [[4., 5.], [14., 15.]]])) def _nd_check_shape(self, kind='linear'): # Check large N-D output shape a = [4, 5, 6, 7] y = np.arange(np.prod(a)).reshape(*a) for n, s in enumerate(a): x = np.arange(s) z = interp1d(x, y, axis=n, kind=kind) assert_array_almost_equal(z(x), y, err_msg=kind) x2 = np.arange(2*3*1).reshape((2,3,1)) / 12. b = list(a) b[n:n+1] = [2,3,1] assert_array_almost_equal(z(x2).shape, b, err_msg=kind) def test_nd(self): for kind in ('linear', 'cubic', 'slinear', 'quadratic', 'nearest', 'zero', 'previous', 'next'): self._nd_check_interp(kind) self._nd_check_shape(kind) def _check_complex(self, dtype=np.complex128, kind='linear'): x = np.array([1, 2.5, 3, 3.1, 4, 6.4, 7.9, 8.0, 9.5, 10]) y = x * x ** (1 + 2j) y = y.astype(dtype) # simple test c = interp1d(x, y, kind=kind) assert_array_almost_equal(y[:-1], c(x)[:-1]) # check against interpolating real+imag separately xi = np.linspace(1, 10, 31) cr = interp1d(x, y.real, kind=kind) ci = interp1d(x, y.imag, kind=kind) assert_array_almost_equal(c(xi).real, cr(xi)) assert_array_almost_equal(c(xi).imag, ci(xi)) def test_complex(self): for kind in ('linear', 'nearest', 'cubic', 'slinear', 'quadratic', 'zero', 'previous', 'next'): self._check_complex(np.complex64, kind) self._check_complex(np.complex128, kind) @pytest.mark.skipif(IS_PYPY, reason="Test not meaningful on PyPy") def test_circular_refs(self): # Test interp1d can be automatically garbage collected x = np.linspace(0, 1) y = np.linspace(0, 1) # Confirm interp can be released from memory after use with assert_deallocated(interp1d, x, y) as interp: interp([0.1, 0.2]) del interp def test_overflow_nearest(self): # Test that the x range doesn't overflow when given integers as input for kind in ('nearest', 'previous', 'next'): x = np.array([0, 50, 127], dtype=np.int8) ii = interp1d(x, x, kind=kind) assert_array_almost_equal(ii(x), x) def test_local_nans(self): # check that for local interpolation kinds (slinear, zero) a single nan # only affects its local neighborhood x = np.arange(10).astype(float) y = x.copy() y[6] = np.nan for kind in ('zero', 'slinear'): ir = interp1d(x, y, kind=kind) vals = ir([4.9, 7.0]) assert_(np.isfinite(vals).all()) def test_spline_nans(self): # Backwards compat: a single nan makes the whole spline interpolation # return nans in an array of the correct shape. And it doesn't raise, # just quiet nans because of backcompat. x = np.arange(8).astype(float) y = x.copy() yn = y.copy() yn[3] = np.nan for kind in ['quadratic', 'cubic']: ir = interp1d(x, y, kind=kind) irn = interp1d(x, yn, kind=kind) for xnew in (6, [1, 6], [[1, 6], [3, 5]]): xnew = np.asarray(xnew) out, outn = ir(x), irn(x) assert_(np.isnan(outn).all()) assert_equal(out.shape, outn.shape) def test_all_nans(self): # regression test for gh-11637: interp1d core dumps with all-nan `x` x = np.ones(10) * np.nan y = np.arange(10) with assert_raises(ValueError): interp1d(x, y, kind='cubic') def test_read_only(self): x = np.arange(0, 10) y = np.exp(-x / 3.0) xnew = np.arange(0, 9, 0.1) # Check both read-only and not read-only: for xnew_writeable in (True, False): xnew.flags.writeable = xnew_writeable x.flags.writeable = False for kind in ('linear', 'nearest', 'zero', 'slinear', 'quadratic', 'cubic'): f = interp1d(x, y, kind=kind) vals = f(xnew) assert_(np.isfinite(vals).all()) @pytest.mark.parametrize( "kind", ("linear", "nearest", "nearest-up", "previous", "next") ) def test_single_value(self, kind): # https://github.com/scipy/scipy/issues/4043 f = interp1d([1.5], [6], kind=kind, bounds_error=False, fill_value=(2, 10)) assert_array_equal(f([1, 1.5, 2]), [2, 6, 10]) # check still error if bounds_error=True f = interp1d([1.5], [6], kind=kind, bounds_error=True) with assert_raises(ValueError, match="x_new is above"): f(2.0) class TestLagrange: def test_lagrange(self): p = poly1d([5,2,1,4,3]) xs = np.arange(len(p.coeffs)) ys = p(xs) pl = lagrange(xs,ys) assert_array_almost_equal(p.coeffs,pl.coeffs) class TestAkima1DInterpolator: def test_eval(self): x = np.arange(0., 11.) y = np.array([0., 2., 1., 3., 2., 6., 5.5, 5.5, 2.7, 5.1, 3.]) ak = Akima1DInterpolator(x, y) xi = np.array([0., 0.5, 1., 1.5, 2.5, 3.5, 4.5, 5.1, 6.5, 7.2, 8.6, 9.9, 10.]) yi = np.array([0., 1.375, 2., 1.5, 1.953125, 2.484375, 4.1363636363636366866103344, 5.9803623910336236590978842, 5.5067291516462386624652936, 5.2031367459745245795943447, 4.1796554159017080820603951, 3.4110386597938129327189927, 3.]) assert_allclose(ak(xi), yi) def test_eval_mod(self): # Reference values generated with the following MATLAB code: # format longG # x = 0:10; y = [0. 2. 1. 3. 2. 6. 5.5 5.5 2.7 5.1 3.]; # xi = [0. 0.5 1. 1.5 2.5 3.5 4.5 5.1 6.5 7.2 8.6 9.9 10.]; # makima(x, y, xi) x = np.arange(0., 11.) y = np.array([0., 2., 1., 3., 2., 6., 5.5, 5.5, 2.7, 5.1, 3.]) ak = Akima1DInterpolator(x, y, method="makima") xi = np.array([0., 0.5, 1., 1.5, 2.5, 3.5, 4.5, 5.1, 6.5, 7.2, 8.6, 9.9, 10.]) yi = np.array([ 0.0, 1.34471153846154, 2.0, 1.44375, 1.94375, 2.51939102564103, 4.10366931918656, 5.98501550899192, 5.51756330960439, 5.1757231914014, 4.12326636931311, 3.32931513157895, 3.0]) assert_allclose(ak(xi), yi) def test_eval_2d(self): x = np.arange(0., 11.) y = np.array([0., 2., 1., 3., 2., 6., 5.5, 5.5, 2.7, 5.1, 3.]) y = np.column_stack((y, 2. * y)) ak = Akima1DInterpolator(x, y) xi = np.array([0., 0.5, 1., 1.5, 2.5, 3.5, 4.5, 5.1, 6.5, 7.2, 8.6, 9.9, 10.]) yi = np.array([0., 1.375, 2., 1.5, 1.953125, 2.484375, 4.1363636363636366866103344, 5.9803623910336236590978842, 5.5067291516462386624652936, 5.2031367459745245795943447, 4.1796554159017080820603951, 3.4110386597938129327189927, 3.]) yi = np.column_stack((yi, 2. * yi)) assert_allclose(ak(xi), yi) def test_eval_3d(self): x = np.arange(0., 11.) y_ = np.array([0., 2., 1., 3., 2., 6., 5.5, 5.5, 2.7, 5.1, 3.]) y = np.empty((11, 2, 2)) y[:, 0, 0] = y_ y[:, 1, 0] = 2. * y_ y[:, 0, 1] = 3. * y_ y[:, 1, 1] = 4. * y_ ak = Akima1DInterpolator(x, y) xi = np.array([0., 0.5, 1., 1.5, 2.5, 3.5, 4.5, 5.1, 6.5, 7.2, 8.6, 9.9, 10.]) yi = np.empty((13, 2, 2)) yi_ = np.array([0., 1.375, 2., 1.5, 1.953125, 2.484375, 4.1363636363636366866103344, 5.9803623910336236590978842, 5.5067291516462386624652936, 5.2031367459745245795943447, 4.1796554159017080820603951, 3.4110386597938129327189927, 3.]) yi[:, 0, 0] = yi_ yi[:, 1, 0] = 2. * yi_ yi[:, 0, 1] = 3. * yi_ yi[:, 1, 1] = 4. * yi_ assert_allclose(ak(xi), yi) def test_degenerate_case_multidimensional(self): # This test is for issue #5683. x = np.array([0, 1, 2]) y = np.vstack((x, x**2)).T ak = Akima1DInterpolator(x, y) x_eval = np.array([0.5, 1.5]) y_eval = ak(x_eval) assert_allclose(y_eval, np.vstack((x_eval, x_eval**2)).T) def test_extend(self): x = np.arange(0., 11.) y = np.array([0., 2., 1., 3., 2., 6., 5.5, 5.5, 2.7, 5.1, 3.]) ak = Akima1DInterpolator(x, y) match = "Extending a 1-D Akima interpolator is not yet implemented" with pytest.raises(NotImplementedError, match=match): ak.extend(None, None) def test_mod_invalid_method(self): x = np.arange(0., 11.) y = np.array([0., 2., 1., 3., 2., 6., 5.5, 5.5, 2.7, 5.1, 3.]) match = "`method`=invalid is unsupported." with pytest.raises(NotImplementedError, match=match): Akima1DInterpolator(x, y, method="invalid") # type: ignore def test_complex(self): # Complex-valued data deprecated x = np.arange(0., 11.) y = np.array([0., 2., 1., 3., 2., 6., 5.5, 5.5, 2.7, 5.1, 3.]) y = y - 2j*y # actually raises ComplexWarning, which subclasses RuntimeWarning, see # https://github.com/numpy/numpy/blob/main/numpy/exceptions.py msg = "Passing an array with a complex.*|Casting complex values to real.*" with pytest.warns((RuntimeWarning, DeprecationWarning), match=msg): Akima1DInterpolator(x, y) class TestPPolyCommon: # test basic functionality for PPoly and BPoly def test_sort_check(self): c = np.array([[1, 4], [2, 5], [3, 6]]) x = np.array([0, 1, 0.5]) assert_raises(ValueError, PPoly, c, x) assert_raises(ValueError, BPoly, c, x) def test_ctor_c(self): # wrong shape: `c` must be at least 2D with assert_raises(ValueError): PPoly([1, 2], [0, 1]) def test_extend(self): # Test adding new points to the piecewise polynomial np.random.seed(1234) order = 3 x = np.unique(np.r_[0, 10 * np.random.rand(30), 10]) c = 2*np.random.rand(order+1, len(x)-1, 2, 3) - 1 for cls in (PPoly, BPoly): pp = cls(c[:,:9], x[:10]) pp.extend(c[:,9:], x[10:]) pp2 = cls(c[:, 10:], x[10:]) pp2.extend(c[:, :10], x[:10]) pp3 = cls(c, x) assert_array_equal(pp.c, pp3.c) assert_array_equal(pp.x, pp3.x) assert_array_equal(pp2.c, pp3.c) assert_array_equal(pp2.x, pp3.x) def test_extend_diff_orders(self): # Test extending polynomial with different order one np.random.seed(1234) x = np.linspace(0, 1, 6) c = np.random.rand(2, 5) x2 = np.linspace(1, 2, 6) c2 = np.random.rand(4, 5) for cls in (PPoly, BPoly): pp1 = cls(c, x) pp2 = cls(c2, x2) pp_comb = cls(c, x) pp_comb.extend(c2, x2[1:]) # NB. doesn't match to pp1 at the endpoint, because pp1 is not # continuous with pp2 as we took random coefs. xi1 = np.linspace(0, 1, 300, endpoint=False) xi2 = np.linspace(1, 2, 300) assert_allclose(pp1(xi1), pp_comb(xi1)) assert_allclose(pp2(xi2), pp_comb(xi2)) def test_extend_descending(self): np.random.seed(0) order = 3 x = np.sort(np.random.uniform(0, 10, 20)) c = np.random.rand(order + 1, x.shape[0] - 1, 2, 3) for cls in (PPoly, BPoly): p = cls(c, x) p1 = cls(c[:, :9], x[:10]) p1.extend(c[:, 9:], x[10:]) p2 = cls(c[:, 10:], x[10:]) p2.extend(c[:, :10], x[:10]) assert_array_equal(p1.c, p.c) assert_array_equal(p1.x, p.x) assert_array_equal(p2.c, p.c) assert_array_equal(p2.x, p.x) def test_shape(self): np.random.seed(1234) c = np.random.rand(8, 12, 5, 6, 7) x = np.sort(np.random.rand(13)) xp = np.random.rand(3, 4) for cls in (PPoly, BPoly): p = cls(c, x) assert_equal(p(xp).shape, (3, 4, 5, 6, 7)) # 'scalars' for cls in (PPoly, BPoly): p = cls(c[..., 0, 0, 0], x) assert_equal(np.shape(p(0.5)), ()) assert_equal(np.shape(p(np.array(0.5))), ()) assert_raises(ValueError, p, np.array([[0.1, 0.2], [0.4]], dtype=object)) 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)) def test_axis(self): np.random.seed(12345) c = np.random.rand(3, 4, 5, 6, 7, 8) c_s = c.shape xp = np.random.random((1, 2)) for axis in (0, 1, 2, 3): m = c.shape[axis+1] x = np.sort(np.random.rand(m+1)) for cls in (PPoly, BPoly): p = cls(c, x, axis=axis) assert_equal(p.c.shape, c_s[axis:axis+2] + c_s[:axis] + c_s[axis+2:]) res = p(xp) targ_shape = c_s[:axis] + xp.shape + c_s[2+axis:] assert_equal(res.shape, targ_shape) # deriv/antideriv does not drop the axis for p1 in [cls(c, x, axis=axis).derivative(), cls(c, x, axis=axis).derivative(2), cls(c, x, axis=axis).antiderivative(), cls(c, x, axis=axis).antiderivative(2)]: assert_equal(p1.axis, p.axis) # c array needs two axes for the coefficients and intervals, so # 0 <= axis < c.ndim-1; raise otherwise for axis in (-1, 4, 5, 6): for cls in (BPoly, PPoly): assert_raises(ValueError, cls, **dict(c=c, x=x, axis=axis)) class TestPolySubclassing: class P(PPoly): pass class B(BPoly): pass def _make_polynomials(self): np.random.seed(1234) x = np.sort(np.random.random(3)) c = np.random.random((4, 2)) return self.P(c, x), self.B(c, x) def test_derivative(self): pp, bp = self._make_polynomials() for p in (pp, bp): pd = p.derivative() assert_equal(p.__class__, pd.__class__) ppa = pp.antiderivative() assert_equal(pp.__class__, ppa.__class__) def test_from_spline(self): np.random.seed(1234) x = np.sort(np.r_[0, np.random.rand(11), 1]) y = np.random.rand(len(x)) spl = splrep(x, y, s=0) pp = self.P.from_spline(spl) assert_equal(pp.__class__, self.P) def test_conversions(self): pp, bp = self._make_polynomials() pp1 = self.P.from_bernstein_basis(bp) assert_equal(pp1.__class__, self.P) bp1 = self.B.from_power_basis(pp) assert_equal(bp1.__class__, self.B) def test_from_derivatives(self): x = [0, 1, 2] y = [[1], [2], [3]] bp = self.B.from_derivatives(x, y) assert_equal(bp.__class__, self.B) class TestPPoly: def test_simple(self): c = np.array([[1, 4], [2, 5], [3, 6]]) x = np.array([0, 0.5, 1]) p = PPoly(c, x) assert_allclose(p(0.3), 1*0.3**2 + 2*0.3 + 3) assert_allclose(p(0.7), 4*(0.7-0.5)**2 + 5*(0.7-0.5) + 6) def test_periodic(self): c = np.array([[1, 4], [2, 5], [3, 6]]) x = np.array([0, 0.5, 1]) p = PPoly(c, x, extrapolate='periodic') assert_allclose(p(1.3), 1 * 0.3 ** 2 + 2 * 0.3 + 3) assert_allclose(p(-0.3), 4 * (0.7 - 0.5) ** 2 + 5 * (0.7 - 0.5) + 6) assert_allclose(p(1.3, 1), 2 * 0.3 + 2) assert_allclose(p(-0.3, 1), 8 * (0.7 - 0.5) + 5) def test_read_only(self): c = np.array([[1, 4], [2, 5], [3, 6]]) x = np.array([0, 0.5, 1]) xnew = np.array([0, 0.1, 0.2]) PPoly(c, x, extrapolate='periodic') for writeable in (True, False): x.flags.writeable = writeable c.flags.writeable = writeable f = PPoly(c, x) vals = f(xnew) assert_(np.isfinite(vals).all()) def test_descending(self): def binom_matrix(power): n = np.arange(power + 1).reshape(-1, 1) k = np.arange(power + 1) B = binom(n, k) return B[::-1, ::-1] np.random.seed(0) power = 3 for m in [10, 20, 30]: x = np.sort(np.random.uniform(0, 10, m + 1)) ca = np.random.uniform(-2, 2, size=(power + 1, m)) h = np.diff(x) h_powers = h[None, :] ** np.arange(power + 1)[::-1, None] B = binom_matrix(power) cap = ca * h_powers cdp = np.dot(B.T, cap) cd = cdp / h_powers pa = PPoly(ca, x, extrapolate=True) pd = PPoly(cd[:, ::-1], x[::-1], extrapolate=True) x_test = np.random.uniform(-10, 20, 100) assert_allclose(pa(x_test), pd(x_test), rtol=1e-13) assert_allclose(pa(x_test, 1), pd(x_test, 1), rtol=1e-13) pa_d = pa.derivative() pd_d = pd.derivative() assert_allclose(pa_d(x_test), pd_d(x_test), rtol=1e-13) # Antiderivatives won't be equal because fixing continuity is # done in the reverse order, but surely the differences should be # equal. pa_i = pa.antiderivative() pd_i = pd.antiderivative() for a, b in np.random.uniform(-10, 20, (5, 2)): int_a = pa.integrate(a, b) int_d = pd.integrate(a, b) assert_allclose(int_a, int_d, rtol=1e-13) assert_allclose(pa_i(b) - pa_i(a), pd_i(b) - pd_i(a), rtol=1e-13) roots_d = pd.roots() roots_a = pa.roots() assert_allclose(roots_a, np.sort(roots_d), rtol=1e-12) def test_multi_shape(self): c = np.random.rand(6, 2, 1, 2, 3) x = np.array([0, 0.5, 1]) p = PPoly(c, x) assert_equal(p.x.shape, x.shape) assert_equal(p.c.shape, c.shape) assert_equal(p(0.3).shape, c.shape[2:]) assert_equal(p(np.random.rand(5, 6)).shape, (5, 6) + c.shape[2:]) dp = p.derivative() assert_equal(dp.c.shape, (5, 2, 1, 2, 3)) ip = p.antiderivative() assert_equal(ip.c.shape, (7, 2, 1, 2, 3)) def test_construct_fast(self): np.random.seed(1234) c = np.array([[1, 4], [2, 5], [3, 6]], dtype=float) x = np.array([0, 0.5, 1]) p = PPoly.construct_fast(c, x) assert_allclose(p(0.3), 1*0.3**2 + 2*0.3 + 3) assert_allclose(p(0.7), 4*(0.7-0.5)**2 + 5*(0.7-0.5) + 6) def test_vs_alternative_implementations(self): np.random.seed(1234) c = np.random.rand(3, 12, 22) x = np.sort(np.r_[0, np.random.rand(11), 1]) p = PPoly(c, x) xp = np.r_[0.3, 0.5, 0.33, 0.6] expected = _ppoly_eval_1(c, x, xp) assert_allclose(p(xp), expected) expected = _ppoly_eval_2(c[:,:,0], x, xp) assert_allclose(p(xp)[:,0], expected) def test_from_spline(self): np.random.seed(1234) x = np.sort(np.r_[0, np.random.rand(11), 1]) y = np.random.rand(len(x)) spl = splrep(x, y, s=0) pp = PPoly.from_spline(spl) xi = np.linspace(0, 1, 200) assert_allclose(pp(xi), splev(xi, spl)) # make sure .from_spline accepts BSpline objects b = BSpline(*spl) ppp = PPoly.from_spline(b) assert_allclose(ppp(xi), b(xi)) # BSpline's extrapolate attribute propagates unless overridden t, c, k = spl for extrap in (None, True, False): b = BSpline(t, c, k, extrapolate=extrap) p = PPoly.from_spline(b) assert_equal(p.extrapolate, b.extrapolate) def test_derivative_simple(self): np.random.seed(1234) c = np.array([[4, 3, 2, 1]]).T dc = np.array([[3*4, 2*3, 2]]).T ddc = np.array([[2*3*4, 1*2*3]]).T x = np.array([0, 1]) pp = PPoly(c, x) dpp = PPoly(dc, x) ddpp = PPoly(ddc, x) assert_allclose(pp.derivative().c, dpp.c) assert_allclose(pp.derivative(2).c, ddpp.c) def test_derivative_eval(self): np.random.seed(1234) x = np.sort(np.r_[0, np.random.rand(11), 1]) y = np.random.rand(len(x)) spl = splrep(x, y, s=0) pp = PPoly.from_spline(spl) xi = np.linspace(0, 1, 200) for dx in range(0, 3): assert_allclose(pp(xi, dx), splev(xi, spl, dx)) def test_derivative(self): np.random.seed(1234) x = np.sort(np.r_[0, np.random.rand(11), 1]) y = np.random.rand(len(x)) spl = splrep(x, y, s=0, k=5) pp = PPoly.from_spline(spl) xi = np.linspace(0, 1, 200) for dx in range(0, 10): assert_allclose(pp(xi, dx), pp.derivative(dx)(xi), err_msg="dx=%d" % (dx,)) def test_antiderivative_of_constant(self): # https://github.com/scipy/scipy/issues/4216 p = PPoly([[1.]], [0, 1]) assert_equal(p.antiderivative().c, PPoly([[1], [0]], [0, 1]).c) assert_equal(p.antiderivative().x, PPoly([[1], [0]], [0, 1]).x) def test_antiderivative_regression_4355(self): # https://github.com/scipy/scipy/issues/4355 p = PPoly([[1., 0.5]], [0, 1, 2]) q = p.antiderivative() assert_equal(q.c, [[1, 0.5], [0, 1]]) assert_equal(q.x, [0, 1, 2]) assert_allclose(p.integrate(0, 2), 1.5) assert_allclose(q(2) - q(0), 1.5) def test_antiderivative_simple(self): np.random.seed(1234) # [ p1(x) = 3*x**2 + 2*x + 1, # p2(x) = 1.6875] c = np.array([[3, 2, 1], [0, 0, 1.6875]]).T # [ pp1(x) = x**3 + x**2 + x, # pp2(x) = 1.6875*(x - 0.25) + pp1(0.25)] ic = np.array([[1, 1, 1, 0], [0, 0, 1.6875, 0.328125]]).T # [ ppp1(x) = (1/4)*x**4 + (1/3)*x**3 + (1/2)*x**2, # ppp2(x) = (1.6875/2)*(x - 0.25)**2 + pp1(0.25)*x + ppp1(0.25)] iic = np.array([[1/4, 1/3, 1/2, 0, 0], [0, 0, 1.6875/2, 0.328125, 0.037434895833333336]]).T x = np.array([0, 0.25, 1]) pp = PPoly(c, x) ipp = pp.antiderivative() iipp = pp.antiderivative(2) iipp2 = ipp.antiderivative() assert_allclose(ipp.x, x) assert_allclose(ipp.c.T, ic.T) assert_allclose(iipp.c.T, iic.T) assert_allclose(iipp2.c.T, iic.T) def test_antiderivative_vs_derivative(self): np.random.seed(1234) x = np.linspace(0, 1, 30)**2 y = np.random.rand(len(x)) spl = splrep(x, y, s=0, k=5) pp = PPoly.from_spline(spl) for dx in range(0, 10): ipp = pp.antiderivative(dx) # check that derivative is inverse op pp2 = ipp.derivative(dx) assert_allclose(pp.c, pp2.c) # check continuity for k in range(dx): pp2 = ipp.derivative(k) r = 1e-13 endpoint = r*pp2.x[:-1] + (1 - r)*pp2.x[1:] assert_allclose(pp2(pp2.x[1:]), pp2(endpoint), rtol=1e-7, err_msg="dx=%d k=%d" % (dx, k)) def test_antiderivative_vs_spline(self): np.random.seed(1234) x = np.sort(np.r_[0, np.random.rand(11), 1]) y = np.random.rand(len(x)) spl = splrep(x, y, s=0, k=5) pp = PPoly.from_spline(spl) for dx in range(0, 10): pp2 = pp.antiderivative(dx) spl2 = splantider(spl, dx) xi = np.linspace(0, 1, 200) assert_allclose(pp2(xi), splev(xi, spl2), rtol=1e-7) def test_antiderivative_continuity(self): c = np.array([[2, 1, 2, 2], [2, 1, 3, 3]]).T x = np.array([0, 0.5, 1]) p = PPoly(c, x) ip = p.antiderivative() # check continuity assert_allclose(ip(0.5 - 1e-9), ip(0.5 + 1e-9), rtol=1e-8) # check that only lowest order coefficients were changed p2 = ip.derivative() assert_allclose(p2.c, p.c) def test_integrate(self): np.random.seed(1234) x = np.sort(np.r_[0, np.random.rand(11), 1]) y = np.random.rand(len(x)) spl = splrep(x, y, s=0, k=5) pp = PPoly.from_spline(spl) a, b = 0.3, 0.9 ig = pp.integrate(a, b) ipp = pp.antiderivative() assert_allclose(ig, ipp(b) - ipp(a)) assert_allclose(ig, splint(a, b, spl)) a, b = -0.3, 0.9 ig = pp.integrate(a, b, extrapolate=True) assert_allclose(ig, ipp(b) - ipp(a)) assert_(np.isnan(pp.integrate(a, b, extrapolate=False)).all()) def test_integrate_readonly(self): x = np.array([1, 2, 4]) c = np.array([[0., 0.], [-1., -1.], [2., -0.], [1., 2.]]) for writeable in (True, False): x.flags.writeable = writeable P = PPoly(c, x) vals = P.integrate(1, 4) assert_(np.isfinite(vals).all()) def test_integrate_periodic(self): x = np.array([1, 2, 4]) c = np.array([[0., 0.], [-1., -1.], [2., -0.], [1., 2.]]) P = PPoly(c, x, extrapolate='periodic') I = P.antiderivative() period_int = I(4) - I(1) assert_allclose(P.integrate(1, 4), period_int) assert_allclose(P.integrate(-10, -7), period_int) assert_allclose(P.integrate(-10, -4), 2 * period_int) assert_allclose(P.integrate(1.5, 2.5), I(2.5) - I(1.5)) assert_allclose(P.integrate(3.5, 5), I(2) - I(1) + I(4) - I(3.5)) assert_allclose(P.integrate(3.5 + 12, 5 + 12), I(2) - I(1) + I(4) - I(3.5)) assert_allclose(P.integrate(3.5, 5 + 12), I(2) - I(1) + I(4) - I(3.5) + 4 * period_int) assert_allclose(P.integrate(0, -1), I(2) - I(3)) assert_allclose(P.integrate(-9, -10), I(2) - I(3)) assert_allclose(P.integrate(0, -10), I(2) - I(3) - 3 * period_int) def test_roots(self): x = np.linspace(0, 1, 31)**2 y = np.sin(30*x) spl = splrep(x, y, s=0, k=3) pp = PPoly.from_spline(spl) r = pp.roots() r = r[(r >= 0 - 1e-15) & (r <= 1 + 1e-15)] assert_allclose(r, sproot(spl), atol=1e-15) def test_roots_idzero(self): # Roots for piecewise polynomials with identically zero # sections. c = np.array([[-1, 0.25], [0, 0], [-1, 0.25]]).T x = np.array([0, 0.4, 0.6, 1.0]) pp = PPoly(c, x) assert_array_equal(pp.roots(), [0.25, 0.4, np.nan, 0.6 + 0.25]) # ditto for p.solve(const) with sections identically equal const const = 2. c1 = c.copy() c1[1, :] += const pp1 = PPoly(c1, x) assert_array_equal(pp1.solve(const), [0.25, 0.4, np.nan, 0.6 + 0.25]) def test_roots_all_zero(self): # test the code path for the polynomial being identically zero everywhere c = [[0], [0]] x = [0, 1] p = PPoly(c, x) assert_array_equal(p.roots(), [0, np.nan]) assert_array_equal(p.solve(0), [0, np.nan]) assert_array_equal(p.solve(1), []) c = [[0, 0], [0, 0]] x = [0, 1, 2] p = PPoly(c, x) assert_array_equal(p.roots(), [0, np.nan, 1, np.nan]) assert_array_equal(p.solve(0), [0, np.nan, 1, np.nan]) assert_array_equal(p.solve(1), []) def test_roots_repeated(self): # Check roots repeated in multiple sections are reported only # once. # [(x + 1)**2 - 1, -x**2] ; x == 0 is a repeated root c = np.array([[1, 0, -1], [-1, 0, 0]]).T x = np.array([-1, 0, 1]) pp = PPoly(c, x) assert_array_equal(pp.roots(), [-2, 0]) assert_array_equal(pp.roots(extrapolate=False), [0]) def test_roots_discont(self): # Check that a discontinuity across zero is reported as root c = np.array([[1], [-1]]).T x = np.array([0, 0.5, 1]) pp = PPoly(c, x) assert_array_equal(pp.roots(), [0.5]) assert_array_equal(pp.roots(discontinuity=False), []) # ditto for a discontinuity across y: assert_array_equal(pp.solve(0.5), [0.5]) assert_array_equal(pp.solve(0.5, discontinuity=False), []) assert_array_equal(pp.solve(1.5), []) assert_array_equal(pp.solve(1.5, discontinuity=False), []) def test_roots_random(self): # Check high-order polynomials with random coefficients np.random.seed(1234) num = 0 for extrapolate in (True, False): for order in range(0, 20): x = np.unique(np.r_[0, 10 * np.random.rand(30), 10]) c = 2*np.random.rand(order+1, len(x)-1, 2, 3) - 1 pp = PPoly(c, x) for y in [0, np.random.random()]: r = pp.solve(y, discontinuity=False, extrapolate=extrapolate) for i in range(2): for j in range(3): rr = r[i,j] if rr.size > 0: # Check that the reported roots indeed are roots num += rr.size val = pp(rr, extrapolate=extrapolate)[:,i,j] cmpval = pp(rr, nu=1, extrapolate=extrapolate)[:,i,j] msg = f"({extrapolate!r}) r = {repr(rr)}" assert_allclose((val-y) / cmpval, 0, atol=1e-7, err_msg=msg) # Check that we checked a number of roots assert_(num > 100, repr(num)) def test_roots_croots(self): # Test the complex root finding algorithm np.random.seed(1234) for k in range(1, 15): c = np.random.rand(k, 1, 130) if k == 3: # add a case with zero discriminant c[:,0,0] = 1, 2, 1 for y in [0, np.random.random()]: w = np.empty(c.shape, dtype=complex) _ppoly._croots_poly1(c, w) if k == 1: assert_(np.isnan(w).all()) continue res = 0 cres = 0 for i in range(k): res += c[i,None] * w**(k-1-i) cres += abs(c[i,None] * w**(k-1-i)) with np.errstate(invalid='ignore'): res /= cres res = res.ravel() res = res[~np.isnan(res)] assert_allclose(res, 0, atol=1e-10) def test_extrapolate_attr(self): # [ 1 - x**2 ] c = np.array([[-1, 0, 1]]).T x = np.array([0, 1]) for extrapolate in [True, False, None]: pp = PPoly(c, x, extrapolate=extrapolate) pp_d = pp.derivative() pp_i = pp.antiderivative() if extrapolate is False: assert_(np.isnan(pp([-0.1, 1.1])).all()) assert_(np.isnan(pp_i([-0.1, 1.1])).all()) assert_(np.isnan(pp_d([-0.1, 1.1])).all()) assert_equal(pp.roots(), [1]) else: assert_allclose(pp([-0.1, 1.1]), [1-0.1**2, 1-1.1**2]) assert_(not np.isnan(pp_i([-0.1, 1.1])).any()) assert_(not np.isnan(pp_d([-0.1, 1.1])).any()) assert_allclose(pp.roots(), [1, -1]) class TestBPoly: def test_simple(self): x = [0, 1] c = [[3]] bp = BPoly(c, x) assert_allclose(bp(0.1), 3.) def test_simple2(self): x = [0, 1] c = [[3], [1]] bp = BPoly(c, x) # 3*(1-x) + 1*x assert_allclose(bp(0.1), 3*0.9 + 1.*0.1) def test_simple3(self): x = [0, 1] c = [[3], [1], [4]] bp = BPoly(c, x) # 3 * (1-x)**2 + 2 * x (1-x) + 4 * x**2 assert_allclose(bp(0.2), 3 * 0.8*0.8 + 1 * 2*0.2*0.8 + 4 * 0.2*0.2) def test_simple4(self): x = [0, 1] c = [[1], [1], [1], [2]] bp = BPoly(c, x) assert_allclose(bp(0.3), 0.7**3 + 3 * 0.7**2 * 0.3 + 3 * 0.7 * 0.3**2 + 2 * 0.3**3) def test_simple5(self): x = [0, 1] c = [[1], [1], [8], [2], [1]] bp = BPoly(c, x) assert_allclose(bp(0.3), 0.7**4 + 4 * 0.7**3 * 0.3 + 8 * 6 * 0.7**2 * 0.3**2 + 2 * 4 * 0.7 * 0.3**3 + 0.3**4) def test_periodic(self): x = [0, 1, 3] c = [[3, 0], [0, 0], [0, 2]] # [3*(1-x)**2, 2*((x-1)/2)**2] bp = BPoly(c, x, extrapolate='periodic') assert_allclose(bp(3.4), 3 * 0.6**2) assert_allclose(bp(-1.3), 2 * (0.7/2)**2) assert_allclose(bp(3.4, 1), -6 * 0.6) assert_allclose(bp(-1.3, 1), 2 * (0.7/2)) def test_descending(self): np.random.seed(0) power = 3 for m in [10, 20, 30]: x = np.sort(np.random.uniform(0, 10, m + 1)) ca = np.random.uniform(-0.1, 0.1, size=(power + 1, m)) # We need only to flip coefficients to get it right! cd = ca[::-1].copy() pa = BPoly(ca, x, extrapolate=True) pd = BPoly(cd[:, ::-1], x[::-1], extrapolate=True) x_test = np.random.uniform(-10, 20, 100) assert_allclose(pa(x_test), pd(x_test), rtol=1e-13) assert_allclose(pa(x_test, 1), pd(x_test, 1), rtol=1e-13) pa_d = pa.derivative() pd_d = pd.derivative() assert_allclose(pa_d(x_test), pd_d(x_test), rtol=1e-13) # Antiderivatives won't be equal because fixing continuity is # done in the reverse order, but surely the differences should be # equal. pa_i = pa.antiderivative() pd_i = pd.antiderivative() for a, b in np.random.uniform(-10, 20, (5, 2)): int_a = pa.integrate(a, b) int_d = pd.integrate(a, b) assert_allclose(int_a, int_d, rtol=1e-12) assert_allclose(pa_i(b) - pa_i(a), pd_i(b) - pd_i(a), rtol=1e-12) def test_multi_shape(self): c = np.random.rand(6, 2, 1, 2, 3) x = np.array([0, 0.5, 1]) p = BPoly(c, x) assert_equal(p.x.shape, x.shape) assert_equal(p.c.shape, c.shape) assert_equal(p(0.3).shape, c.shape[2:]) assert_equal(p(np.random.rand(5,6)).shape, (5,6)+c.shape[2:]) dp = p.derivative() assert_equal(dp.c.shape, (5, 2, 1, 2, 3)) def test_interval_length(self): x = [0, 2] c = [[3], [1], [4]] bp = BPoly(c, x) xval = 0.1 s = xval / 2 # s = (x - xa) / (xb - xa) assert_allclose(bp(xval), 3 * (1-s)*(1-s) + 1 * 2*s*(1-s) + 4 * s*s) def test_two_intervals(self): x = [0, 1, 3] c = [[3, 0], [0, 0], [0, 2]] bp = BPoly(c, x) # [3*(1-x)**2, 2*((x-1)/2)**2] assert_allclose(bp(0.4), 3 * 0.6*0.6) assert_allclose(bp(1.7), 2 * (0.7/2)**2) def test_extrapolate_attr(self): x = [0, 2] c = [[3], [1], [4]] bp = BPoly(c, x) for extrapolate in (True, False, None): bp = BPoly(c, x, extrapolate=extrapolate) bp_d = bp.derivative() if extrapolate is False: assert_(np.isnan(bp([-0.1, 2.1])).all()) assert_(np.isnan(bp_d([-0.1, 2.1])).all()) else: assert_(not np.isnan(bp([-0.1, 2.1])).any()) assert_(not np.isnan(bp_d([-0.1, 2.1])).any()) class TestBPolyCalculus: def test_derivative(self): x = [0, 1, 3] c = [[3, 0], [0, 0], [0, 2]] bp = BPoly(c, x) # [3*(1-x)**2, 2*((x-1)/2)**2] bp_der = bp.derivative() 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) 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) pp = PPoly.from_bernstein_basis(bp) for d in range(k): bp = bp.derivative() pp = pp.derivative() 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)) # test both real and complex coefficients for cc in [c.copy(), c*(1. + 2.j)]: bp = BPoly(cc, x) xp = np.linspace(x[0], x[-1], 21) for i in range(k): assert_allclose(bp(xp, i), bp.derivative(i)(xp)) def test_antiderivative_simple(self): # f(x) = x for x \in [0, 1), # (x-1)/2 for x \in [1, 3] # # antiderivative is then # F(x) = x**2 / 2 for x \in [0, 1), # 0.5*x*(x/2 - 1) + A for x \in [1, 3] # where A = 3/4 for continuity at x = 1. x = [0, 1, 3] c = [[0, 0], [1, 1]] bp = BPoly(c, x) bi = bp.antiderivative() xx = np.linspace(0, 3, 11) assert_allclose(bi(xx), np.where(xx < 1, xx**2 / 2., 0.5 * xx * (xx/2. - 1) + 3./4), atol=1e-12, rtol=1e-12) def test_der_antider(self): np.random.seed(1234) x = np.sort(np.random.random(11)) c = np.random.random((4, 10, 2, 3)) bp = BPoly(c, x) xx = np.linspace(x[0], x[-1], 100) assert_allclose(bp.antiderivative().derivative()(xx), bp(xx), atol=1e-12, rtol=1e-12) def test_antider_ppoly(self): np.random.seed(1234) x = np.sort(np.random.random(11)) c = np.random.random((4, 10, 2, 3)) bp = BPoly(c, x) pp = PPoly.from_bernstein_basis(bp) xx = np.linspace(x[0], x[-1], 10) assert_allclose(bp.antiderivative(2)(xx), pp.antiderivative(2)(xx), atol=1e-12, rtol=1e-12) def test_antider_continuous(self): np.random.seed(1234) x = np.sort(np.random.random(11)) c = np.random.random((4, 10)) bp = BPoly(c, x).antiderivative() xx = bp.x[1:-1] assert_allclose(bp(xx - 1e-14), bp(xx + 1e-14), atol=1e-12, rtol=1e-12) def test_integrate(self): np.random.seed(1234) x = np.sort(np.random.random(11)) c = np.random.random((4, 10)) bp = BPoly(c, x) pp = PPoly.from_bernstein_basis(bp) assert_allclose(bp.integrate(0, 1), pp.integrate(0, 1), atol=1e-12, rtol=1e-12) def test_integrate_extrap(self): c = [[1]] x = [0, 1] b = BPoly(c, x) # default is extrapolate=True assert_allclose(b.integrate(0, 2), 2., atol=1e-14) # .integrate argument overrides self.extrapolate b1 = BPoly(c, x, extrapolate=False) assert_(np.isnan(b1.integrate(0, 2))) assert_allclose(b1.integrate(0, 2, extrapolate=True), 2., atol=1e-14) def test_integrate_periodic(self): x = np.array([1, 2, 4]) c = np.array([[0., 0.], [-1., -1.], [2., -0.], [1., 2.]]) P = BPoly.from_power_basis(PPoly(c, x), extrapolate='periodic') I = P.antiderivative() period_int = I(4) - I(1) assert_allclose(P.integrate(1, 4), period_int) assert_allclose(P.integrate(-10, -7), period_int) assert_allclose(P.integrate(-10, -4), 2 * period_int) assert_allclose(P.integrate(1.5, 2.5), I(2.5) - I(1.5)) assert_allclose(P.integrate(3.5, 5), I(2) - I(1) + I(4) - I(3.5)) assert_allclose(P.integrate(3.5 + 12, 5 + 12), I(2) - I(1) + I(4) - I(3.5)) assert_allclose(P.integrate(3.5, 5 + 12), I(2) - I(1) + I(4) - I(3.5) + 4 * period_int) assert_allclose(P.integrate(0, -1), I(2) - I(3)) assert_allclose(P.integrate(-9, -10), I(2) - I(3)) assert_allclose(P.integrate(0, -10), I(2) - I(3) - 3 * period_int) def test_antider_neg(self): # .derivative(-nu) ==> .andiderivative(nu) and vice versa c = [[1]] x = [0, 1] b = BPoly(c, x) xx = np.linspace(0, 1, 21) assert_allclose(b.derivative(-1)(xx), b.antiderivative()(xx), atol=1e-12, rtol=1e-12) assert_allclose(b.derivative(1)(xx), b.antiderivative(-1)(xx), atol=1e-12, rtol=1e-12) class TestPolyConversions: def test_bp_from_pp(self): x = [0, 1, 3] c = [[3, 2], [1, 8], [4, 3]] pp = PPoly(c, x) bp = BPoly.from_power_basis(pp) pp1 = PPoly.from_bernstein_basis(bp) xp = [0.1, 1.4] assert_allclose(pp(xp), bp(xp)) assert_allclose(pp(xp), pp1(xp)) def test_bp_from_pp_random(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)) pp = PPoly(c, x) bp = BPoly.from_power_basis(pp) pp1 = PPoly.from_bernstein_basis(bp) xp = np.linspace(x[0], x[-1], 21) assert_allclose(pp(xp), bp(xp)) assert_allclose(pp(xp), pp1(xp)) def test_pp_from_bp(self): x = [0, 1, 3] c = [[3, 3], [1, 1], [4, 2]] bp = BPoly(c, x) pp = PPoly.from_bernstein_basis(bp) bp1 = BPoly.from_power_basis(pp) xp = [0.1, 1.4] assert_allclose(bp(xp), pp(xp)) assert_allclose(bp(xp), bp1(xp)) def test_broken_conversions(self): # regression test for gh-10597: from_power_basis only accepts PPoly etc. x = [0, 1, 3] c = [[3, 3], [1, 1], [4, 2]] pp = PPoly(c, x) with assert_raises(TypeError): PPoly.from_bernstein_basis(pp) bp = BPoly(c, x) with assert_raises(TypeError): BPoly.from_power_basis(bp) class TestBPolyFromDerivatives: def test_make_poly_1(self): c1 = BPoly._construct_from_derivatives(0, 1, [2], [3]) assert_allclose(c1, [2., 3.]) def test_make_poly_2(self): c1 = BPoly._construct_from_derivatives(0, 1, [1, 0], [1]) assert_allclose(c1, [1., 1., 1.]) # f'(0) = 3 c2 = BPoly._construct_from_derivatives(0, 1, [2, 3], [1]) assert_allclose(c2, [2., 7./2, 1.]) # f'(1) = 3 c3 = BPoly._construct_from_derivatives(0, 1, [2], [1, 3]) assert_allclose(c3, [2., -0.5, 1.]) def test_make_poly_3(self): # f'(0)=2, f''(0)=3 c1 = BPoly._construct_from_derivatives(0, 1, [1, 2, 3], [4]) assert_allclose(c1, [1., 5./3, 17./6, 4.]) # f'(1)=2, f''(1)=3 c2 = BPoly._construct_from_derivatives(0, 1, [1], [4, 2, 3]) assert_allclose(c2, [1., 19./6, 10./3, 4.]) # f'(0)=2, f'(1)=3 c3 = BPoly._construct_from_derivatives(0, 1, [1, 2], [4, 3]) assert_allclose(c3, [1., 5./3, 3., 4.]) def test_make_poly_12(self): np.random.seed(12345) ya = np.r_[0, np.random.random(5)] yb = np.r_[0, np.random.random(5)] c = BPoly._construct_from_derivatives(0, 1, ya, yb) pp = BPoly(c[:, None], [0, 1]) for j in range(6): assert_allclose([pp(0.), pp(1.)], [ya[j], yb[j]]) pp = pp.derivative() def test_raise_degree(self): np.random.seed(12345) x = [0, 1] k, d = 8, 5 c = np.random.random((k, 1, 2, 3, 4)) bp = BPoly(c, x) c1 = BPoly._raise_degree(c, d) bp1 = BPoly(c1, x) xp = np.linspace(0, 1, 11) assert_allclose(bp(xp), bp1(xp)) def test_xi_yi(self): assert_raises(ValueError, BPoly.from_derivatives, [0, 1], [0]) def test_coords_order(self): xi = [0, 0, 1] yi = [[0], [0], [0]] assert_raises(ValueError, BPoly.from_derivatives, xi, yi) def test_zeros(self): xi = [0, 1, 2, 3] yi = [[0, 0], [0], [0, 0], [0, 0]] # NB: will have to raise the degree pp = BPoly.from_derivatives(xi, yi) assert_(pp.c.shape == (4, 3)) ppd = pp.derivative() for xp in [0., 0.1, 1., 1.1, 1.9, 2., 2.5]: assert_allclose([pp(xp), ppd(xp)], [0., 0.]) def _make_random_mk(self, m, k): # k derivatives at each breakpoint np.random.seed(1234) xi = np.asarray([1. * j**2 for j in range(m+1)]) yi = [np.random.random(k) for j in range(m+1)] return xi, yi def test_random_12(self): m, k = 5, 12 xi, yi = self._make_random_mk(m, k) pp = BPoly.from_derivatives(xi, yi) for order in range(k//2): assert_allclose(pp(xi), [yy[order] for yy in yi]) pp = pp.derivative() def test_order_zero(self): m, k = 5, 12 xi, yi = self._make_random_mk(m, k) assert_raises(ValueError, BPoly.from_derivatives, **dict(xi=xi, yi=yi, orders=0)) def test_orders_too_high(self): m, k = 5, 12 xi, yi = self._make_random_mk(m, k) BPoly.from_derivatives(xi, yi, orders=2*k-1) # this is still ok assert_raises(ValueError, BPoly.from_derivatives, # but this is not **dict(xi=xi, yi=yi, orders=2*k)) def test_orders_global(self): m, k = 5, 12 xi, yi = self._make_random_mk(m, k) # ok, this is confusing. Local polynomials will be of the order 5 # which means that up to the 2nd derivatives will be used at each point order = 5 pp = BPoly.from_derivatives(xi, yi, orders=order) for j in range(order//2+1): assert_allclose(pp(xi[1:-1] - 1e-12), pp(xi[1:-1] + 1e-12)) pp = pp.derivative() assert_(not np.allclose(pp(xi[1:-1] - 1e-12), pp(xi[1:-1] + 1e-12))) # now repeat with `order` being even: on each interval, it uses # order//2 'derivatives' @ the right-hand endpoint and # order//2+1 @ 'derivatives' the left-hand endpoint order = 6 pp = BPoly.from_derivatives(xi, yi, orders=order) for j in range(order//2): assert_allclose(pp(xi[1:-1] - 1e-12), pp(xi[1:-1] + 1e-12)) pp = pp.derivative() assert_(not np.allclose(pp(xi[1:-1] - 1e-12), pp(xi[1:-1] + 1e-12))) def test_orders_local(self): m, k = 7, 12 xi, yi = self._make_random_mk(m, k) orders = [o + 1 for o in range(m)] for i, x in enumerate(xi[1:-1]): pp = BPoly.from_derivatives(xi, yi, orders=orders) for j in range(orders[i] // 2 + 1): assert_allclose(pp(x - 1e-12), pp(x + 1e-12)) pp = pp.derivative() assert_(not np.allclose(pp(x - 1e-12), pp(x + 1e-12))) def test_yi_trailing_dims(self): m, k = 7, 5 xi = np.sort(np.random.random(m+1)) yi = np.random.random((m+1, k, 6, 7, 8)) pp = BPoly.from_derivatives(xi, yi) assert_equal(pp.c.shape, (2*k, m, 6, 7, 8)) def test_gh_5430(self): # At least one of these raises an error unless gh-5430 is # fixed. In py2k an int is implemented using a C long, so # which one fails depends on your system. In py3k there is only # one arbitrary precision integer type, so both should fail. orders = np.int32(1) p = BPoly.from_derivatives([0, 1], [[0], [0]], orders=orders) assert_almost_equal(p(0), 0) orders = np.int64(1) p = BPoly.from_derivatives([0, 1], [[0], [0]], orders=orders) assert_almost_equal(p(0), 0) orders = 1 # This worked before; make sure it still works p = BPoly.from_derivatives([0, 1], [[0], [0]], orders=orders) assert_almost_equal(p(0), 0) orders = 1 class TestNdPPoly: def test_simple_1d(self): np.random.seed(1234) c = np.random.rand(4, 5) x = np.linspace(0, 1, 5+1) xi = np.random.rand(200) p = NdPPoly(c, (x,)) v1 = p((xi,)) v2 = _ppoly_eval_1(c[:,:,None], x, xi).ravel() assert_allclose(v1, v2) def test_simple_2d(self): np.random.seed(1234) c = np.random.rand(4, 5, 6, 7) x = np.linspace(0, 1, 6+1) y = np.linspace(0, 1, 7+1)**2 xi = np.random.rand(200) yi = np.random.rand(200) v1 = np.empty([len(xi), 1], dtype=c.dtype) v1.fill(np.nan) _ppoly.evaluate_nd(c.reshape(4*5, 6*7, 1), (x, y), np.array([4, 5], dtype=np.intc), np.c_[xi, yi], np.array([0, 0], dtype=np.intc), 1, v1) v1 = v1.ravel() v2 = _ppoly2d_eval(c, (x, y), xi, yi) assert_allclose(v1, v2) p = NdPPoly(c, (x, y)) for nu in (None, (0, 0), (0, 1), (1, 0), (2, 3), (9, 2)): v1 = p(np.c_[xi, yi], nu=nu) v2 = _ppoly2d_eval(c, (x, y), xi, yi, nu=nu) assert_allclose(v1, v2, err_msg=repr(nu)) def test_simple_3d(self): np.random.seed(1234) c = np.random.rand(4, 5, 6, 7, 8, 9) x = np.linspace(0, 1, 7+1) y = np.linspace(0, 1, 8+1)**2 z = np.linspace(0, 1, 9+1)**3 xi = np.random.rand(40) yi = np.random.rand(40) zi = np.random.rand(40) p = NdPPoly(c, (x, y, z)) for nu in (None, (0, 0, 0), (0, 1, 0), (1, 0, 0), (2, 3, 0), (6, 0, 2)): v1 = p((xi, yi, zi), nu=nu) v2 = _ppoly3d_eval(c, (x, y, z), xi, yi, zi, nu=nu) assert_allclose(v1, v2, err_msg=repr(nu)) def test_simple_4d(self): np.random.seed(1234) c = np.random.rand(4, 5, 6, 7, 8, 9, 10, 11) x = np.linspace(0, 1, 8+1) y = np.linspace(0, 1, 9+1)**2 z = np.linspace(0, 1, 10+1)**3 u = np.linspace(0, 1, 11+1)**4 xi = np.random.rand(20) yi = np.random.rand(20) zi = np.random.rand(20) ui = np.random.rand(20) p = NdPPoly(c, (x, y, z, u)) v1 = p((xi, yi, zi, ui)) v2 = _ppoly4d_eval(c, (x, y, z, u), xi, yi, zi, ui) assert_allclose(v1, v2) def test_deriv_1d(self): np.random.seed(1234) c = np.random.rand(4, 5) x = np.linspace(0, 1, 5+1) p = NdPPoly(c, (x,)) # derivative dp = p.derivative(nu=[1]) p1 = PPoly(c, x) dp1 = p1.derivative() assert_allclose(dp.c, dp1.c) # antiderivative dp = p.antiderivative(nu=[2]) p1 = PPoly(c, x) dp1 = p1.antiderivative(2) assert_allclose(dp.c, dp1.c) def test_deriv_3d(self): np.random.seed(1234) c = np.random.rand(4, 5, 6, 7, 8, 9) x = np.linspace(0, 1, 7+1) y = np.linspace(0, 1, 8+1)**2 z = np.linspace(0, 1, 9+1)**3 p = NdPPoly(c, (x, y, z)) # differentiate vs x p1 = PPoly(c.transpose(0, 3, 1, 2, 4, 5), x) dp = p.derivative(nu=[2]) dp1 = p1.derivative(2) assert_allclose(dp.c, dp1.c.transpose(0, 2, 3, 1, 4, 5)) # antidifferentiate vs y p1 = PPoly(c.transpose(1, 4, 0, 2, 3, 5), y) dp = p.antiderivative(nu=[0, 1, 0]) dp1 = p1.antiderivative(1) assert_allclose(dp.c, dp1.c.transpose(2, 0, 3, 4, 1, 5)) # differentiate vs z p1 = PPoly(c.transpose(2, 5, 0, 1, 3, 4), z) dp = p.derivative(nu=[0, 0, 3]) dp1 = p1.derivative(3) assert_allclose(dp.c, dp1.c.transpose(2, 3, 0, 4, 5, 1)) def test_deriv_3d_simple(self): # Integrate to obtain function x y**2 z**4 / (2! 4!) c = np.ones((1, 1, 1, 3, 4, 5)) x = np.linspace(0, 1, 3+1)**1 y = np.linspace(0, 1, 4+1)**2 z = np.linspace(0, 1, 5+1)**3 p = NdPPoly(c, (x, y, z)) ip = p.antiderivative((1, 0, 4)) ip = ip.antiderivative((0, 2, 0)) xi = np.random.rand(20) yi = np.random.rand(20) zi = np.random.rand(20) assert_allclose(ip((xi, yi, zi)), xi * yi**2 * zi**4 / (gamma(3)*gamma(5))) def test_integrate_2d(self): np.random.seed(1234) c = np.random.rand(4, 5, 16, 17) x = np.linspace(0, 1, 16+1)**1 y = np.linspace(0, 1, 17+1)**2 # make continuously differentiable so that nquad() has an # easier time c = c.transpose(0, 2, 1, 3) cx = c.reshape(c.shape[0], c.shape[1], -1).copy() _ppoly.fix_continuity(cx, x, 2) c = cx.reshape(c.shape) c = c.transpose(0, 2, 1, 3) c = c.transpose(1, 3, 0, 2) cx = c.reshape(c.shape[0], c.shape[1], -1).copy() _ppoly.fix_continuity(cx, y, 2) c = cx.reshape(c.shape) c = c.transpose(2, 0, 3, 1).copy() # Check integration p = NdPPoly(c, (x, y)) for ranges in [[(0, 1), (0, 1)], [(0, 0.5), (0, 1)], [(0, 1), (0, 0.5)], [(0.3, 0.7), (0.6, 0.2)]]: ig = p.integrate(ranges) ig2, err2 = nquad(lambda x, y: p((x, y)), ranges, opts=[dict(epsrel=1e-5, epsabs=1e-5)]*2) assert_allclose(ig, ig2, rtol=1e-5, atol=1e-5, err_msg=repr(ranges)) def test_integrate_1d(self): np.random.seed(1234) c = np.random.rand(4, 5, 6, 16, 17, 18) x = np.linspace(0, 1, 16+1)**1 y = np.linspace(0, 1, 17+1)**2 z = np.linspace(0, 1, 18+1)**3 # Check 1-D integration p = NdPPoly(c, (x, y, z)) u = np.random.rand(200) v = np.random.rand(200) a, b = 0.2, 0.7 px = p.integrate_1d(a, b, axis=0) pax = p.antiderivative((1, 0, 0)) assert_allclose(px((u, v)), pax((b, u, v)) - pax((a, u, v))) py = p.integrate_1d(a, b, axis=1) pay = p.antiderivative((0, 1, 0)) assert_allclose(py((u, v)), pay((u, b, v)) - pay((u, a, v))) pz = p.integrate_1d(a, b, axis=2) paz = p.antiderivative((0, 0, 1)) assert_allclose(pz((u, v)), paz((u, v, b)) - paz((u, v, a))) def _ppoly_eval_1(c, x, xps): """Evaluate piecewise polynomial manually""" out = np.zeros((len(xps), c.shape[2])) for i, xp in enumerate(xps): if xp < 0 or xp > 1: out[i,:] = np.nan continue j = np.searchsorted(x, xp) - 1 d = xp - x[j] assert_(x[j] <= xp < x[j+1]) r = sum(c[k,j] * d**(c.shape[0]-k-1) for k in range(c.shape[0])) out[i,:] = r return out def _ppoly_eval_2(coeffs, breaks, xnew, fill=np.nan): """Evaluate piecewise polynomial manually (another way)""" a = breaks[0] b = breaks[-1] K = coeffs.shape[0] saveshape = np.shape(xnew) xnew = np.ravel(xnew) res = np.empty_like(xnew) mask = (xnew >= a) & (xnew <= b) res[~mask] = fill xx = xnew.compress(mask) indxs = np.searchsorted(breaks, xx)-1 indxs = indxs.clip(0, len(breaks)) pp = coeffs diff = xx - breaks.take(indxs) V = np.vander(diff, N=K) values = np.array([np.dot(V[k, :], pp[:, indxs[k]]) for k in range(len(xx))]) res[mask] = values res.shape = saveshape return res def _dpow(x, y, n): """ d^n (x**y) / dx^n """ if n < 0: raise ValueError("invalid derivative order") elif n > y: return 0 else: return poch(y - n + 1, n) * x**(y - n) def _ppoly2d_eval(c, xs, xnew, ynew, nu=None): """ Straightforward evaluation of 2-D piecewise polynomial """ if nu is None: nu = (0, 0) out = np.empty((len(xnew),), dtype=c.dtype) nx, ny = c.shape[:2] for jout, (x, y) in enumerate(zip(xnew, ynew)): if not ((xs[0][0] <= x <= xs[0][-1]) and (xs[1][0] <= y <= xs[1][-1])): out[jout] = np.nan continue j1 = np.searchsorted(xs[0], x) - 1 j2 = np.searchsorted(xs[1], y) - 1 s1 = x - xs[0][j1] s2 = y - xs[1][j2] val = 0 for k1 in range(c.shape[0]): for k2 in range(c.shape[1]): val += (c[nx-k1-1,ny-k2-1,j1,j2] * _dpow(s1, k1, nu[0]) * _dpow(s2, k2, nu[1])) out[jout] = val return out def _ppoly3d_eval(c, xs, xnew, ynew, znew, nu=None): """ Straightforward evaluation of 3-D piecewise polynomial """ if nu is None: nu = (0, 0, 0) out = np.empty((len(xnew),), dtype=c.dtype) nx, ny, nz = c.shape[:3] for jout, (x, y, z) in enumerate(zip(xnew, ynew, znew)): if not ((xs[0][0] <= x <= xs[0][-1]) and (xs[1][0] <= y <= xs[1][-1]) and (xs[2][0] <= z <= xs[2][-1])): out[jout] = np.nan continue j1 = np.searchsorted(xs[0], x) - 1 j2 = np.searchsorted(xs[1], y) - 1 j3 = np.searchsorted(xs[2], z) - 1 s1 = x - xs[0][j1] s2 = y - xs[1][j2] s3 = z - xs[2][j3] val = 0 for k1 in range(c.shape[0]): for k2 in range(c.shape[1]): for k3 in range(c.shape[2]): val += (c[nx-k1-1,ny-k2-1,nz-k3-1,j1,j2,j3] * _dpow(s1, k1, nu[0]) * _dpow(s2, k2, nu[1]) * _dpow(s3, k3, nu[2])) out[jout] = val return out def _ppoly4d_eval(c, xs, xnew, ynew, znew, unew, nu=None): """ Straightforward evaluation of 4-D piecewise polynomial """ if nu is None: nu = (0, 0, 0, 0) out = np.empty((len(xnew),), dtype=c.dtype) mx, my, mz, mu = c.shape[:4] for jout, (x, y, z, u) in enumerate(zip(xnew, ynew, znew, unew)): if not ((xs[0][0] <= x <= xs[0][-1]) and (xs[1][0] <= y <= xs[1][-1]) and (xs[2][0] <= z <= xs[2][-1]) and (xs[3][0] <= u <= xs[3][-1])): out[jout] = np.nan continue j1 = np.searchsorted(xs[0], x) - 1 j2 = np.searchsorted(xs[1], y) - 1 j3 = np.searchsorted(xs[2], z) - 1 j4 = np.searchsorted(xs[3], u) - 1 s1 = x - xs[0][j1] s2 = y - xs[1][j2] s3 = z - xs[2][j3] s4 = u - xs[3][j4] val = 0 for k1 in range(c.shape[0]): for k2 in range(c.shape[1]): for k3 in range(c.shape[2]): for k4 in range(c.shape[3]): val += (c[mx-k1-1,my-k2-1,mz-k3-1,mu-k4-1,j1,j2,j3,j4] * _dpow(s1, k1, nu[0]) * _dpow(s2, k2, nu[1]) * _dpow(s3, k3, nu[2]) * _dpow(s4, k4, nu[3])) out[jout] = val return out