image-match/python/py/Lib/site-packages/scipy/optimize/tests/test_isotonic_regression.py

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))