Vehicle-cpp/hkpc/software/scipy/optimize/_trustregion_constr/tests/test_projections.py

import numpy as np
import scipy.linalg
from scipy.sparse import csc_matrix
from scipy.optimize._trustregion_constr.projections \
    import projections, orthogonality
from numpy.testing import (TestCase, assert_array_almost_equal,
                           assert_equal, assert_allclose)

try:
    from sksparse.cholmod import cholesky_AAt
    sksparse_available = True
    available_sparse_methods = ("NormalEquation", "AugmentedSystem")
except ImportError:
    sksparse_available = False
    available_sparse_methods = ("AugmentedSystem",)
available_dense_methods = ('QRFactorization', 'SVDFactorization')


class TestProjections(TestCase):

    def test_nullspace_and_least_squares_sparse(self):
        A_dense = np.array([[1, 2, 3, 4, 0, 5, 0, 7],
                            [0, 8, 7, 0, 1, 5, 9, 0],
                            [1, 0, 0, 0, 0, 1, 2, 3]])
        At_dense = A_dense.T
        A = csc_matrix(A_dense)
        test_points = ([1, 2, 3, 4, 5, 6, 7, 8],
                       [1, 10, 3, 0, 1, 6, 7, 8],
                       [1.12, 10, 0, 0, 100000, 6, 0.7, 8])

        for method in available_sparse_methods:
            Z, LS, _ = projections(A, method)
            for z in test_points:
                # Test if x is in the null_space
                x = Z.matvec(z)
                assert_array_almost_equal(A.dot(x), 0)
                # Test orthogonality
                assert_array_almost_equal(orthogonality(A, x), 0)
                # Test if x is the least square solution
                x = LS.matvec(z)
                x2 = scipy.linalg.lstsq(At_dense, z)[0]
                assert_array_almost_equal(x, x2)

    def test_iterative_refinements_sparse(self):
        A_dense = np.array([[1, 2, 3, 4, 0, 5, 0, 7],
                            [0, 8, 7, 0, 1, 5, 9, 0],
                            [1, 0, 0, 0, 0, 1, 2, 3]])
        A = csc_matrix(A_dense)
        test_points = ([1, 2, 3, 4, 5, 6, 7, 8],
                       [1, 10, 3, 0, 1, 6, 7, 8],
                       [1.12, 10, 0, 0, 100000, 6, 0.7, 8],
                       [1, 0, 0, 0, 0, 1, 2, 3+1e-10])

        for method in available_sparse_methods:
            Z, LS, _ = projections(A, method, orth_tol=1e-18, max_refin=100)
            for z in test_points:
                # Test if x is in the null_space
                x = Z.matvec(z)
                atol = 1e-13 * abs(x).max()
                assert_allclose(A.dot(x), 0, atol=atol)
                # Test orthogonality
                assert_allclose(orthogonality(A, x), 0, atol=1e-13)

    def test_rowspace_sparse(self):
        A_dense = np.array([[1, 2, 3, 4, 0, 5, 0, 7],
                            [0, 8, 7, 0, 1, 5, 9, 0],
                            [1, 0, 0, 0, 0, 1, 2, 3]])
        A = csc_matrix(A_dense)
        test_points = ([1, 2, 3],
                       [1, 10, 3],
                       [1.12, 10, 0])

        for method in available_sparse_methods:
            _, _, Y = projections(A, method)
            for z in test_points:
                # Test if x is solution of A x = z
                x = Y.matvec(z)
                assert_array_almost_equal(A.dot(x), z)
                # Test if x is in the return row space of A
                A_ext = np.vstack((A_dense, x))
                assert_equal(np.linalg.matrix_rank(A_dense),
                             np.linalg.matrix_rank(A_ext))

    def test_nullspace_and_least_squares_dense(self):
        A = np.array([[1, 2, 3, 4, 0, 5, 0, 7],
                      [0, 8, 7, 0, 1, 5, 9, 0],
                      [1, 0, 0, 0, 0, 1, 2, 3]])
        At = A.T
        test_points = ([1, 2, 3, 4, 5, 6, 7, 8],
                       [1, 10, 3, 0, 1, 6, 7, 8],
                       [1.12, 10, 0, 0, 100000, 6, 0.7, 8])

        for method in available_dense_methods:
            Z, LS, _ = projections(A, method)
            for z in test_points:
                # Test if x is in the null_space
                x = Z.matvec(z)
                assert_array_almost_equal(A.dot(x), 0)
                # Test orthogonality
                assert_array_almost_equal(orthogonality(A, x), 0)
                # Test if x is the least square solution
                x = LS.matvec(z)
                x2 = scipy.linalg.lstsq(At, z)[0]
                assert_array_almost_equal(x, x2)

    def test_compare_dense_and_sparse(self):
        D = np.diag(range(1, 101))
        A = np.hstack([D, D, D, D])
        A_sparse = csc_matrix(A)
        np.random.seed(0)

        Z, LS, Y = projections(A)
        Z_sparse, LS_sparse, Y_sparse = projections(A_sparse)
        for k in range(20):
            z = np.random.normal(size=(400,))
            assert_array_almost_equal(Z.dot(z), Z_sparse.dot(z))
            assert_array_almost_equal(LS.dot(z), LS_sparse.dot(z))
            x = np.random.normal(size=(100,))
            assert_array_almost_equal(Y.dot(x), Y_sparse.dot(x))

    def test_compare_dense_and_sparse2(self):
        D1 = np.diag([-1.7, 1, 0.5])
        D2 = np.diag([1, -0.6, -0.3])
        D3 = np.diag([-0.3, -1.5, 2])
        A = np.hstack([D1, D2, D3])
        A_sparse = csc_matrix(A)
        np.random.seed(0)

        Z, LS, Y = projections(A)
        Z_sparse, LS_sparse, Y_sparse = projections(A_sparse)
        for k in range(1):
            z = np.random.normal(size=(9,))
            assert_array_almost_equal(Z.dot(z), Z_sparse.dot(z))
            assert_array_almost_equal(LS.dot(z), LS_sparse.dot(z))
            x = np.random.normal(size=(3,))
            assert_array_almost_equal(Y.dot(x), Y_sparse.dot(x))

    def test_iterative_refinements_dense(self):
        A = np.array([[1, 2, 3, 4, 0, 5, 0, 7],
                            [0, 8, 7, 0, 1, 5, 9, 0],
                            [1, 0, 0, 0, 0, 1, 2, 3]])
        test_points = ([1, 2, 3, 4, 5, 6, 7, 8],
                       [1, 10, 3, 0, 1, 6, 7, 8],
                       [1, 0, 0, 0, 0, 1, 2, 3+1e-10])

        for method in available_dense_methods:
            Z, LS, _ = projections(A, method, orth_tol=1e-18, max_refin=10)
            for z in test_points:
                # Test if x is in the null_space
                x = Z.matvec(z)
                assert_allclose(A.dot(x), 0, rtol=0, atol=2.5e-14)
                # Test orthogonality
                assert_allclose(orthogonality(A, x), 0, rtol=0, atol=5e-16)

    def test_rowspace_dense(self):
        A = np.array([[1, 2, 3, 4, 0, 5, 0, 7],
                      [0, 8, 7, 0, 1, 5, 9, 0],
                      [1, 0, 0, 0, 0, 1, 2, 3]])
        test_points = ([1, 2, 3],
                       [1, 10, 3],
                       [1.12, 10, 0])

        for method in available_dense_methods:
            _, _, Y = projections(A, method)
            for z in test_points:
                # Test if x is solution of A x = z
                x = Y.matvec(z)
                assert_array_almost_equal(A.dot(x), z)
                # Test if x is in the return row space of A
                A_ext = np.vstack((A, x))
                assert_equal(np.linalg.matrix_rank(A),
                             np.linalg.matrix_rank(A_ext))


