import numpy as np from numpy.testing import assert_allclose, assert_equal import pytest from scipy.optimize._pava_pybind import pava from scipy.optimize import isotonic_regression class TestIsotonicRegression: @pytest.mark.parametrize( ("y", "w", "msg"), [ ([[0, 1]], None, "array has incorrect number of dimensions: 2; expected 1"), ([0, 1], [[1, 2]], "Input arrays y and w must have one dimension of equal length"), ([0, 1], [1], "Input arrays y and w must have one dimension of equal length"), (1, [1, 2], "Input arrays y and w must have one dimension of equal length"), ([1, 2], 1, "Input arrays y and w must have one dimension of equal length"), ([0, 1], [0, 1], "Weights w must be strictly positive"), ] ) def test_raise_error(self, y, w, msg): with pytest.raises(ValueError, match=msg): isotonic_regression(y=y, weights=w) def test_simple_pava(self): # Test case of Busing 2020 # https://doi.org/10.18637/jss.v102.c01 y = np.array([8, 4, 8, 2, 2, 0, 8], dtype=np.float64) w = np.ones_like(y) r = np.full(shape=y.shape[0] + 1, fill_value=-1, dtype=np.intp) pava(y, w, r) assert_allclose(y, [4, 4, 4, 4, 4, 4, 8]) # Only first 2 elements of w are changed. assert_allclose(w, [6, 1, 1, 1, 1, 1, 1]) # Only first 3 elements of r are changed. assert_allclose(r, [0, 6, 7, -1, -1, -1, -1, -1]) @pytest.mark.parametrize("y_dtype", [np.float64, np.float32, np.int64, np.int32]) @pytest.mark.parametrize("w_dtype", [np.float64, np.float32, np.int64, np.int32]) @pytest.mark.parametrize("w", [None, "ones"]) def test_simple_isotonic_regression(self, w, w_dtype, y_dtype): # Test case of Busing 2020 # https://doi.org/10.18637/jss.v102.c01 y = np.array([8, 4, 8, 2, 2, 0, 8], dtype=y_dtype) if w is not None: w = np.ones_like(y, dtype=w_dtype) res = isotonic_regression(y, weights=w) assert res.x.dtype == np.float64 assert res.weights.dtype == np.float64 assert_allclose(res.x, [4, 4, 4, 4, 4, 4, 8]) assert_allclose(res.weights, [6, 1]) assert_allclose(res.blocks, [0, 6, 7]) # Assert that y was not overwritten assert_equal(y, np.array([8, 4, 8, 2, 2, 0, 8], dtype=np.float64)) @pytest.mark.parametrize("increasing", [True, False]) def test_linspace(self, increasing): n = 10 y = np.linspace(0, 1, n) if increasing else np.linspace(1, 0, n) res = isotonic_regression(y, increasing=increasing) assert_allclose(res.x, y) assert_allclose(res.blocks, np.arange(n + 1)) def test_weights(self): w = np.array([1, 2, 5, 0.5, 0.5, 0.5, 1, 3]) y = np.array([3, 2, 1, 10, 9, 8, 20, 10]) res = isotonic_regression(y, weights=w) assert_allclose(res.x, [12/8, 12/8, 12/8, 9, 9, 9, 50/4, 50/4]) assert_allclose(res.weights, [8, 1.5, 4]) assert_allclose(res.blocks, [0, 3, 6, 8]) # weights are like repeated observations, we repeat the 3rd element 5 # times. w2 = np.array([1, 2, 1, 1, 1, 1, 1, 0.5, 0.5, 0.5, 1, 3]) y2 = np.array([3, 2, 1, 1, 1, 1, 1, 10, 9, 8, 20, 10]) res2 = isotonic_regression(y2, weights=w2) assert_allclose(np.diff(res2.x[0:7]), 0) assert_allclose(res2.x[4:], res.x) assert_allclose(res2.weights, res.weights) assert_allclose(res2.blocks[1:] - 4, res.blocks[1:]) def test_against_R_monotone(self): y = [0, 6, 8, 3, 5, 2, 1, 7, 9, 4] res = isotonic_regression(y) # R code # library(monotone) # options(digits=8) # monotone(c(0, 6, 8, 3, 5, 2, 1, 7, 9, 4)) x_R = [ 0, 4.1666667, 4.1666667, 4.1666667, 4.1666667, 4.1666667, 4.1666667, 6.6666667, 6.6666667, 6.6666667, ] assert_allclose(res.x, x_R) assert_equal(res.blocks, [0, 1, 7, 10]) n = 100 y = np.linspace(0, 1, num=n, endpoint=False) y = 5 * y + np.sin(10 * y) res = isotonic_regression(y) # R code # library(monotone) # n <- 100 # y <- 5 * ((1:n)-1)/n + sin(10 * ((1:n)-1)/n) # options(digits=8) # monotone(y) x_R = [ 0.00000000, 0.14983342, 0.29866933, 0.44552021, 0.58941834, 0.72942554, 0.86464247, 0.99421769, 1.11735609, 1.23332691, 1.34147098, 1.44120736, 1.53203909, 1.57081100, 1.57081100, 1.57081100, 1.57081100, 1.57081100, 1.57081100, 1.57081100, 1.57081100, 1.57081100, 1.57081100, 1.57081100, 1.57081100, 1.57081100, 1.57081100, 1.57081100, 1.57081100, 1.57081100, 1.57081100, 1.57081100, 1.57081100, 1.57081100, 1.57081100, 1.57081100, 1.57081100, 1.57081100, 1.57081100, 1.57081100, 1.57081100, 1.57081100, 1.57081100, 1.57081100, 1.57081100, 1.57081100, 1.57081100, 1.57081100, 1.57081100, 1.57081100, 1.57081100, 1.62418532, 1.71654534, 1.81773256, 1.92723551, 2.04445967, 2.16873336, 2.29931446, 2.43539782, 2.57612334, 2.72058450, 2.86783750, 3.01691060, 3.16681390, 3.31654920, 3.46511999, 3.61154136, 3.75484992, 3.89411335, 4.02843976, 4.15698660, 4.27896904, 4.39366786, 4.50043662, 4.59870810, 4.68799998, 4.76791967, 4.83816823, 4.86564130, 4.86564130, 4.86564130, 4.86564130, 4.86564130, 4.86564130, 4.86564130, 4.86564130, 4.86564130, 4.86564130, 4.86564130, 4.86564130, 4.86564130, 4.86564130, 4.86564130, 4.86564130, 4.86564130, 4.86564130, 4.86564130, 4.86564130, 4.86564130, 4.86564130, ] assert_allclose(res.x, x_R) # Test increasing assert np.all(np.diff(res.x) >= 0) # Test balance property: sum(y) == sum(x) assert_allclose(np.sum(res.x), np.sum(y)) # Reverse order res_inv = isotonic_regression(-y, increasing=False) assert_allclose(-res_inv.x, res.x) assert_equal(res_inv.blocks, res.blocks) def test_readonly(self): x = np.arange(3, dtype=float) w = np.ones(3, dtype=float) x.flags.writeable = False w.flags.writeable = False res = isotonic_regression(x, weights=w) assert np.all(np.isfinite(res.x)) assert np.all(np.isfinite(res.weights)) assert np.all(np.isfinite(res.blocks)) def test_non_contiguous_arrays(self): x = np.arange(10, dtype=float)[::3] w = np.ones(10, dtype=float)[::3] assert not x.flags.c_contiguous assert not x.flags.f_contiguous assert not w.flags.c_contiguous assert not w.flags.f_contiguous res = isotonic_regression(x, weights=w) assert np.all(np.isfinite(res.x)) assert np.all(np.isfinite(res.weights)) assert np.all(np.isfinite(res.blocks))