""" Test functions for linalg module """ from __future__ import division, absolute_import, print_function import os import sys import numpy as np from numpy.testing import (TestCase, assert_, assert_equal, assert_raises, assert_array_equal, assert_almost_equal, run_module_suite, dec) from numpy import array, single, double, csingle, cdouble, dot, identity from numpy import multiply, atleast_2d, inf, asarray, matrix from numpy import linalg from numpy.linalg import matrix_power, norm, matrix_rank def ifthen(a, b): return not a or b old_assert_almost_equal = assert_almost_equal def imply(a, b): return not a or b def assert_almost_equal(a, b, **kw): if asarray(a).dtype.type in (single, csingle): decimal = 6 else: decimal = 12 old_assert_almost_equal(a, b, decimal=decimal, **kw) def get_real_dtype(dtype): return {single: single, double: double, csingle: single, cdouble: double}[dtype] def get_complex_dtype(dtype): return {single: csingle, double: cdouble, csingle: csingle, cdouble: cdouble}[dtype] class LinalgTestCase(object): def test_single(self): a = array([[1.,2.], [3.,4.]], dtype=single) b = array([2., 1.], dtype=single) self.do(a, b) def test_double(self): a = array([[1.,2.], [3.,4.]], dtype=double) b = array([2., 1.], dtype=double) self.do(a, b) def test_double_2(self): a = array([[1.,2.], [3.,4.]], dtype=double) b = array([[2., 1., 4.], [3., 4., 6.]], dtype=double) self.do(a, b) def test_csingle(self): a = array([[1.+2j,2+3j], [3+4j,4+5j]], dtype=csingle) b = array([2.+1j, 1.+2j], dtype=csingle) self.do(a, b) def test_cdouble(self): a = array([[1.+2j,2+3j], [3+4j,4+5j]], dtype=cdouble) b = array([2.+1j, 1.+2j], dtype=cdouble) self.do(a, b) def test_cdouble_2(self): a = array([[1.+2j,2+3j], [3+4j,4+5j]], dtype=cdouble) b = array([[2.+1j, 1.+2j, 1+3j], [1-2j, 1-3j, 1-6j]], dtype=cdouble) self.do(a, b) def test_empty(self): a = atleast_2d(array([], dtype = double)) b = atleast_2d(array([], dtype = double)) try: self.do(a, b) raise AssertionError("%s should fail with empty matrices", self.__name__[5:]) except linalg.LinAlgError as e: pass def test_nonarray(self): a = [[1,2], [3,4]] b = [2, 1] self.do(a,b) def test_matrix_b_only(self): """Check that matrix type is preserved.""" a = array([[1.,2.], [3.,4.]]) b = matrix([2., 1.]).T self.do(a, b) def test_matrix_a_and_b(self): """Check that matrix type is preserved.""" a = matrix([[1.,2.], [3.,4.]]) b = matrix([2., 1.]).T self.do(a, b) class LinalgNonsquareTestCase(object): def test_single_nsq_1(self): a = array([[1.,2.,3.], [3.,4.,6.]], dtype=single) b = array([2., 1.], dtype=single) self.do(a, b) def test_single_nsq_2(self): a = array([[1.,2.], [3.,4.], [5.,6.]], dtype=single) b = array([2., 1., 3.], dtype=single) self.do(a, b) def test_double_nsq_1(self): a = array([[1.,2.,3.], [3.,4.,6.]], dtype=double) b = array([2., 1.], dtype=double) self.do(a, b) def test_double_nsq_2(self): a = array([[1.,2.], [3.,4.], [5.,6.]], dtype=double) b = array([2., 1., 3.], dtype=double) self.do(a, b) def test_csingle_nsq_1(self): a = array([[1.+1j,2.+2j,3.-3j], [3.-5j,4.+9j,6.+2j]], dtype=csingle) b = array([2.+1j, 1.+2j], dtype=csingle) self.do(a, b) def test_csingle_nsq_2(self): a = array([[1.+1j,2.+2j], [3.-3j,4.-9j], [5.-4j,6.+8j]], dtype=csingle) b = array([2.+1j, 1.+2j, 3.-3j], dtype=csingle) self.do(a, b) def test_cdouble_nsq_1(self): a = array([[1.+1j,2.+2j,3.-3j], [3.-5j,4.+9j,6.+2j]], dtype=cdouble) b = array([2.+1j, 1.+2j], dtype=cdouble) self.do(a, b) def test_cdouble_nsq_2(self): a = array([[1.+1j,2.+2j], [3.-3j,4.-9j], [5.-4j,6.+8j]], dtype=cdouble) b = array([2.+1j, 1.+2j, 3.-3j], dtype=cdouble) self.do(a, b) def test_cdouble_nsq_1_2(self): a = array([[1.+1j,2.+2j,3.-3j], [3.-5j,4.+9j,6.+2j]], dtype=cdouble) b = array([[2.+1j, 1.+2j], [1-1j, 2-2j]], dtype=cdouble) self.do(a, b) def test_cdouble_nsq_2_2(self): a = array([[1.+1j,2.+2j], [3.-3j,4.-9j], [5.-4j,6.+8j]], dtype=cdouble) b = array([[2.+1j, 1.+2j], [1-1j, 2-2j], [1-1j, 2-2j]], dtype=cdouble) self.do(a, b) def _generalized_testcase(new_cls_name, old_cls): def get_method(old_name, new_name): def method(self): base = old_cls() def do(a, b): a = np.array([a, a, a]) b = np.array([b, b, b]) self.do(a, b) base.do = do getattr(base, old_name)() method.__name__ = new_name return method dct = dict() for old_name in dir(old_cls): if old_name.startswith('test_'): new_name = old_name + '_generalized' dct[new_name] = get_method(old_name, new_name) return type(new_cls_name, (object,), dct) LinalgGeneralizedTestCase = _generalized_testcase( 'LinalgGeneralizedTestCase', LinalgTestCase) LinalgGeneralizedNonsquareTestCase = _generalized_testcase( 'LinalgGeneralizedNonsquareTestCase', LinalgNonsquareTestCase) def dot_generalized(a, b): a = asarray(a) if a.ndim == 3: return np.array([dot(ax, bx) for ax, bx in zip(a, b)]) elif a.ndim > 3: raise ValueError("Not implemented...") return dot(a, b) def identity_like_generalized(a): a = asarray(a) if a.ndim == 3: return np.array([identity(a.shape[-2]) for ax in a]) elif a.ndim > 3: raise ValueError("Not implemented...") return identity(a.shape[0]) class TestSolve(LinalgTestCase, LinalgGeneralizedTestCase, TestCase): def do(self, a, b): x = linalg.solve(a, b) assert_almost_equal(b, dot_generalized(a, x)) assert_(imply(isinstance(b, matrix), isinstance(x, matrix))) def test_types(self): def check(dtype): x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype) assert_equal(linalg.solve(x, x).dtype, dtype) for dtype in [single, double, csingle, cdouble]: yield check, dtype class TestInv(LinalgTestCase, LinalgGeneralizedTestCase, TestCase): def do(self, a, b): a_inv = linalg.inv(a) assert_almost_equal(dot_generalized(a, a_inv), identity_like_generalized(a)) assert_(imply(isinstance(a, matrix), isinstance(a_inv, matrix))) def test_types(self): def check(dtype): x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype) assert_equal(linalg.inv(x).dtype, dtype) for dtype in [single, double, csingle, cdouble]: yield check, dtype class TestEigvals(LinalgTestCase, LinalgGeneralizedTestCase, TestCase): def do(self, a, b): ev = linalg.eigvals(a) evalues, evectors = linalg.eig(a) assert_almost_equal(ev, evalues) def test_types(self): def check(dtype): x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype) assert_equal(linalg.eigvals(x).dtype, dtype) x = np.array([[1, 0.5], [-1, 1]], dtype=dtype) assert_equal(linalg.eigvals(x).dtype, get_complex_dtype(dtype)) for dtype in [single, double, csingle, cdouble]: yield check, dtype class TestEig(LinalgTestCase, LinalgGeneralizedTestCase, TestCase): def do(self, a, b): evalues, evectors = linalg.eig(a) if evectors.ndim == 3: assert_almost_equal(dot_generalized(a, evectors), evectors * evalues[:,None,:]) else: assert_almost_equal(dot(a, evectors), multiply(evectors, evalues)) assert_(imply(isinstance(a, matrix), isinstance(evectors, matrix))) def test_types(self): def check(dtype): x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype) w, v = np.linalg.eig(x) assert_equal(w.dtype, dtype) assert_equal(v.dtype, dtype) x = np.array([[1, 0.5], [-1, 1]], dtype=dtype) w, v = np.linalg.eig(x) assert_equal(w.dtype, get_complex_dtype(dtype)) assert_equal(v.dtype, get_complex_dtype(dtype)) for dtype in [single, double, csingle, cdouble]: yield dtype class TestSVD(LinalgTestCase, LinalgGeneralizedTestCase, TestCase): def do(self, a, b): u, s, vt = linalg.svd(a, 0) if u.ndim == 3: assert_almost_equal(a, dot_generalized(u * s[:,None,:], vt)) else: assert_almost_equal(a, dot(multiply(u, s), vt)) assert_(imply(isinstance(a, matrix), isinstance(u, matrix))) assert_(imply(isinstance(a, matrix), isinstance(vt, matrix))) def test_types(self): def check(dtype): x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype) u, s, vh = linalg.svd(x) assert_equal(u.dtype, dtype) assert_equal(s.dtype, get_real_dtype(dtype)) assert_equal(vh.dtype, dtype) s = linalg.svd(x, compute_uv=False) assert_equal(s.dtype, get_real_dtype(dtype)) for dtype in [single, double, csingle, cdouble]: yield check, dtype class TestCondSVD(LinalgTestCase, LinalgGeneralizedTestCase, TestCase): def do(self, a, b): c = asarray(a) # a might be a matrix s = linalg.svd(c, compute_uv=False) old_assert_almost_equal(s[0]/s[-1], linalg.cond(a), decimal=5) class TestCond2(LinalgTestCase, TestCase): def do(self, a, b): c = asarray(a) # a might be a matrix s = linalg.svd(c, compute_uv=False) old_assert_almost_equal(s[0]/s[-1], linalg.cond(a,2), decimal=5) class TestCondInf(TestCase): def test(self): A = array([[1.,0,0],[0,-2.,0],[0,0,3.]]) assert_almost_equal(linalg.cond(A,inf),3.) class TestPinv(LinalgTestCase, TestCase): def do(self, a, b): a_ginv = linalg.pinv(a) assert_almost_equal(dot(a, a_ginv), identity(asarray(a).shape[0])) assert_(imply(isinstance(a, matrix), isinstance(a_ginv, matrix))) class TestDet(LinalgTestCase, LinalgGeneralizedTestCase, TestCase): def do(self, a, b): d = linalg.det(a) (s, ld) = linalg.slogdet(a) if asarray(a).dtype.type in (single, double): ad = asarray(a).astype(double) else: ad = asarray(a).astype(cdouble) ev = linalg.eigvals(ad) assert_almost_equal(d, multiply.reduce(ev, axis=-1)) assert_almost_equal(s * np.exp(ld), multiply.reduce(ev, axis=-1)) s = np.atleast_1d(s) ld = np.atleast_1d(ld) m = (s != 0) assert_almost_equal(np.abs(s[m]), 1) assert_equal(ld[~m], -inf) def test_zero(self): assert_equal(linalg.det([[0.0]]), 0.0) assert_equal(type(linalg.det([[0.0]])), double) assert_equal(linalg.det([[0.0j]]), 0.0) assert_equal(type(linalg.det([[0.0j]])), cdouble) assert_equal(linalg.slogdet([[0.0]]), (0.0, -inf)) assert_equal(type(linalg.slogdet([[0.0]])[0]), double) assert_equal(type(linalg.slogdet([[0.0]])[1]), double) assert_equal(linalg.slogdet([[0.0j]]), (0.0j, -inf)) assert_equal(type(linalg.slogdet([[0.0j]])[0]), cdouble) assert_equal(type(linalg.slogdet([[0.0j]])[1]), double) def test_types(self): def check(dtype): x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype) assert_equal(np.linalg.det(x), get_real_dtype(dtype)) for dtype in [single, double, csingle, cdouble]: yield check, dtype class TestLstsq(LinalgTestCase, LinalgNonsquareTestCase, TestCase): def do(self, a, b): arr = np.asarray(a) m, n = arr.shape u, s, vt = linalg.svd(a, 0) x, residuals, rank, sv = linalg.lstsq(a, b) if m <= n: assert_almost_equal(b, dot(a, x)) assert_equal(rank, m) else: assert_equal(rank, n) assert_almost_equal(sv, sv.__array_wrap__(s)) if rank == n and m > n: expect_resids = (np.asarray(abs(np.dot(a, x) - b))**2).sum(axis=0) expect_resids = np.asarray(expect_resids) if len(np.asarray(b).shape) == 1: expect_resids.shape = (1,) assert_equal(residuals.shape, expect_resids.shape) else: expect_resids = type(x)([]) assert_almost_equal(residuals, expect_resids) assert_(np.issubdtype(residuals.dtype, np.floating)) assert_(imply(isinstance(b, matrix), isinstance(x, matrix))) assert_(imply(isinstance(b, matrix), isinstance(residuals, matrix))) class TestMatrixPower(object): R90 = array([[0,1],[-1,0]]) Arb22 = array([[4,-7],[-2,10]]) noninv = array([[1,0],[0,0]]) arbfloat = array([[0.1,3.2],[1.2,0.7]]) large = identity(10) t = large[1,:].copy() large[1,:] = large[0,:] large[0,:] = t def test_large_power(self): assert_equal(matrix_power(self.R90,2**100+2**10+2**5+1),self.R90) def test_large_power_trailing_zero(self): assert_equal(matrix_power(self.R90,2**100+2**10+2**5),identity(2)) def testip_zero(self): def tz(M): mz = matrix_power(M,0) assert_equal(mz, identity(M.shape[0])) assert_equal(mz.dtype, M.dtype) for M in [self.Arb22, self.arbfloat, self.large]: yield tz, M def testip_one(self): def tz(M): mz = matrix_power(M,1) assert_equal(mz, M) assert_equal(mz.dtype, M.dtype) for M in [self.Arb22, self.arbfloat, self.large]: yield tz, M def testip_two(self): def tz(M): mz = matrix_power(M,2) assert_equal(mz, dot(M,M)) assert_equal(mz.dtype, M.dtype) for M in [self.Arb22, self.arbfloat, self.large]: yield tz, M def testip_invert(self): def tz(M): mz = matrix_power(M,-1) assert_almost_equal(identity(M.shape[0]), dot(mz,M)) for M in [self.R90, self.Arb22, self.arbfloat, self.large]: yield tz, M def test_invert_noninvertible(self): import numpy.linalg assert_raises(numpy.linalg.linalg.LinAlgError, lambda: matrix_power(self.noninv,-1)) class TestBoolPower(TestCase): def test_square(self): A = array([[True,False],[True,True]]) assert_equal(matrix_power(A,2),A) class HermitianTestCase(object): def test_single(self): a = array([[1.,2.], [2.,1.]], dtype=single) self.do(a, None) def test_double(self): a = array([[1.,2.], [2.,1.]], dtype=double) self.do(a, None) def test_csingle(self): a = array([[1.,2+3j], [2-3j,1]], dtype=csingle) self.do(a, None) def test_cdouble(self): a = array([[1.,2+3j], [2-3j,1]], dtype=cdouble) self.do(a, None) def test_empty(self): a = atleast_2d(array([], dtype = double)) assert_raises(linalg.LinAlgError, self.do, a, None) def test_nonarray(self): a = [[1,2], [2,1]] self.do(a, None) def test_matrix_b_only(self): """Check that matrix type is preserved.""" a = array([[1.,2.], [2.,1.]]) self.do(a, None) def test_matrix_a_and_b(self): """Check that matrix type is preserved.""" a = matrix([[1.,2.], [2.,1.]]) self.do(a, None) HermitianGeneralizedTestCase = _generalized_testcase( 'HermitianGeneralizedTestCase', HermitianTestCase) class TestEigvalsh(HermitianTestCase, HermitianGeneralizedTestCase, TestCase): def do(self, a, b): # note that eigenvalue arrays must be sorted since # their order isn't guaranteed. ev = linalg.eigvalsh(a) evalues, evectors = linalg.eig(a) ev.sort(axis=-1) evalues.sort(axis=-1) assert_almost_equal(ev, evalues) def test_types(self): def check(dtype): x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype) assert_equal(np.linalg.eigvalsh(x), get_real_dtype(dtype)) for dtype in [single, double, csingle, cdouble]: yield check, dtype class TestEigh(HermitianTestCase, HermitianGeneralizedTestCase, TestCase): def do(self, a, b): # note that eigenvalue arrays must be sorted since # their order isn't guaranteed. ev, evc = linalg.eigh(a) evalues, evectors = linalg.eig(a) ev.sort(axis=-1) evalues.sort(axis=-1) assert_almost_equal(ev, evalues) def test_types(self): def check(dtype): x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype) w, v = np.linalg.eig(x) assert_equal(w, get_real_dtype(dtype)) assert_equal(v, dtype) for dtype in [single, double, csingle, cdouble]: yield check, dtype class _TestNorm(TestCase): dt = None dec = None def test_empty(self): assert_equal(norm([]), 0.0) assert_equal(norm(array([], dtype=self.dt)), 0.0) assert_equal(norm(atleast_2d(array([], dtype=self.dt))), 0.0) def test_vector(self): a = [1, 2, 3, 4] b = [-1, -2, -3, -4] c = [-1, 2, -3, 4] def _test(v): np.testing.assert_almost_equal(norm(v), 30**0.5, decimal=self.dec) np.testing.assert_almost_equal(norm(v, inf), 4.0, decimal=self.dec) np.testing.assert_almost_equal(norm(v, -inf), 1.0, decimal=self.dec) np.testing.assert_almost_equal(norm(v, 1), 10.0, decimal=self.dec) np.testing.assert_almost_equal(norm(v, -1), 12.0/25, decimal=self.dec) np.testing.assert_almost_equal(norm(v, 2), 30**0.5, decimal=self.dec) np.testing.assert_almost_equal(norm(v, -2), ((205./144)**-0.5), decimal=self.dec) np.testing.assert_almost_equal(norm(v, 0), 4, decimal=self.dec) for v in (a, b, c,): _test(v) for v in (array(a, dtype=self.dt), array(b, dtype=self.dt), array(c, dtype=self.dt)): _test(v) def test_matrix(self): A = matrix([[1, 3], [5, 7]], dtype=self.dt) assert_almost_equal(norm(A), 84**0.5) assert_almost_equal(norm(A, 'fro'), 84**0.5) assert_almost_equal(norm(A, inf), 12.0) assert_almost_equal(norm(A, -inf), 4.0) assert_almost_equal(norm(A, 1), 10.0) assert_almost_equal(norm(A, -1), 6.0) assert_almost_equal(norm(A, 2), 9.1231056256176615) assert_almost_equal(norm(A, -2), 0.87689437438234041) self.assertRaises(ValueError, norm, A, 'nofro') self.assertRaises(ValueError, norm, A, -3) self.assertRaises(ValueError, norm, A, 0) def test_axis(self): # Vector norms. # Compare the use of `axis` with computing the norm of each row # or column separately. A = array([[1, 2, 3], [4, 5, 6]], dtype=self.dt) for order in [None, -1, 0, 1, 2, 3, np.Inf, -np.Inf]: expected0 = [norm(A[:,k], ord=order) for k in range(A.shape[1])] assert_almost_equal(norm(A, ord=order, axis=0), expected0) expected1 = [norm(A[k,:], ord=order) for k in range(A.shape[0])] assert_almost_equal(norm(A, ord=order, axis=1), expected1) # Matrix norms. B = np.arange(1, 25, dtype=self.dt).reshape(2, 3, 4) for order in [None, -2, 2, -1, 1, np.Inf, -np.Inf, 'fro']: assert_almost_equal(norm(A, ord=order), norm(A, ord=order, axis=(0, 1))) n = norm(B, ord=order, axis=(1, 2)) expected = [norm(B[k], ord=order) for k in range(B.shape[0])] assert_almost_equal(n, expected) n = norm(B, ord=order, axis=(2, 1)) expected = [norm(B[k].T, ord=order) for k in range(B.shape[0])] assert_almost_equal(n, expected) n = norm(B, ord=order, axis=(0, 2)) expected = [norm(B[:,k,:], ord=order) for k in range(B.shape[1])] assert_almost_equal(n, expected) n = norm(B, ord=order, axis=(0, 1)) expected = [norm(B[:,:,k], ord=order) for k in range(B.shape[2])] assert_almost_equal(n, expected) def test_bad_args(self): # Check that bad arguments raise the appropriate exceptions. A = array([[1, 2, 3], [4, 5, 6]], dtype=self.dt) B = np.arange(1, 25, dtype=self.dt).reshape(2, 3, 4) # Using `axis=` or passing in a 1-D array implies vector # norms are being computed, so also using `ord='fro'` raises a # ValueError. self.assertRaises(ValueError, norm, A, 'fro', 0) self.assertRaises(ValueError, norm, [3, 4], 'fro', None) # Similarly, norm should raise an exception when ord is any finite # number other than 1, 2, -1 or -2 when computing matrix norms. for order in [0, 3]: self.assertRaises(ValueError, norm, A, order, None) self.assertRaises(ValueError, norm, A, order, (0, 1)) self.assertRaises(ValueError, norm, B, order, (1, 2)) # Invalid axis self.assertRaises(ValueError, norm, B, None, 3) self.assertRaises(ValueError, norm, B, None, (2, 3)) self.assertRaises(ValueError, norm, B, None, (0, 1, 2)) class TestNormDouble(_TestNorm): dt = np.double dec = 12 class TestNormSingle(_TestNorm): dt = np.float32 dec = 6 class TestNormInt64(_TestNorm): dt = np.int64 dec = 12 class TestMatrixRank(object): def test_matrix_rank(self): # Full rank matrix yield assert_equal, 4, matrix_rank(np.eye(4)) # rank deficient matrix I=np.eye(4); I[-1,-1] = 0. yield assert_equal, matrix_rank(I), 3 # All zeros - zero rank yield assert_equal, matrix_rank(np.zeros((4,4))), 0 # 1 dimension - rank 1 unless all 0 yield assert_equal, matrix_rank([1, 0, 0, 0]), 1 yield assert_equal, matrix_rank(np.zeros((4,))), 0 # accepts array-like yield assert_equal, matrix_rank([1]), 1 # greater than 2 dimensions raises error yield assert_raises, TypeError, matrix_rank, np.zeros((2,2,2)) # works on scalar yield assert_equal, matrix_rank(1), 1 def test_reduced_rank(): # Test matrices with reduced rank rng = np.random.RandomState(20120714) for i in range(100): # Make a rank deficient matrix X = rng.normal(size=(40, 10)) X[:, 0] = X[:, 1] + X[:, 2] # Assert that matrix_rank detected deficiency assert_equal(matrix_rank(X), 9) X[:, 3] = X[:, 4] + X[:, 5] assert_equal(matrix_rank(X), 8) class TestQR(TestCase): def check_qr(self, a): # This test expects the argument `a` to be an ndarray or # a subclass of an ndarray of inexact type. a_type = type(a) a_dtype = a.dtype m, n = a.shape k = min(m, n) # mode == 'complete' q, r = linalg.qr(a, mode='complete') assert_(q.dtype == a_dtype) assert_(r.dtype == a_dtype) assert_(isinstance(q, a_type)) assert_(isinstance(r, a_type)) assert_(q.shape == (m, m)) assert_(r.shape == (m, n)) assert_almost_equal(dot(q, r), a) assert_almost_equal(dot(q.T.conj(), q), np.eye(m)) assert_almost_equal(np.triu(r), r) # mode == 'reduced' q1, r1 = linalg.qr(a, mode='reduced') assert_(q1.dtype == a_dtype) assert_(r1.dtype == a_dtype) assert_(isinstance(q1, a_type)) assert_(isinstance(r1, a_type)) assert_(q1.shape == (m, k)) assert_(r1.shape == (k, n)) assert_almost_equal(dot(q1, r1), a) assert_almost_equal(dot(q1.T.conj(), q1), np.eye(k)) assert_almost_equal(np.triu(r1), r1) # mode == 'r' r2 = linalg.qr(a, mode='r') assert_(r2.dtype == a_dtype) assert_(isinstance(r2, a_type)) assert_almost_equal(r2, r1) def test_qr_empty(self): a = np.zeros((0,2)) self.assertRaises(linalg.LinAlgError, linalg.qr, a) def test_mode_raw(self): a = array([[1, 2], [3, 4], [5, 6]], dtype=np.double) b = a.astype(np.single) # m > n h1, tau1 = ( array([[-5.91607978, 0.43377175, 0.72295291], [-7.43735744, 0.82807867, 0.89262383]]), array([ 1.16903085, 1.113104 ]) ) # m > n h2, tau2 = ( array([[-2.23606798, 0.61803399], [-4.91934955, -0.89442719], [-7.60263112, -1.78885438]]), array([ 1.4472136, 0. ]) ) # Test double h, tau = linalg.qr(a, mode='raw') assert_(h.dtype == np.double) assert_(tau.dtype == np.double) old_assert_almost_equal(h, h1, decimal=8) old_assert_almost_equal(tau, tau1, decimal=8) h, tau = linalg.qr(a.T, mode='raw') assert_(h.dtype == np.double) assert_(tau.dtype == np.double) old_assert_almost_equal(h, h2, decimal=8) old_assert_almost_equal(tau, tau2, decimal=8) # Test single h, tau = linalg.qr(b, mode='raw') assert_(h.dtype == np.double) assert_(tau.dtype == np.double) old_assert_almost_equal(h, h1, decimal=8) old_assert_almost_equal(tau, tau1, decimal=8) def test_mode_all_but_economic(self): a = array([[1, 2], [3, 4]]) b = array([[1, 2], [3, 4], [5, 6]]) for dt in "fd": m1 = a.astype(dt) m2 = b.astype(dt) self.check_qr(m1) self.check_qr(m2) self.check_qr(m2.T) self.check_qr(matrix(m1)) for dt in "fd": m1 = 1 + 1j * a.astype(dt) m2 = 1 + 1j * b.astype(dt) self.check_qr(m1) self.check_qr(m2) self.check_qr(m2.T) self.check_qr(matrix(m1)) def test_byteorder_check(): # Byte order check should pass for native order if sys.byteorder == 'little': native = '<' else: native = '>' for dtt in (np.float32, np.float64): arr = np.eye(4, dtype=dtt) n_arr = arr.newbyteorder(native) sw_arr = arr.newbyteorder('S').byteswap() assert_equal(arr.dtype.byteorder, '=') for routine in (linalg.inv, linalg.det, linalg.pinv): # Normal call res = routine(arr) # Native but not '=' assert_array_equal(res, routine(n_arr)) # Swapped assert_array_equal(res, routine(sw_arr)) def test_generalized_raise_multiloop(): # It should raise an error even if the error doesn't occur in the # last iteration of the ufunc inner loop invertible = np.array([[1, 2], [3, 4]]) non_invertible = np.array([[1, 1], [1, 1]]) x = np.zeros([4, 4, 2, 2])[1::2] x[...] = invertible x[0,0] = non_invertible assert_raises(np.linalg.LinAlgError, np.linalg.inv, x) def _is_xerbla_safe(): """ Check that running the xerbla test is safe --- if python_xerbla is not successfully linked in, the standard xerbla routine is called, which aborts the process. """ try: pid = os.fork() except (OSError, AttributeError): # fork failed, or not running on POSIX return False if pid == 0: # child; close i/o file handles os.close(1) os.close(0) # avoid producing core files import resource resource.setrlimit(resource.RLIMIT_CORE, (0, 0)) # these calls may abort try: a = np.array([[1]]) np.linalg.lapack_lite.dgetrf( 1, 1, a.astype(np.double), 0, # <- invalid value a.astype(np.intc), 0) except: pass try: np.linalg.lapack_lite.xerbla() except: pass os._exit(111) else: # parent pid, status = os.wait() if os.WEXITSTATUS(status) == 111 and not os.WIFSIGNALED(status): return True return False @dec.skipif(not _is_xerbla_safe(), "python_xerbla not found") def test_xerbla(): # Test that xerbla works (with GIL) a = np.array([[1]]) try: np.linalg.lapack_lite.dgetrf( 1, 1, a.astype(np.double), 0, # <- invalid value a.astype(np.intc), 0) except ValueError as e: assert_("DGETRF parameter number 4" in str(e)) else: assert_(False) # Test that xerbla works (without GIL) assert_raises(ValueError, np.linalg.lapack_lite.xerbla) if __name__ == "__main__": run_module_suite()