#!/usr/bin/env python # # Created by: Pearu Peterson, March 2002 # """ Test functions for linalg.basic module """ from __future__ import division, print_function, absolute_import """ Bugs: 1) solve.check_random_sym_complex fails if a is complex and transpose(a) = conjugate(a) (a is Hermitian). """ __usage__ = """ Build linalg: python setup_linalg.py build Run tests if scipy is installed: python -c 'import scipy;scipy.linalg.test()' Run tests if linalg is not installed: python tests/test_basic.py """ import numpy as np from numpy import (arange, array, dot, zeros, identity, conjugate, transpose, float32) import numpy.linalg as linalg from numpy.random import random from numpy.testing import (TestCase, run_module_suite, assert_raises, assert_equal, assert_almost_equal, assert_array_almost_equal, assert_, assert_allclose, assert_array_equal, dec) from scipy.linalg import (solve, inv, det, lstsq, pinv, pinv2, pinvh, norm, solve_banded, solveh_banded, solve_triangular, solve_circulant, circulant, LinAlgError) from scipy.linalg.basic import LstsqLapackError from scipy.linalg._testutils import assert_no_overwrite from scipy._lib._version import NumpyVersion REAL_DTYPES = [np.float32, np.float64, np.longdouble] COMPLEX_DTYPES = [np.complex64, np.complex128, np.clongdouble] DTYPES = REAL_DTYPES + COMPLEX_DTYPES def _eps_cast(dtyp): """Get the epsilon for dtype, possibly downcast to BLAS types.""" dt = dtyp if dt == np.longdouble: dt = np.float64 elif dt == np.clongdouble: dt = np.complex128 return np.finfo(dt).eps class TestSolveBanded(TestCase): def test_real(self): a = array([[1.0, 20, 0, 0], [-30, 4, 6, 0], [2, 1, 20, 2], [0, -1, 7, 14]]) ab = array([[0.0, 20, 6, 2], [1, 4, 20, 14], [-30, 1, 7, 0], [2, -1, 0, 0]]) l,u = 2,1 b4 = array([10.0, 0.0, 2.0, 14.0]) b4by1 = b4.reshape(-1,1) b4by2 = array([[2, 1], [-30, 4], [2, 3], [1, 3]]) b4by4 = array([[1, 0, 0, 0], [0, 0, 0, 1], [0, 1, 0, 0], [0, 1, 0, 0]]) for b in [b4, b4by1, b4by2, b4by4]: x = solve_banded((l, u), ab, b) assert_array_almost_equal(dot(a, x), b) def test_complex(self): a = array([[1.0, 20, 0, 0], [-30, 4, 6, 0], [2j, 1, 20, 2j], [0, -1, 7, 14]]) ab = array([[0.0, 20, 6, 2j], [1, 4, 20, 14], [-30, 1, 7, 0], [2j, -1, 0, 0]]) l,u = 2,1 b4 = array([10.0, 0.0, 2.0, 14.0j]) b4by1 = b4.reshape(-1,1) b4by2 = array([[2, 1], [-30, 4], [2, 3], [1, 3]]) b4by4 = array([[1, 0, 0, 0], [0, 0, 0,1j], [0, 1, 0, 0], [0, 1, 0, 0]]) for b in [b4, b4by1, b4by2, b4by4]: x = solve_banded((l, u), ab, b) assert_array_almost_equal(dot(a, x), b) def test_tridiag_real(self): ab = array([[0.0, 20, 6, 2], [1, 4, 20, 14], [-30, 1, 7, 0]]) a = np.diag(ab[0,1:], 1) + np.diag(ab[1,:], 0) + np.diag(ab[2,:-1], -1) b4 = array([10.0, 0.0, 2.0, 14.0]) b4by1 = b4.reshape(-1,1) b4by2 = array([[2, 1], [-30, 4], [2, 3], [1, 3]]) b4by4 = array([[1, 0, 0, 0], [0, 0, 0, 1], [0, 1, 0, 0], [0, 1, 0, 0]]) for b in [b4, b4by1, b4by2, b4by4]: x = solve_banded((1, 1), ab, b) assert_array_almost_equal(dot(a, x), b) def test_tridiag_complex(self): ab = array([[0.0, 20, 6, 2j], [1, 4, 20, 14], [-30, 1, 7, 0]]) a = np.diag(ab[0,1:], 1) + np.diag(ab[1,:], 0) + np.diag(ab[2,:-1], -1) b4 = array([10.0, 0.0, 2.0, 14.0j]) b4by1 = b4.reshape(-1,1) b4by2 = array([[2, 1], [-30, 4], [2, 3], [1, 3]]) b4by4 = array([[1, 0, 0, 0], [0, 0, 0, 1], [0, 1, 0, 0], [0, 1, 0, 0]]) for b in [b4, b4by1, b4by2, b4by4]: x = solve_banded((1, 1), ab, b) assert_array_almost_equal(dot(a, x), b) def test_check_finite(self): a = array([[1.0, 20, 0, 0], [-30, 4, 6, 0], [2, 1, 20, 2], [0, -1, 7, 14]]) ab = array([[0.0, 20, 6, 2], [1, 4, 20, 14], [-30, 1, 7, 0], [2, -1, 0, 0]]) l,u = 2,1 b4 = array([10.0, 0.0, 2.0, 14.0]) x = solve_banded((l, u), ab, b4, check_finite=False) assert_array_almost_equal(dot(a, x), b4) def test_bad_shape(self): ab = array([[0.0, 20, 6, 2], [1, 4, 20, 14], [-30, 1, 7, 0], [2, -1, 0, 0]]) l,u = 2,1 bad = array([1.0, 2.0, 3.0, 4.0]).reshape(-1,4) assert_raises(ValueError, solve_banded, (l, u), ab, bad) assert_raises(ValueError, solve_banded, (l, u), ab, [1.0, 2.0]) # Values of (l,u) are not compatible with ab. assert_raises(ValueError, solve_banded, (1, 1), ab, [1.0, 2.0]) def test_1x1(self): b = array([[1., 2., 3.]]) x = solve_banded((1, 1), [[0], [2], [0]], b) assert_array_equal(x, [[0.5, 1.0, 1.5]]) assert_equal(x.dtype, np.dtype('f8')) assert_array_equal(b, [[1.0, 2.0, 3.0]]) class TestSolveHBanded(TestCase): def test_01_upper(self): # Solve # [ 4 1 2 0] [1] # [ 1 4 1 2] X = [4] # [ 2 1 4 1] [1] # [ 0 2 1 4] [2] # with the RHS as a 1D array. ab = array([[0.0, 0.0, 2.0, 2.0], [-99, 1.0, 1.0, 1.0], [4.0, 4.0, 4.0, 4.0]]) b = array([1.0, 4.0, 1.0, 2.0]) x = solveh_banded(ab, b) assert_array_almost_equal(x, [0.0, 1.0, 0.0, 0.0]) def test_02_upper(self): # Solve # [ 4 1 2 0] [1 6] # [ 1 4 1 2] X = [4 2] # [ 2 1 4 1] [1 6] # [ 0 2 1 4] [2 1] # ab = array([[0.0, 0.0, 2.0, 2.0], [-99, 1.0, 1.0, 1.0], [4.0, 4.0, 4.0, 4.0]]) b = array([[1.0, 6.0], [4.0, 2.0], [1.0, 6.0], [2.0, 1.0]]) x = solveh_banded(ab, b) expected = array([[0.0, 1.0], [1.0, 0.0], [0.0, 1.0], [0.0, 0.0]]) assert_array_almost_equal(x, expected) def test_03_upper(self): # Solve # [ 4 1 2 0] [1] # [ 1 4 1 2] X = [4] # [ 2 1 4 1] [1] # [ 0 2 1 4] [2] # with the RHS as a 2D array with shape (3,1). ab = array([[0.0, 0.0, 2.0, 2.0], [-99, 1.0, 1.0, 1.0], [4.0, 4.0, 4.0, 4.0]]) b = array([1.0, 4.0, 1.0, 2.0]).reshape(-1,1) x = solveh_banded(ab, b) assert_array_almost_equal(x, array([0.0, 1.0, 0.0, 0.0]).reshape(-1,1)) def test_01_lower(self): # Solve # [ 4 1 2 0] [1] # [ 1 4 1 2] X = [4] # [ 2 1 4 1] [1] # [ 0 2 1 4] [2] # ab = array([[4.0, 4.0, 4.0, 4.0], [1.0, 1.0, 1.0, -99], [2.0, 2.0, 0.0, 0.0]]) b = array([1.0, 4.0, 1.0, 2.0]) x = solveh_banded(ab, b, lower=True) assert_array_almost_equal(x, [0.0, 1.0, 0.0, 0.0]) def test_02_lower(self): # Solve # [ 4 1 2 0] [1 6] # [ 1 4 1 2] X = [4 2] # [ 2 1 4 1] [1 6] # [ 0 2 1 4] [2 1] # ab = array([[4.0, 4.0, 4.0, 4.0], [1.0, 1.0, 1.0, -99], [2.0, 2.0, 0.0, 0.0]]) b = array([[1.0, 6.0], [4.0, 2.0], [1.0, 6.0], [2.0, 1.0]]) x = solveh_banded(ab, b, lower=True) expected = array([[0.0, 1.0], [1.0, 0.0], [0.0, 1.0], [0.0, 0.0]]) assert_array_almost_equal(x, expected) def test_01_float32(self): # Solve # [ 4 1 2 0] [1] # [ 1 4 1 2] X = [4] # [ 2 1 4 1] [1] # [ 0 2 1 4] [2] # ab = array([[0.0, 0.0, 2.0, 2.0], [-99, 1.0, 1.0, 1.0], [4.0, 4.0, 4.0, 4.0]], dtype=float32) b = array([1.0, 4.0, 1.0, 2.0], dtype=float32) x = solveh_banded(ab, b) assert_array_almost_equal(x, [0.0, 1.0, 0.0, 0.0]) def test_02_float32(self): # Solve # [ 4 1 2 0] [1 6] # [ 1 4 1 2] X = [4 2] # [ 2 1 4 1] [1 6] # [ 0 2 1 4] [2 1] # ab = array([[0.0, 0.0, 2.0, 2.0], [-99, 1.0, 1.0, 1.0], [4.0, 4.0, 4.0, 4.0]], dtype=float32) b = array([[1.0, 6.0], [4.0, 2.0], [1.0, 6.0], [2.0, 1.0]], dtype=float32) x = solveh_banded(ab, b) expected = array([[0.0, 1.0], [1.0, 0.0], [0.0, 1.0], [0.0, 0.0]]) assert_array_almost_equal(x, expected) def test_01_complex(self): # Solve # [ 4 -j 2 0] [2-j] # [ j 4 -j 2] X = [4-j] # [ 2 j 4 -j] [4+j] # [ 0 2 j 4] [2+j] # ab = array([[0.0, 0.0, 2.0, 2.0], [-99, -1.0j, -1.0j, -1.0j], [4.0, 4.0, 4.0, 4.0]]) b = array([2-1.0j, 4.0-1j, 4+1j, 2+1j]) x = solveh_banded(ab, b) assert_array_almost_equal(x, [0.0, 1.0, 1.0, 0.0]) def test_02_complex(self): # Solve # [ 4 -j 2 0] [2-j 2+4j] # [ j 4 -j 2] X = [4-j -1-j] # [ 2 j 4 -j] [4+j 4+2j] # [ 0 2 j 4] [2+j j] # ab = array([[0.0, 0.0, 2.0, 2.0], [-99, -1.0j, -1.0j, -1.0j], [4.0, 4.0, 4.0, 4.0]]) b = array([[2-1j, 2+4j], [4.0-1j, -1-1j], [4.0+1j, 4+2j], [2+1j, 1j]]) x = solveh_banded(ab, b) expected = array([[0.0, 1.0j], [1.0, 0.0], [1.0, 1.0], [0.0, 0.0]]) assert_array_almost_equal(x, expected) def test_tridiag_01_upper(self): # Solve # [ 4 1 0] [1] # [ 1 4 1] X = [4] # [ 0 1 4] [1] # with the RHS as a 1D array. ab = array([[-99, 1.0, 1.0], [4.0, 4.0, 4.0]]) b = array([1.0, 4.0, 1.0]) x = solveh_banded(ab, b) assert_array_almost_equal(x, [0.0, 1.0, 0.0]) def test_tridiag_02_upper(self): # Solve # [ 4 1 0] [1 4] # [ 1 4 1] X = [4 2] # [ 0 1 4] [1 4] # ab = array([[-99, 1.0, 1.0], [4.0, 4.0, 4.0]]) b = array([[1.0, 4.0], [4.0, 2.0], [1.0, 4.0]]) x = solveh_banded(ab, b) expected = array([[0.0, 1.0], [1.0, 0.0], [0.0, 1.0]]) assert_array_almost_equal(x, expected) def test_tridiag_03_upper(self): # Solve # [ 4 1 0] [1] # [ 1 4 1] X = [4] # [ 0 1 4] [1] # with the RHS as a 2D array with shape (3,1). ab = array([[-99, 1.0, 1.0], [4.0, 4.0, 4.0]]) b = array([1.0, 4.0, 1.0]).reshape(-1,1) x = solveh_banded(ab, b) assert_array_almost_equal(x, array([0.0, 1.0, 0.0]).reshape(-1,1)) def test_tridiag_01_lower(self): # Solve # [ 4 1 0] [1] # [ 1 4 1] X = [4] # [ 0 1 4] [1] # ab = array([[4.0, 4.0, 4.0], [1.0, 1.0, -99]]) b = array([1.0, 4.0, 1.0]) x = solveh_banded(ab, b, lower=True) assert_array_almost_equal(x, [0.0, 1.0, 0.0]) def test_tridiag_02_lower(self): # Solve # [ 4 1 0] [1 4] # [ 1 4 1] X = [4 2] # [ 0 1 4] [1 4] # ab = array([[4.0, 4.0, 4.0], [1.0, 1.0, -99]]) b = array([[1.0, 4.0], [4.0, 2.0], [1.0, 4.0]]) x = solveh_banded(ab, b, lower=True) expected = array([[0.0, 1.0], [1.0, 0.0], [0.0, 1.0]]) assert_array_almost_equal(x, expected) def test_tridiag_01_float32(self): # Solve # [ 4 1 0] [1] # [ 1 4 1] X = [4] # [ 0 1 4] [1] # ab = array([[-99, 1.0, 1.0], [4.0, 4.0, 4.0]], dtype=float32) b = array([1.0, 4.0, 1.0], dtype=float32) x = solveh_banded(ab, b) assert_array_almost_equal(x, [0.0, 1.0, 0.0]) def test_tridiag_02_float32(self): # Solve # [ 4 1 0] [1 4] # [ 1 4 1] X = [4 2] # [ 0 1 4] [1 4] # ab = array([[-99, 1.0, 1.0], [4.0, 4.0, 4.0]], dtype=float32) b = array([[1.0, 4.0], [4.0, 2.0], [1.0, 4.0]], dtype=float32) x = solveh_banded(ab, b) expected = array([[0.0, 1.0], [1.0, 0.0], [0.0, 1.0]]) assert_array_almost_equal(x, expected) def test_tridiag_01_complex(self): # Solve # [ 4 -j 0] [ -j] # [ j 4 -j] X = [4-j] # [ 0 j 4] [4+j] # ab = array([[-99, -1.0j, -1.0j], [4.0, 4.0, 4.0]]) b = array([-1.0j, 4.0-1j, 4+1j]) x = solveh_banded(ab, b) assert_array_almost_equal(x, [0.0, 1.0, 1.0]) def test_tridiag_02_complex(self): # Solve # [ 4 -j 0] [ -j 4j] # [ j 4 -j] X = [4-j -1-j] # [ 0 j 4] [4+j 4 ] # ab = array([[-99, -1.0j, -1.0j], [4.0, 4.0, 4.0]]) b = array([[-1j, 4.0j], [4.0-1j, -1.0-1j], [4.0+1j, 4.0]]) x = solveh_banded(ab, b) expected = array([[0.0, 1.0j], [1.0, 0.0], [1.0, 1.0]]) assert_array_almost_equal(x, expected) def test_check_finite(self): # Solve # [ 4 1 0] [1] # [ 1 4 1] X = [4] # [ 0 1 4] [1] # with the RHS as a 1D array. ab = array([[-99, 1.0, 1.0], [4.0, 4.0, 4.0]]) b = array([1.0, 4.0, 1.0]) x = solveh_banded(ab, b, check_finite=False) assert_array_almost_equal(x, [0.0, 1.0, 0.0]) def test_bad_shapes(self): ab = array([[-99, 1.0, 1.0], [4.0, 4.0, 4.0]]) b = array([[1.0, 4.0], [4.0, 2.0]]) assert_raises(ValueError, solveh_banded, ab, b) assert_raises(ValueError, solveh_banded, ab, [1.0, 2.0]) assert_raises(ValueError, solveh_banded, ab, [1.0]) def test_1x1(self): x = solveh_banded([[1]], [[1, 2, 3]]) assert_array_equal(x, [[1.0, 2.0, 3.0]]) assert_equal(x.dtype, np.dtype('f8')) class TestSolve(TestCase): def setUp(self): np.random.seed(1234) def test_20Feb04_bug(self): a = [[1,1],[1.0,0]] # ok x0 = solve(a,[1,0j]) assert_array_almost_equal(dot(a,x0),[1,0]) a = [[1,1],[1.2,0]] # gives failure with clapack.zgesv(..,rowmajor=0) b = [1,0j] x0 = solve(a,b) assert_array_almost_equal(dot(a,x0),[1,0]) def test_simple(self): a = [[1,20],[-30,4]] for b in ([[1,0],[0,1]],[1,0], [[2,1],[-30,4]]): x = solve(a,b) assert_array_almost_equal(dot(a,x),b) def test_simple_sym(self): a = [[2,3],[3,5]] for lower in [0,1]: for b in ([[1,0],[0,1]],[1,0]): x = solve(a,b,sym_pos=1,lower=lower) assert_array_almost_equal(dot(a,x),b) def test_simple_sym_complex(self): a = [[5,2],[2,4]] for b in [[1j,0], [[1j,1j], [0,2]], ]: x = solve(a,b,sym_pos=1) assert_array_almost_equal(dot(a,x),b) def test_simple_complex(self): a = array([[5,2],[2j,4]],'D') for b in [[1j,0], [[1j,1j], [0,2]], [1,0j], array([1,0],'D'), ]: x = solve(a,b) assert_array_almost_equal(dot(a,x),b) def test_nils_20Feb04(self): n = 2 A = random([n,n])+random([n,n])*1j X = zeros((n,n),'D') Ainv = inv(A) R = identity(n)+identity(n)*0j for i in arange(0,n): r = R[:,i] X[:,i] = solve(A,r) assert_array_almost_equal(X,Ainv) def test_random(self): n = 20 a = random([n,n]) for i in range(n): a[i,i] = 20*(.1+a[i,i]) for i in range(4): b = random([n,3]) x = solve(a,b) assert_array_almost_equal(dot(a,x),b) def test_random_complex(self): n = 20 a = random([n,n]) + 1j * random([n,n]) for i in range(n): a[i,i] = 20*(.1+a[i,i]) for i in range(2): b = random([n,3]) x = solve(a,b) assert_array_almost_equal(dot(a,x),b) def test_random_sym(self): n = 20 a = random([n,n]) for i in range(n): a[i,i] = abs(20*(.1+a[i,i])) for j in range(i): a[i,j] = a[j,i] for i in range(4): b = random([n]) x = solve(a,b,sym_pos=1) assert_array_almost_equal(dot(a,x),b) def test_random_sym_complex(self): n = 20 a = random([n,n]) # a = a + 1j*random([n,n]) # XXX: with this the accuracy will be very low for i in range(n): a[i,i] = abs(20*(.1+a[i,i])) for j in range(i): a[i,j] = conjugate(a[j,i]) b = random([n])+2j*random([n]) for i in range(2): x = solve(a,b,sym_pos=1) assert_array_almost_equal(dot(a,x),b) def test_check_finite(self): a = [[1,20],[-30,4]] for b in ([[1,0],[0,1]],[1,0], [[2,1],[-30,4]]): x = solve(a,b, check_finite=False) assert_array_almost_equal(dot(a,x),b) class TestSolveTriangular(TestCase): def test_simple(self): """ solve_triangular on a simple 2x2 matrix. """ A = array([[1,0], [1,2]]) b = [1, 1] sol = solve_triangular(A, b, lower=True) assert_array_almost_equal(sol, [1, 0]) # check that it works also for non-contiguous matrices sol = solve_triangular(A.T, b, lower=False) assert_array_almost_equal(sol, [.5, .5]) # and that it gives the same result as trans=1 sol = solve_triangular(A, b, lower=True, trans=1) assert_array_almost_equal(sol, [.5, .5]) b = identity(2) sol = solve_triangular(A, b, lower=True, trans=1) assert_array_almost_equal(sol, [[1., -.5], [0, 0.5]]) def test_simple_complex(self): """ solve_triangular on a simple 2x2 complex matrix """ A = array([[1+1j, 0], [1j, 2]]) b = identity(2) sol = solve_triangular(A, b, lower=True, trans=1) assert_array_almost_equal(sol, [[.5-.5j, -.25-.25j], [0, 0.5]]) def test_check_finite(self): """ solve_triangular on a simple 2x2 matrix. """ A = array([[1,0], [1,2]]) b = [1, 1] sol = solve_triangular(A, b, lower=True, check_finite=False) assert_array_almost_equal(sol, [1, 0]) class TestInv(TestCase): def setUp(self): np.random.seed(1234) def test_simple(self): a = [[1,2],[3,4]] a_inv = inv(a) assert_array_almost_equal(dot(a,a_inv), [[1,0],[0,1]]) a = [[1,2,3],[4,5,6],[7,8,10]] a_inv = inv(a) assert_array_almost_equal(dot(a,a_inv), [[1,0,0],[0,1,0],[0,0,1]]) def test_random(self): n = 20 for i in range(4): a = random([n,n]) for i in range(n): a[i,i] = 20*(.1+a[i,i]) a_inv = inv(a) assert_array_almost_equal(dot(a,a_inv), identity(n)) def test_simple_complex(self): a = [[1,2],[3,4j]] a_inv = inv(a) assert_array_almost_equal(dot(a,a_inv), [[1,0],[0,1]]) def test_random_complex(self): n = 20 for i in range(4): a = random([n,n])+2j*random([n,n]) for i in range(n): a[i,i] = 20*(.1+a[i,i]) a_inv = inv(a) assert_array_almost_equal(dot(a,a_inv), identity(n)) def test_check_finite(self): a = [[1,2],[3,4]] a_inv = inv(a, check_finite=False) assert_array_almost_equal(dot(a,a_inv), [[1,0],[0,1]]) class TestDet(TestCase): def setUp(self): np.random.seed(1234) def test_simple(self): a = [[1,2],[3,4]] a_det = det(a) assert_almost_equal(a_det,-2.0) def test_simple_complex(self): a = [[1,2],[3,4j]] a_det = det(a) assert_almost_equal(a_det,-6+4j) def test_random(self): basic_det = linalg.det n = 20 for i in range(4): a = random([n,n]) d1 = det(a) d2 = basic_det(a) assert_almost_equal(d1,d2) def test_random_complex(self): basic_det = linalg.det n = 20 for i in range(4): a = random([n,n]) + 2j*random([n,n]) d1 = det(a) d2 = basic_det(a) assert_allclose(d1, d2, rtol=1e-13) def test_check_finite(self): a = [[1,2],[3,4]] a_det = det(a, check_finite=False) assert_almost_equal(a_det,-2.0) def direct_lstsq(a,b,cmplx=0): at = transpose(a) if cmplx: at = conjugate(at) a1 = dot(at, a) b1 = dot(at, b) return solve(a1, b1) class TestLstsq(TestCase): lapack_drivers = ('gelsd', 'gelss', 'gelsy', None) def setUp(self): np.random.seed(1234) def test_simple_exact(self): for dtype in REAL_DTYPES: a = np.array([[1,20], [-30,4]], dtype=dtype) for lapack_driver in TestLstsq.lapack_drivers: for overwrite in (True, False): for bt in (((1,0),(0,1)), (1,0),((2,1),(-30,4))): # Store values in case they are overwritten # later a1 = a.copy() b = np.array(bt, dtype=dtype) b1 = b.copy() try: out = lstsq(a1, b1, lapack_driver=lapack_driver, overwrite_a=overwrite, overwrite_b=overwrite) except LstsqLapackError: if lapack_driver is None: mesg = ('LstsqLapackError raised with ' 'lapack_driver being None.') raise AssertionError(mesg) else: # can't proceed, skip to the next iteration continue x = out[0] r = out[2] assert_(r == 2, 'expected efficient rank 2, got %s' % r) assert_allclose(dot(a, x), b, atol=25 * _eps_cast(a1.dtype), rtol=25 * _eps_cast(a1.dtype), err_msg="driver: %s" % lapack_driver) def test_simple_overdet(self): for dtype in REAL_DTYPES: a = np.array([[1,2], [4,5], [3,4]], dtype=dtype) b = np.array([1,2,3], dtype=dtype) for lapack_driver in TestLstsq.lapack_drivers: for overwrite in (True, False): # Store values in case they are overwritten later a1 = a.copy() b1 = b.copy() try: out = lstsq(a1, b1, lapack_driver=lapack_driver, overwrite_a=overwrite, overwrite_b=overwrite) except LstsqLapackError: if lapack_driver is None: mesg = ('LstsqLapackError raised with ' 'lapack_driver being None.') raise AssertionError(mesg) else: # can't proceed, skip to the next iteration continue x = out[0] if lapack_driver == 'gelsy': residuals = np.sum((b - a.dot(x))**2) else: residuals = out[1] r = out[2] assert_(r == 2, 'expected efficient rank 2, got %s' % r) assert_allclose(abs((dot(a,x) - b)**2).sum(axis=0), residuals, rtol=25 * _eps_cast(a1.dtype), atol=25 * _eps_cast(a1.dtype), err_msg="driver: %s" % lapack_driver) assert_allclose(x, (-0.428571428571429, 0.85714285714285), rtol=25 * _eps_cast(a1.dtype), atol=25 * _eps_cast(a1.dtype), err_msg="driver: %s" % lapack_driver) def test_simple_overdet_complex(self): for dtype in COMPLEX_DTYPES: a = np.array([[1+2j,2], [4,5], [3,4]], dtype=dtype) b = np.array([1, 2+4j, 3], dtype=dtype) for lapack_driver in TestLstsq.lapack_drivers: for overwrite in (True, False): # Store values in case they are overwritten later a1 = a.copy() b1 = b.copy() try: out = lstsq(a1, b1, lapack_driver=lapack_driver, overwrite_a=overwrite, overwrite_b=overwrite) except LstsqLapackError: if lapack_driver is None: mesg = ('LstsqLapackError raised with ' 'lapack_driver being None.') raise AssertionError(mesg) else: # can't proceed, skip to the next iteration continue x = out[0] if lapack_driver == 'gelsy': res = b - a.dot(x) residuals = np.sum(res * res.conj()) else: residuals = out[1] r = out[2] assert_(r == 2, 'expected efficient rank 2, got %s' % r) assert_allclose(abs((dot(a,x) - b)**2).sum(axis=0), residuals, rtol=25 * _eps_cast(a1.dtype), atol=25 * _eps_cast(a1.dtype), err_msg="driver: %s" % lapack_driver) assert_allclose(x, (-0.4831460674157303 + 0.258426966292135j, 0.921348314606741 + 0.292134831460674j), rtol=25 * _eps_cast(a1.dtype), atol=25 * _eps_cast(a1.dtype), err_msg="driver: %s" % lapack_driver) def test_simple_underdet(self): for dtype in REAL_DTYPES: a = np.array([[1,2,3], [4,5,6]], dtype=dtype) b = np.array([1,2], dtype=dtype) for lapack_driver in TestLstsq.lapack_drivers: for overwrite in (True, False): # Store values in case they are overwritten later a1 = a.copy() b1 = b.copy() try: out = lstsq(a1, b1, lapack_driver=lapack_driver, overwrite_a=overwrite, overwrite_b=overwrite) except LstsqLapackError: if lapack_driver is None: mesg = ('LstsqLapackError raised with ' 'lapack_driver being None.') raise AssertionError(mesg) else: # can't proceed, skip to the next iteration continue x = out[0] r = out[2] assert_(r == 2, 'expected efficient rank 2, got %s' % r) assert_allclose(x, (-0.055555555555555, 0.111111111111111, 0.277777777777777), rtol=25 * _eps_cast(a1.dtype), atol=25 * _eps_cast(a1.dtype), err_msg="driver: %s" % lapack_driver) def test_random_exact(self): for dtype in REAL_DTYPES: for n in (20, 200): for lapack_driver in TestLstsq.lapack_drivers: for overwrite in (True, False): a = np.asarray(random([n, n]), dtype=dtype) for i in range(n): a[i, i] = 20 * (0.1 + a[i, i]) for i in range(4): b = np.asarray(random([n, 3]), dtype=dtype) # Store values in case they are overwritten later a1 = a.copy() b1 = b.copy() try: out = lstsq(a1, b1, lapack_driver=lapack_driver, overwrite_a=overwrite, overwrite_b=overwrite) except LstsqLapackError: if lapack_driver is None: mesg = ('LstsqLapackError raised with ' 'lapack_driver being None.') raise AssertionError(mesg) else: # can't proceed, skip to the next iteration continue x = out[0] r = out[2] assert_(r == n, 'expected efficient rank %s, got ' '%s' % (n, r)) if dtype is np.float32: assert_allclose(dot(a, x), b, rtol=500 * _eps_cast(a1.dtype), atol=500 * _eps_cast(a1.dtype), err_msg="driver: %s" % lapack_driver) else: assert_allclose(dot(a, x), b, rtol=1000 * _eps_cast(a1.dtype), atol=1000 * _eps_cast(a1.dtype), err_msg="driver: %s" % lapack_driver) def test_random_complex_exact(self): for dtype in COMPLEX_DTYPES: for n in (20, 200): for lapack_driver in TestLstsq.lapack_drivers: for overwrite in (True, False): a = np.asarray(random([n, n]) + 1j*random([n, n]), dtype=dtype) for i in range(n): a[i, i] = 20 * (0.1 + a[i, i]) for i in range(2): b = np.asarray(random([n, 3]), dtype=dtype) # Store values in case they are overwritten later a1 = a.copy() b1 = b.copy() out = lstsq(a1, b1, lapack_driver=lapack_driver, overwrite_a=overwrite, overwrite_b=overwrite) x = out[0] r = out[2] assert_(r == n, 'expected efficient rank %s, got ' '%s' % (n, r)) if dtype is np.complex64: assert_allclose(dot(a, x), b, rtol=400 * _eps_cast(a1.dtype), atol=400 * _eps_cast(a1.dtype), err_msg="driver: %s" % lapack_driver) else: assert_allclose(dot(a, x), b, rtol=1000 * _eps_cast(a1.dtype), atol=1000 * _eps_cast(a1.dtype), err_msg="driver: %s" % lapack_driver) def test_random_overdet(self): for dtype in REAL_DTYPES: for (n,m) in ((20,15), (200,2)): for lapack_driver in TestLstsq.lapack_drivers: for overwrite in (True, False): a = np.asarray(random([n,m]), dtype=dtype) for i in range(m): a[i,i] = 20 * (0.1 + a[i,i]) for i in range(4): b = np.asarray(random([n,3]), dtype=dtype) # Store values in case they are overwritten later a1 = a.copy() b1 = b.copy() try: out = lstsq(a1, b1, lapack_driver=lapack_driver, overwrite_a=overwrite, overwrite_b=overwrite) except LstsqLapackError: if lapack_driver is None: mesg = ('LstsqLapackError raised with ' 'lapack_driver being None.') raise AssertionError(mesg) else: # can't proceed, skip to the next iteration continue x = out[0] r = out[2] assert_(r == m, 'expected efficient rank %s, got ' '%s' % (m, r)) assert_allclose(x, direct_lstsq(a, b, cmplx=0), rtol=25 * _eps_cast(a1.dtype), atol=25 * _eps_cast(a1.dtype), err_msg="driver: %s" % lapack_driver) def test_random_complex_overdet(self): for dtype in COMPLEX_DTYPES: for (n, m) in ((20, 15), (200, 2)): for lapack_driver in TestLstsq.lapack_drivers: for overwrite in (True, False): a = np.asarray(random([n, m]) + 1j*random([n, m]), dtype=dtype) for i in range(m): a[i, i] = 20 * (0.1 + a[i, i]) for i in range(2): b = np.asarray(random([n, 3]), dtype=dtype) # Store values in case they are overwritten # later a1 = a.copy() b1 = b.copy() out = lstsq(a1, b1, lapack_driver=lapack_driver, overwrite_a=overwrite, overwrite_b=overwrite) x = out[0] r = out[2] assert_(r == m, 'expected efficient rank %s, got ' '%s' % (m, r)) assert_allclose(x, direct_lstsq(a, b, cmplx=1), rtol=25 * _eps_cast(a1.dtype), atol=25 * _eps_cast(a1.dtype), err_msg="driver: %s" % lapack_driver) def test_check_finite(self): for dtype in REAL_DTYPES: a = np.array(((1, 20), (-30, 4)), dtype=dtype) for bt in (((1, 0), (0, 1)), (1, 0), ((2, 1), (-30, 4))): for lapack_driver in TestLstsq.lapack_drivers: for overwrite in (True, False): for check_finite in (True, False): b = np.array(bt, dtype=dtype) # Store values in case they are overwritten # later a1 = a.copy() b1 = b.copy() try: out = lstsq(a1, b1, lapack_driver=lapack_driver, check_finite=check_finite, overwrite_a=overwrite, overwrite_b=overwrite) except LstsqLapackError: if lapack_driver is None: mesg = ('LstsqLapackError raised with ' 'lapack_driver being None.') raise AssertionError(mesg) else: # can't proceed, skip to the next iteration continue x = out[0] r = out[2] assert_(r == 2, 'expected efficient rank 2, ' 'got %s' % r) assert_allclose(dot(a, x), b, rtol=25 * _eps_cast(a.dtype), atol=25 * _eps_cast(a.dtype), err_msg="driver: %s" % lapack_driver) class TestPinv(TestCase): def test_simple_real(self): a = array([[1, 2, 3], [4, 5, 6], [7, 8, 10]], dtype=float) a_pinv = pinv(a) assert_array_almost_equal(dot(a,a_pinv), np.eye(3)) a_pinv = pinv2(a) assert_array_almost_equal(dot(a,a_pinv), np.eye(3)) def test_simple_complex(self): a = (array([[1, 2, 3], [4, 5, 6], [7, 8, 10]], dtype=float) + 1j * array([[10, 8, 7], [6, 5, 4], [3, 2, 1]], dtype=float)) a_pinv = pinv(a) assert_array_almost_equal(dot(a, a_pinv), np.eye(3)) a_pinv = pinv2(a) assert_array_almost_equal(dot(a, a_pinv), np.eye(3)) def test_simple_singular(self): a = array([[1, 2, 3], [4, 5, 6], [7, 8, 9]], dtype=float) a_pinv = pinv(a) a_pinv2 = pinv2(a) assert_array_almost_equal(a_pinv,a_pinv2) def test_simple_cols(self): a = array([[1, 2, 3], [4, 5, 6]], dtype=float) a_pinv = pinv(a) a_pinv2 = pinv2(a) assert_array_almost_equal(a_pinv,a_pinv2) def test_simple_rows(self): a = array([[1, 2], [3, 4], [5, 6]], dtype=float) a_pinv = pinv(a) a_pinv2 = pinv2(a) assert_array_almost_equal(a_pinv,a_pinv2) def test_check_finite(self): a = array([[1,2,3],[4,5,6.],[7,8,10]]) a_pinv = pinv(a, check_finite=False) assert_array_almost_equal(dot(a,a_pinv),[[1,0,0],[0,1,0],[0,0,1]]) a_pinv = pinv2(a, check_finite=False) assert_array_almost_equal(dot(a,a_pinv),[[1,0,0],[0,1,0],[0,0,1]]) class TestPinvSymmetric(TestCase): def test_simple_real(self): a = array([[1, 2, 3], [4, 5, 6], [7, 8, 10]], dtype=float) a = np.dot(a, a.T) a_pinv = pinvh(a) assert_array_almost_equal(np.dot(a, a_pinv), np.eye(3)) def test_nonpositive(self): a = array([[1, 2, 3], [4, 5, 6], [7, 8, 9]], dtype=float) a = np.dot(a, a.T) u, s, vt = np.linalg.svd(a) s[0] *= -1 a = np.dot(u * s, vt) # a is now symmetric non-positive and singular a_pinv = pinv2(a) a_pinvh = pinvh(a) assert_array_almost_equal(a_pinv, a_pinvh) def test_simple_complex(self): a = (array([[1, 2, 3], [4, 5, 6], [7, 8, 10]], dtype=float) + 1j * array([[10, 8, 7], [6, 5, 4], [3, 2, 1]], dtype=float)) a = np.dot(a, a.conj().T) a_pinv = pinvh(a) assert_array_almost_equal(np.dot(a, a_pinv), np.eye(3)) class TestVectorNorms(object): def test_types(self): for dtype in np.typecodes['AllFloat']: x = np.array([1,2,3], dtype=dtype) tol = max(1e-15, np.finfo(dtype).eps.real * 20) assert_allclose(norm(x), np.sqrt(14), rtol=tol) assert_allclose(norm(x, 2), np.sqrt(14), rtol=tol) for dtype in np.typecodes['Complex']: x = np.array([1j,2j,3j], dtype=dtype) tol = max(1e-15, np.finfo(dtype).eps.real * 20) assert_allclose(norm(x), np.sqrt(14), rtol=tol) assert_allclose(norm(x, 2), np.sqrt(14), rtol=tol) def test_overflow(self): # unlike numpy's norm, this one is # safer on overflow a = array([1e20], dtype=float32) assert_almost_equal(norm(a), a) def test_stable(self): # more stable than numpy's norm a = array([1e4] + [1]*10000, dtype=float32) try: # snrm in double precision; we obtain the same as for float64 # -- large atol needed due to varying blas implementations assert_allclose(norm(a) - 1e4, 0.5, atol=1e-2) except AssertionError: # snrm implemented in single precision, == np.linalg.norm result msg = ": Result should equal either 0.0 or 0.5 (depending on " \ "implementation of snrm2)." assert_almost_equal(norm(a) - 1e4, 0.0, err_msg=msg) def test_zero_norm(self): assert_equal(norm([1,0,3], 0), 2) assert_equal(norm([1,2,3], 0), 3) @dec.skipif(NumpyVersion(np.__version__) < '1.8.0') def test_axis_kwd(self): a = np.array([[[2, 1], [3, 4]]] * 2, 'd') assert_allclose(norm(a, axis=1), [[3.60555128, 4.12310563]] * 2) assert_allclose(norm(a, 1, axis=1), [[5.] * 2] * 2) @dec.skipif(NumpyVersion(np.__version__) < '1.10.0') def test_keepdims_kwd(self): a = np.array([[[2, 1], [3, 4]]] * 2, 'd') b = norm(a, axis=1, keepdims=True) assert_allclose(b, [[[3.60555128, 4.12310563]]] * 2) assert_(b.shape == (2, 1, 2)) assert_allclose(norm(a, 1, axis=2, keepdims=True), [[[3.],[7.]]] * 2) class TestMatrixNorms(object): def test_matrix_norms(self): # Not all of these are matrix norms in the most technical sense. np.random.seed(1234) for n, m in (1, 1), (1, 3), (3, 1), (4, 4), (4, 5), (5, 4): for t in np.single, np.double, np.csingle, np.cdouble, np.int64: A = 10 * np.random.randn(n, m).astype(t) if np.issubdtype(A.dtype, np.complexfloating): A = (A + 10j * np.random.randn(n, m)).astype(t) t_high = np.cdouble else: t_high = np.double for order in (None, 'fro', 1, -1, 2, -2, np.inf, -np.inf): actual = norm(A, ord=order) desired = np.linalg.norm(A, ord=order) # SciPy may return higher precision matrix norms. # This is a consequence of using LAPACK. if not np.allclose(actual, desired): desired = np.linalg.norm(A.astype(t_high), ord=order) np.assert_allclose(actual, desired) @dec.skipif(NumpyVersion(np.__version__) < '1.8.0') def test_axis_kwd(self): a = np.array([[[2, 1], [3, 4]]] * 2, 'd') b = norm(a, ord=np.inf, axis=(1, 0)) c = norm(np.swapaxes(a, 0, 1), ord=np.inf, axis=(0, 1)) d = norm(a, ord=1, axis=(0, 1)) assert_allclose(b, c) assert_allclose(c, d) assert_allclose(b, d) assert_(b.shape == c.shape == d.shape) b = norm(a, ord=1, axis=(1, 0)) c = norm(np.swapaxes(a, 0, 1), ord=1, axis=(0, 1)) d = norm(a, ord=np.inf, axis=(0, 1)) assert_allclose(b, c) assert_allclose(c, d) assert_allclose(b, d) assert_(b.shape == c.shape == d.shape) @dec.skipif(NumpyVersion(np.__version__) < '1.10.0') def test_keepdims_kwd(self): a = np.arange(120, dtype='d').reshape(2, 3, 4, 5) b = norm(a, ord=np.inf, axis=(1, 0), keepdims=True) c = norm(a, ord=1, axis=(0, 1), keepdims=True) assert_allclose(b, c) assert_(b.shape == c.shape) class TestOverwrite(object): def test_solve(self): assert_no_overwrite(solve, [(3,3), (3,)]) def test_solve_triangular(self): assert_no_overwrite(solve_triangular, [(3,3), (3,)]) def test_solve_banded(self): assert_no_overwrite(lambda ab, b: solve_banded((2,1), ab, b), [(4,6), (6,)]) def test_solveh_banded(self): assert_no_overwrite(solveh_banded, [(2,6), (6,)]) def test_inv(self): assert_no_overwrite(inv, [(3,3)]) def test_det(self): assert_no_overwrite(det, [(3,3)]) def test_lstsq(self): assert_no_overwrite(lstsq, [(3,2), (3,)]) def test_pinv(self): assert_no_overwrite(pinv, [(3,3)]) def test_pinv2(self): assert_no_overwrite(pinv2, [(3,3)]) def test_pinvh(self): assert_no_overwrite(pinvh, [(3,3)]) class TestSolveCirculant(TestCase): def test_basic1(self): c = np.array([1, 2, 3, 5]) b = np.array([1, -1, 1, 0]) x = solve_circulant(c, b) y = solve(circulant(c), b) assert_allclose(x, y) def test_basic2(self): # b is a 2-d matrix. c = np.array([1, 2, -3, -5]) b = np.arange(12).reshape(4, 3) x = solve_circulant(c, b) y = solve(circulant(c), b) assert_allclose(x, y) def test_basic3(self): # b is a 3-d matrix. c = np.array([1, 2, -3, -5]) b = np.arange(24).reshape(4, 3, 2) x = solve_circulant(c, b) y = solve(circulant(c), b) assert_allclose(x, y) def test_complex(self): # Complex b and c c = np.array([1+2j, -3, 4j, 5]) b = np.arange(8).reshape(4, 2) + 0.5j x = solve_circulant(c, b) y = solve(circulant(c), b) assert_allclose(x, y) def test_random_b_and_c(self): # Random b and c np.random.seed(54321) c = np.random.randn(50) b = np.random.randn(50) x = solve_circulant(c, b) y = solve(circulant(c), b) assert_allclose(x, y) def test_singular(self): # c gives a singular circulant matrix. c = np.array([1, 1, 0, 0]) b = np.array([1, 2, 3, 4]) x = solve_circulant(c, b, singular='lstsq') y, res, rnk, s = lstsq(circulant(c), b) assert_allclose(x, y) assert_raises(LinAlgError, solve_circulant, x, y) def test_axis_args(self): # Test use of caxis, baxis and outaxis. # c has shape (2, 1, 4) c = np.array([[[-1, 2.5, 3, 3.5]], [[1, 6, 6, 6.5]]]) # b has shape (3, 4) b = np.array([[0, 0, 1, 1], [1, 1, 0, 0], [1, -1, 0, 0]]) x = solve_circulant(c, b, baxis=1) assert_equal(x.shape, (4, 2, 3)) expected = np.empty_like(x) expected[:, 0, :] = solve(circulant(c[0]), b.T) expected[:, 1, :] = solve(circulant(c[1]), b.T) assert_allclose(x, expected) x = solve_circulant(c, b, baxis=1, outaxis=-1) assert_equal(x.shape, (2, 3, 4)) assert_allclose(np.rollaxis(x, -1), expected) # np.swapaxes(c, 1, 2) has shape (2, 4, 1); b.T has shape (4, 3). x = solve_circulant(np.swapaxes(c, 1, 2), b.T, caxis=1) assert_equal(x.shape, (4, 2, 3)) assert_allclose(x, expected) if __name__ == "__main__": run_module_suite()