class TestOrthogonality(TestCase):

    def test_dense_matrix(self):
        A = np.array([[1, 2, 3, 4, 0, 5, 0, 7],
                      [0, 8, 7, 0, 1, 5, 9, 0],
                      [1, 0, 0, 0, 0, 1, 2, 3]])
        test_vectors = ([-1.98931144, -1.56363389,
                         -0.84115584, 2.2864762,
                         5.599141, 0.09286976,
                         1.37040802, -0.28145812],
                        [697.92794044, -4091.65114008,
                         -3327.42316335, 836.86906951,
                         99434.98929065, -1285.37653682,
                         -4109.21503806, 2935.29289083])
        test_expected_orth = (0, 0)

        for i in range(len(test_vectors)):
            x = test_vectors[i]
            orth = test_expected_orth[i]
            assert_array_almost_equal(orthogonality(A, x), orth)

    def test_sparse_matrix(self):
        A = np.array([[1, 2, 3, 4, 0, 5, 0, 7],
                      [0, 8, 7, 0, 1, 5, 9, 0],
                      [1, 0, 0, 0, 0, 1, 2, 3]])
        A = csc_matrix(A)
        test_vectors = ([-1.98931144, -1.56363389,
                         -0.84115584, 2.2864762,
                         5.599141, 0.09286976,
                         1.37040802, -0.28145812],
                        [697.92794044, -4091.65114008,
                         -3327.42316335, 836.86906951,
                         99434.98929065, -1285.37653682,
                         -4109.21503806, 2935.29289083])
        test_expected_orth = (0, 0)

        for i in range(len(test_vectors)):
            x = test_vectors[i]
            orth = test_expected_orth[i]
            assert_array_almost_equal(orthogonality(A, x), orth)