Chris@87: """ Test functions for linalg module Chris@87: Chris@87: """ Chris@87: from __future__ import division, absolute_import, print_function Chris@87: Chris@87: import os Chris@87: import sys Chris@87: import itertools Chris@87: import traceback Chris@87: Chris@87: import numpy as np Chris@87: from numpy import array, single, double, csingle, cdouble, dot, identity Chris@87: from numpy import multiply, atleast_2d, inf, asarray, matrix Chris@87: from numpy import linalg Chris@87: from numpy.linalg import matrix_power, norm, matrix_rank Chris@87: from numpy.testing import ( Chris@87: assert_, assert_equal, assert_raises, assert_array_equal, Chris@87: assert_almost_equal, assert_allclose, run_module_suite, Chris@87: dec Chris@87: ) Chris@87: Chris@87: Chris@87: def ifthen(a, b): Chris@87: return not a or b Chris@87: Chris@87: Chris@87: def imply(a, b): Chris@87: return not a or b Chris@87: Chris@87: Chris@87: old_assert_almost_equal = assert_almost_equal Chris@87: Chris@87: Chris@87: def assert_almost_equal(a, b, **kw): Chris@87: if asarray(a).dtype.type in (single, csingle): Chris@87: decimal = 6 Chris@87: else: Chris@87: decimal = 12 Chris@87: old_assert_almost_equal(a, b, decimal=decimal, **kw) Chris@87: Chris@87: Chris@87: def get_real_dtype(dtype): Chris@87: return {single: single, double: double, Chris@87: csingle: single, cdouble: double}[dtype] Chris@87: Chris@87: Chris@87: def get_complex_dtype(dtype): Chris@87: return {single: csingle, double: cdouble, Chris@87: csingle: csingle, cdouble: cdouble}[dtype] Chris@87: Chris@87: def get_rtol(dtype): Chris@87: # Choose a safe rtol Chris@87: if dtype in (single, csingle): Chris@87: return 1e-5 Chris@87: else: Chris@87: return 1e-11 Chris@87: Chris@87: class LinalgCase(object): Chris@87: def __init__(self, name, a, b, exception_cls=None): Chris@87: assert isinstance(name, str) Chris@87: self.name = name Chris@87: self.a = a Chris@87: self.b = b Chris@87: self.exception_cls = exception_cls Chris@87: Chris@87: def check(self, do): Chris@87: if self.exception_cls is None: Chris@87: do(self.a, self.b) Chris@87: else: Chris@87: assert_raises(self.exception_cls, do, self.a, self.b) Chris@87: Chris@87: def __repr__(self): Chris@87: return "" % (self.name,) Chris@87: Chris@87: Chris@87: # Chris@87: # Base test cases Chris@87: # Chris@87: Chris@87: np.random.seed(1234) Chris@87: Chris@87: SQUARE_CASES = [ Chris@87: LinalgCase("single", Chris@87: array([[1., 2.], [3., 4.]], dtype=single), Chris@87: array([2., 1.], dtype=single)), Chris@87: LinalgCase("double", Chris@87: array([[1., 2.], [3., 4.]], dtype=double), Chris@87: array([2., 1.], dtype=double)), Chris@87: LinalgCase("double_2", Chris@87: array([[1., 2.], [3., 4.]], dtype=double), Chris@87: array([[2., 1., 4.], [3., 4., 6.]], dtype=double)), Chris@87: LinalgCase("csingle", Chris@87: array([[1.+2j, 2+3j], [3+4j, 4+5j]], dtype=csingle), Chris@87: array([2.+1j, 1.+2j], dtype=csingle)), Chris@87: LinalgCase("cdouble", Chris@87: array([[1.+2j, 2+3j], [3+4j, 4+5j]], dtype=cdouble), Chris@87: array([2.+1j, 1.+2j], dtype=cdouble)), Chris@87: LinalgCase("cdouble_2", Chris@87: array([[1.+2j, 2+3j], [3+4j, 4+5j]], dtype=cdouble), Chris@87: array([[2.+1j, 1.+2j, 1+3j], [1-2j, 1-3j, 1-6j]], dtype=cdouble)), Chris@87: LinalgCase("empty", Chris@87: atleast_2d(array([], dtype = double)), Chris@87: atleast_2d(array([], dtype = double)), Chris@87: linalg.LinAlgError), Chris@87: LinalgCase("8x8", Chris@87: np.random.rand(8, 8), Chris@87: np.random.rand(8)), Chris@87: LinalgCase("1x1", Chris@87: np.random.rand(1, 1), Chris@87: np.random.rand(1)), Chris@87: LinalgCase("nonarray", Chris@87: [[1, 2], [3, 4]], Chris@87: [2, 1]), Chris@87: LinalgCase("matrix_b_only", Chris@87: array([[1., 2.], [3., 4.]]), Chris@87: matrix([2., 1.]).T), Chris@87: LinalgCase("matrix_a_and_b", Chris@87: matrix([[1., 2.], [3., 4.]]), Chris@87: matrix([2., 1.]).T), Chris@87: ] Chris@87: Chris@87: NONSQUARE_CASES = [ Chris@87: LinalgCase("single_nsq_1", Chris@87: array([[1., 2., 3.], [3., 4., 6.]], dtype=single), Chris@87: array([2., 1.], dtype=single)), Chris@87: LinalgCase("single_nsq_2", Chris@87: array([[1., 2.], [3., 4.], [5., 6.]], dtype=single), Chris@87: array([2., 1., 3.], dtype=single)), Chris@87: LinalgCase("double_nsq_1", Chris@87: array([[1., 2., 3.], [3., 4., 6.]], dtype=double), Chris@87: array([2., 1.], dtype=double)), Chris@87: LinalgCase("double_nsq_2", Chris@87: array([[1., 2.], [3., 4.], [5., 6.]], dtype=double), Chris@87: array([2., 1., 3.], dtype=double)), Chris@87: LinalgCase("csingle_nsq_1", Chris@87: array([[1.+1j, 2.+2j, 3.-3j], [3.-5j, 4.+9j, 6.+2j]], dtype=csingle), Chris@87: array([2.+1j, 1.+2j], dtype=csingle)), Chris@87: LinalgCase("csingle_nsq_2", Chris@87: array([[1.+1j, 2.+2j], [3.-3j, 4.-9j], [5.-4j, 6.+8j]], dtype=csingle), Chris@87: array([2.+1j, 1.+2j, 3.-3j], dtype=csingle)), Chris@87: LinalgCase("cdouble_nsq_1", Chris@87: array([[1.+1j, 2.+2j, 3.-3j], [3.-5j, 4.+9j, 6.+2j]], dtype=cdouble), Chris@87: array([2.+1j, 1.+2j], dtype=cdouble)), Chris@87: LinalgCase("cdouble_nsq_2", Chris@87: array([[1.+1j, 2.+2j], [3.-3j, 4.-9j], [5.-4j, 6.+8j]], dtype=cdouble), Chris@87: array([2.+1j, 1.+2j, 3.-3j], dtype=cdouble)), Chris@87: LinalgCase("cdouble_nsq_1_2", Chris@87: array([[1.+1j, 2.+2j, 3.-3j], [3.-5j, 4.+9j, 6.+2j]], dtype=cdouble), Chris@87: array([[2.+1j, 1.+2j], [1-1j, 2-2j]], dtype=cdouble)), Chris@87: LinalgCase("cdouble_nsq_2_2", Chris@87: array([[1.+1j, 2.+2j], [3.-3j, 4.-9j], [5.-4j, 6.+8j]], dtype=cdouble), Chris@87: array([[2.+1j, 1.+2j], [1-1j, 2-2j], [1-1j, 2-2j]], dtype=cdouble)), Chris@87: LinalgCase("8x11", Chris@87: np.random.rand(8, 11), Chris@87: np.random.rand(11)), Chris@87: LinalgCase("1x5", Chris@87: np.random.rand(1, 5), Chris@87: np.random.rand(5)), Chris@87: LinalgCase("5x1", Chris@87: np.random.rand(5, 1), Chris@87: np.random.rand(1)), Chris@87: ] Chris@87: Chris@87: HERMITIAN_CASES = [ Chris@87: LinalgCase("hsingle", Chris@87: array([[1., 2.], [2., 1.]], dtype=single), Chris@87: None), Chris@87: LinalgCase("hdouble", Chris@87: array([[1., 2.], [2., 1.]], dtype=double), Chris@87: None), Chris@87: LinalgCase("hcsingle", Chris@87: array([[1., 2+3j], [2-3j, 1]], dtype=csingle), Chris@87: None), Chris@87: LinalgCase("hcdouble", Chris@87: array([[1., 2+3j], [2-3j, 1]], dtype=cdouble), Chris@87: None), Chris@87: LinalgCase("hempty", Chris@87: atleast_2d(array([], dtype = double)), Chris@87: None, Chris@87: linalg.LinAlgError), Chris@87: LinalgCase("hnonarray", Chris@87: [[1, 2], [2, 1]], Chris@87: None), Chris@87: LinalgCase("matrix_b_only", Chris@87: array([[1., 2.], [2., 1.]]), Chris@87: None), Chris@87: LinalgCase("hmatrix_a_and_b", Chris@87: matrix([[1., 2.], [2., 1.]]), Chris@87: None), Chris@87: LinalgCase("hmatrix_1x1", Chris@87: np.random.rand(1, 1), Chris@87: None), Chris@87: ] Chris@87: Chris@87: Chris@87: # Chris@87: # Gufunc test cases Chris@87: # Chris@87: Chris@87: GENERALIZED_SQUARE_CASES = [] Chris@87: GENERALIZED_NONSQUARE_CASES = [] Chris@87: GENERALIZED_HERMITIAN_CASES = [] Chris@87: Chris@87: for tgt, src in ((GENERALIZED_SQUARE_CASES, SQUARE_CASES), Chris@87: (GENERALIZED_NONSQUARE_CASES, NONSQUARE_CASES), Chris@87: (GENERALIZED_HERMITIAN_CASES, HERMITIAN_CASES)): Chris@87: for case in src: Chris@87: if not isinstance(case.a, np.ndarray): Chris@87: continue Chris@87: Chris@87: a = np.array([case.a, 2*case.a, 3*case.a]) Chris@87: if case.b is None: Chris@87: b = None Chris@87: else: Chris@87: b = np.array([case.b, 7*case.b, 6*case.b]) Chris@87: new_case = LinalgCase(case.name + "_tile3", a, b, Chris@87: case.exception_cls) Chris@87: tgt.append(new_case) Chris@87: Chris@87: a = np.array([case.a]*2*3).reshape((3, 2) + case.a.shape) Chris@87: if case.b is None: Chris@87: b = None Chris@87: else: Chris@87: b = np.array([case.b]*2*3).reshape((3, 2) + case.b.shape) Chris@87: new_case = LinalgCase(case.name + "_tile213", a, b, Chris@87: case.exception_cls) Chris@87: tgt.append(new_case) Chris@87: Chris@87: # Chris@87: # Generate stride combination variations of the above Chris@87: # Chris@87: Chris@87: def _stride_comb_iter(x): Chris@87: """ Chris@87: Generate cartesian product of strides for all axes Chris@87: """ Chris@87: Chris@87: if not isinstance(x, np.ndarray): Chris@87: yield x, "nop" Chris@87: return Chris@87: Chris@87: stride_set = [(1,)]*x.ndim Chris@87: stride_set[-1] = (1, 3, -4) Chris@87: if x.ndim > 1: Chris@87: stride_set[-2] = (1, 3, -4) Chris@87: if x.ndim > 2: Chris@87: stride_set[-3] = (1, -4) Chris@87: Chris@87: for repeats in itertools.product(*tuple(stride_set)): Chris@87: new_shape = [abs(a*b) for a, b in zip(x.shape, repeats)] Chris@87: slices = tuple([slice(None, None, repeat) for repeat in repeats]) Chris@87: Chris@87: # new array with different strides, but same data Chris@87: xi = np.empty(new_shape, dtype=x.dtype) Chris@87: xi.view(np.uint32).fill(0xdeadbeef) Chris@87: xi = xi[slices] Chris@87: xi[...] = x Chris@87: xi = xi.view(x.__class__) Chris@87: assert np.all(xi == x) Chris@87: yield xi, "stride_" + "_".join(["%+d" % j for j in repeats]) Chris@87: Chris@87: # generate also zero strides if possible Chris@87: if x.ndim >= 1 and x.shape[-1] == 1: Chris@87: s = list(x.strides) Chris@87: s[-1] = 0 Chris@87: xi = np.lib.stride_tricks.as_strided(x, strides=s) Chris@87: yield xi, "stride_xxx_0" Chris@87: if x.ndim >= 2 and x.shape[-2] == 1: Chris@87: s = list(x.strides) Chris@87: s[-2] = 0 Chris@87: xi = np.lib.stride_tricks.as_strided(x, strides=s) Chris@87: yield xi, "stride_xxx_0_x" Chris@87: if x.ndim >= 2 and x.shape[:-2] == (1, 1): Chris@87: s = list(x.strides) Chris@87: s[-1] = 0 Chris@87: s[-2] = 0 Chris@87: xi = np.lib.stride_tricks.as_strided(x, strides=s) Chris@87: yield xi, "stride_xxx_0_0" Chris@87: Chris@87: for src in (SQUARE_CASES, Chris@87: NONSQUARE_CASES, Chris@87: HERMITIAN_CASES, Chris@87: GENERALIZED_SQUARE_CASES, Chris@87: GENERALIZED_NONSQUARE_CASES, Chris@87: GENERALIZED_HERMITIAN_CASES): Chris@87: Chris@87: new_cases = [] Chris@87: for case in src: Chris@87: for a, a_tag in _stride_comb_iter(case.a): Chris@87: for b, b_tag in _stride_comb_iter(case.b): Chris@87: new_case = LinalgCase(case.name + "_" + a_tag + "_" + b_tag, a, b, Chris@87: exception_cls=case.exception_cls) Chris@87: new_cases.append(new_case) Chris@87: src.extend(new_cases) Chris@87: Chris@87: Chris@87: # Chris@87: # Test different routines against the above cases Chris@87: # Chris@87: Chris@87: def _check_cases(func, cases): Chris@87: for case in cases: Chris@87: try: Chris@87: case.check(func) Chris@87: except Exception: Chris@87: msg = "In test case: %r\n\n" % case Chris@87: msg += traceback.format_exc() Chris@87: raise AssertionError(msg) Chris@87: Chris@87: class LinalgTestCase(object): Chris@87: def test_sq_cases(self): Chris@87: _check_cases(self.do, SQUARE_CASES) Chris@87: Chris@87: Chris@87: class LinalgNonsquareTestCase(object): Chris@87: def test_sq_cases(self): Chris@87: _check_cases(self.do, NONSQUARE_CASES) Chris@87: Chris@87: Chris@87: class LinalgGeneralizedTestCase(object): Chris@87: @dec.slow Chris@87: def test_generalized_sq_cases(self): Chris@87: _check_cases(self.do, GENERALIZED_SQUARE_CASES) Chris@87: Chris@87: Chris@87: class LinalgGeneralizedNonsquareTestCase(object): Chris@87: @dec.slow Chris@87: def test_generalized_nonsq_cases(self): Chris@87: _check_cases(self.do, GENERALIZED_NONSQUARE_CASES) Chris@87: Chris@87: Chris@87: class HermitianTestCase(object): Chris@87: def test_herm_cases(self): Chris@87: _check_cases(self.do, HERMITIAN_CASES) Chris@87: Chris@87: Chris@87: class HermitianGeneralizedTestCase(object): Chris@87: @dec.slow Chris@87: def test_generalized_herm_cases(self): Chris@87: _check_cases(self.do, GENERALIZED_HERMITIAN_CASES) Chris@87: Chris@87: Chris@87: def dot_generalized(a, b): Chris@87: a = asarray(a) Chris@87: if a.ndim >= 3: Chris@87: if a.ndim == b.ndim: Chris@87: # matrix x matrix Chris@87: new_shape = a.shape[:-1] + b.shape[-1:] Chris@87: elif a.ndim == b.ndim + 1: Chris@87: # matrix x vector Chris@87: new_shape = a.shape[:-1] Chris@87: else: Chris@87: raise ValueError("Not implemented...") Chris@87: r = np.empty(new_shape, dtype=np.common_type(a, b)) Chris@87: for c in itertools.product(*map(range, a.shape[:-2])): Chris@87: r[c] = dot(a[c], b[c]) Chris@87: return r Chris@87: else: Chris@87: return dot(a, b) Chris@87: Chris@87: Chris@87: def identity_like_generalized(a): Chris@87: a = asarray(a) Chris@87: if a.ndim >= 3: Chris@87: r = np.empty(a.shape, dtype=a.dtype) Chris@87: for c in itertools.product(*map(range, a.shape[:-2])): Chris@87: r[c] = identity(a.shape[-2]) Chris@87: return r Chris@87: else: Chris@87: return identity(a.shape[0]) Chris@87: Chris@87: Chris@87: class TestSolve(LinalgTestCase, LinalgGeneralizedTestCase): Chris@87: def do(self, a, b): Chris@87: x = linalg.solve(a, b) Chris@87: assert_almost_equal(b, dot_generalized(a, x)) Chris@87: assert_(imply(isinstance(b, matrix), isinstance(x, matrix))) Chris@87: Chris@87: def test_types(self): Chris@87: def check(dtype): Chris@87: x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype) Chris@87: assert_equal(linalg.solve(x, x).dtype, dtype) Chris@87: for dtype in [single, double, csingle, cdouble]: Chris@87: yield check, dtype Chris@87: Chris@87: def test_0_size(self): Chris@87: class ArraySubclass(np.ndarray): Chris@87: pass Chris@87: # Test system of 0x0 matrices Chris@87: a = np.arange(8).reshape(2, 2, 2) Chris@87: b = np.arange(6).reshape(1, 2, 3).view(ArraySubclass) Chris@87: Chris@87: expected = linalg.solve(a, b)[:, 0:0,:] Chris@87: result = linalg.solve(a[:, 0:0, 0:0], b[:, 0:0,:]) Chris@87: assert_array_equal(result, expected) Chris@87: assert_(isinstance(result, ArraySubclass)) Chris@87: Chris@87: # Test errors for non-square and only b's dimension being 0 Chris@87: assert_raises(linalg.LinAlgError, linalg.solve, a[:, 0:0, 0:1], b) Chris@87: assert_raises(ValueError, linalg.solve, a, b[:, 0:0,:]) Chris@87: Chris@87: # Test broadcasting error Chris@87: b = np.arange(6).reshape(1, 3, 2) # broadcasting error Chris@87: assert_raises(ValueError, linalg.solve, a, b) Chris@87: assert_raises(ValueError, linalg.solve, a[0:0], b[0:0]) Chris@87: Chris@87: # Test zero "single equations" with 0x0 matrices. Chris@87: b = np.arange(2).reshape(1, 2).view(ArraySubclass) Chris@87: expected = linalg.solve(a, b)[:, 0:0] Chris@87: result = linalg.solve(a[:, 0:0, 0:0], b[:, 0:0]) Chris@87: assert_array_equal(result, expected) Chris@87: assert_(isinstance(result, ArraySubclass)) Chris@87: Chris@87: b = np.arange(3).reshape(1, 3) Chris@87: assert_raises(ValueError, linalg.solve, a, b) Chris@87: assert_raises(ValueError, linalg.solve, a[0:0], b[0:0]) Chris@87: assert_raises(ValueError, linalg.solve, a[:, 0:0, 0:0], b) Chris@87: Chris@87: def test_0_size_k(self): Chris@87: # test zero multiple equation (K=0) case. Chris@87: class ArraySubclass(np.ndarray): Chris@87: pass Chris@87: a = np.arange(4).reshape(1, 2, 2) Chris@87: b = np.arange(6).reshape(3, 2, 1).view(ArraySubclass) Chris@87: Chris@87: expected = linalg.solve(a, b)[:,:, 0:0] Chris@87: result = linalg.solve(a, b[:,:, 0:0]) Chris@87: assert_array_equal(result, expected) Chris@87: assert_(isinstance(result, ArraySubclass)) Chris@87: Chris@87: # test both zero. Chris@87: expected = linalg.solve(a, b)[:, 0:0, 0:0] Chris@87: result = linalg.solve(a[:, 0:0, 0:0], b[:,0:0, 0:0]) Chris@87: assert_array_equal(result, expected) Chris@87: assert_(isinstance(result, ArraySubclass)) Chris@87: Chris@87: Chris@87: class TestInv(LinalgTestCase, LinalgGeneralizedTestCase): Chris@87: def do(self, a, b): Chris@87: a_inv = linalg.inv(a) Chris@87: assert_almost_equal(dot_generalized(a, a_inv), Chris@87: identity_like_generalized(a)) Chris@87: assert_(imply(isinstance(a, matrix), isinstance(a_inv, matrix))) Chris@87: Chris@87: def test_types(self): Chris@87: def check(dtype): Chris@87: x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype) Chris@87: assert_equal(linalg.inv(x).dtype, dtype) Chris@87: for dtype in [single, double, csingle, cdouble]: Chris@87: yield check, dtype Chris@87: Chris@87: def test_0_size(self): Chris@87: # Check that all kinds of 0-sized arrays work Chris@87: class ArraySubclass(np.ndarray): Chris@87: pass Chris@87: a = np.zeros((0, 1, 1), dtype=np.int_).view(ArraySubclass) Chris@87: res = linalg.inv(a) Chris@87: assert_(res.dtype.type is np.float64) Chris@87: assert_equal(a.shape, res.shape) Chris@87: assert_(isinstance(a, ArraySubclass)) Chris@87: Chris@87: a = np.zeros((0, 0), dtype=np.complex64).view(ArraySubclass) Chris@87: res = linalg.inv(a) Chris@87: assert_(res.dtype.type is np.complex64) Chris@87: assert_equal(a.shape, res.shape) Chris@87: Chris@87: Chris@87: class TestEigvals(LinalgTestCase, LinalgGeneralizedTestCase): Chris@87: def do(self, a, b): Chris@87: ev = linalg.eigvals(a) Chris@87: evalues, evectors = linalg.eig(a) Chris@87: assert_almost_equal(ev, evalues) Chris@87: Chris@87: def test_types(self): Chris@87: def check(dtype): Chris@87: x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype) Chris@87: assert_equal(linalg.eigvals(x).dtype, dtype) Chris@87: x = np.array([[1, 0.5], [-1, 1]], dtype=dtype) Chris@87: assert_equal(linalg.eigvals(x).dtype, get_complex_dtype(dtype)) Chris@87: for dtype in [single, double, csingle, cdouble]: Chris@87: yield check, dtype Chris@87: Chris@87: Chris@87: class TestEig(LinalgTestCase, LinalgGeneralizedTestCase): Chris@87: def do(self, a, b): Chris@87: evalues, evectors = linalg.eig(a) Chris@87: assert_allclose(dot_generalized(a, evectors), Chris@87: np.asarray(evectors) * np.asarray(evalues)[...,None,:], Chris@87: rtol=get_rtol(evalues.dtype)) Chris@87: assert_(imply(isinstance(a, matrix), isinstance(evectors, matrix))) Chris@87: Chris@87: def test_types(self): Chris@87: def check(dtype): Chris@87: x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype) Chris@87: w, v = np.linalg.eig(x) Chris@87: assert_equal(w.dtype, dtype) Chris@87: assert_equal(v.dtype, dtype) Chris@87: Chris@87: x = np.array([[1, 0.5], [-1, 1]], dtype=dtype) Chris@87: w, v = np.linalg.eig(x) Chris@87: assert_equal(w.dtype, get_complex_dtype(dtype)) Chris@87: assert_equal(v.dtype, get_complex_dtype(dtype)) Chris@87: Chris@87: for dtype in [single, double, csingle, cdouble]: Chris@87: yield check, dtype Chris@87: Chris@87: Chris@87: class TestSVD(LinalgTestCase, LinalgGeneralizedTestCase): Chris@87: def do(self, a, b): Chris@87: u, s, vt = linalg.svd(a, 0) Chris@87: assert_allclose(a, dot_generalized(np.asarray(u) * np.asarray(s)[...,None,:], Chris@87: np.asarray(vt)), Chris@87: rtol=get_rtol(u.dtype)) Chris@87: assert_(imply(isinstance(a, matrix), isinstance(u, matrix))) Chris@87: assert_(imply(isinstance(a, matrix), isinstance(vt, matrix))) Chris@87: Chris@87: def test_types(self): Chris@87: def check(dtype): Chris@87: x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype) Chris@87: u, s, vh = linalg.svd(x) Chris@87: assert_equal(u.dtype, dtype) Chris@87: assert_equal(s.dtype, get_real_dtype(dtype)) Chris@87: assert_equal(vh.dtype, dtype) Chris@87: s = linalg.svd(x, compute_uv=False) Chris@87: assert_equal(s.dtype, get_real_dtype(dtype)) Chris@87: Chris@87: for dtype in [single, double, csingle, cdouble]: Chris@87: yield check, dtype Chris@87: Chris@87: Chris@87: class TestCondSVD(LinalgTestCase, LinalgGeneralizedTestCase): Chris@87: def do(self, a, b): Chris@87: c = asarray(a) # a might be a matrix Chris@87: s = linalg.svd(c, compute_uv=False) Chris@87: old_assert_almost_equal(s[0]/s[-1], linalg.cond(a), decimal=5) Chris@87: Chris@87: Chris@87: class TestCond2(LinalgTestCase): Chris@87: def do(self, a, b): Chris@87: c = asarray(a) # a might be a matrix Chris@87: s = linalg.svd(c, compute_uv=False) Chris@87: old_assert_almost_equal(s[0]/s[-1], linalg.cond(a, 2), decimal=5) Chris@87: Chris@87: Chris@87: class TestCondInf(object): Chris@87: def test(self): Chris@87: A = array([[1., 0, 0], [0, -2., 0], [0, 0, 3.]]) Chris@87: assert_almost_equal(linalg.cond(A, inf), 3.) Chris@87: Chris@87: Chris@87: class TestPinv(LinalgTestCase): Chris@87: def do(self, a, b): Chris@87: a_ginv = linalg.pinv(a) Chris@87: assert_almost_equal(dot(a, a_ginv), identity(asarray(a).shape[0])) Chris@87: assert_(imply(isinstance(a, matrix), isinstance(a_ginv, matrix))) Chris@87: Chris@87: Chris@87: class TestDet(LinalgTestCase, LinalgGeneralizedTestCase): Chris@87: def do(self, a, b): Chris@87: d = linalg.det(a) Chris@87: (s, ld) = linalg.slogdet(a) Chris@87: if asarray(a).dtype.type in (single, double): Chris@87: ad = asarray(a).astype(double) Chris@87: else: Chris@87: ad = asarray(a).astype(cdouble) Chris@87: ev = linalg.eigvals(ad) Chris@87: assert_almost_equal(d, multiply.reduce(ev, axis=-1)) Chris@87: assert_almost_equal(s * np.exp(ld), multiply.reduce(ev, axis=-1)) Chris@87: Chris@87: s = np.atleast_1d(s) Chris@87: ld = np.atleast_1d(ld) Chris@87: m = (s != 0) Chris@87: assert_almost_equal(np.abs(s[m]), 1) Chris@87: assert_equal(ld[~m], -inf) Chris@87: Chris@87: def test_zero(self): Chris@87: assert_equal(linalg.det([[0.0]]), 0.0) Chris@87: assert_equal(type(linalg.det([[0.0]])), double) Chris@87: assert_equal(linalg.det([[0.0j]]), 0.0) Chris@87: assert_equal(type(linalg.det([[0.0j]])), cdouble) Chris@87: Chris@87: assert_equal(linalg.slogdet([[0.0]]), (0.0, -inf)) Chris@87: assert_equal(type(linalg.slogdet([[0.0]])[0]), double) Chris@87: assert_equal(type(linalg.slogdet([[0.0]])[1]), double) Chris@87: assert_equal(linalg.slogdet([[0.0j]]), (0.0j, -inf)) Chris@87: assert_equal(type(linalg.slogdet([[0.0j]])[0]), cdouble) Chris@87: assert_equal(type(linalg.slogdet([[0.0j]])[1]), double) Chris@87: Chris@87: def test_types(self): Chris@87: def check(dtype): Chris@87: x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype) Chris@87: assert_equal(np.linalg.det(x).dtype, dtype) Chris@87: ph, s = np.linalg.slogdet(x) Chris@87: assert_equal(s.dtype, get_real_dtype(dtype)) Chris@87: assert_equal(ph.dtype, dtype) Chris@87: for dtype in [single, double, csingle, cdouble]: Chris@87: yield check, dtype Chris@87: Chris@87: Chris@87: class TestLstsq(LinalgTestCase, LinalgNonsquareTestCase): Chris@87: def do(self, a, b): Chris@87: arr = np.asarray(a) Chris@87: m, n = arr.shape Chris@87: u, s, vt = linalg.svd(a, 0) Chris@87: x, residuals, rank, sv = linalg.lstsq(a, b) Chris@87: if m <= n: Chris@87: assert_almost_equal(b, dot(a, x)) Chris@87: assert_equal(rank, m) Chris@87: else: Chris@87: assert_equal(rank, n) Chris@87: assert_almost_equal(sv, sv.__array_wrap__(s)) Chris@87: if rank == n and m > n: Chris@87: expect_resids = (np.asarray(abs(np.dot(a, x) - b))**2).sum(axis=0) Chris@87: expect_resids = np.asarray(expect_resids) Chris@87: if len(np.asarray(b).shape) == 1: Chris@87: expect_resids.shape = (1,) Chris@87: assert_equal(residuals.shape, expect_resids.shape) Chris@87: else: Chris@87: expect_resids = np.array([]).view(type(x)) Chris@87: assert_almost_equal(residuals, expect_resids) Chris@87: assert_(np.issubdtype(residuals.dtype, np.floating)) Chris@87: assert_(imply(isinstance(b, matrix), isinstance(x, matrix))) Chris@87: assert_(imply(isinstance(b, matrix), isinstance(residuals, matrix))) Chris@87: Chris@87: Chris@87: class TestMatrixPower(object): Chris@87: R90 = array([[0, 1], [-1, 0]]) Chris@87: Arb22 = array([[4, -7], [-2, 10]]) Chris@87: noninv = array([[1, 0], [0, 0]]) Chris@87: arbfloat = array([[0.1, 3.2], [1.2, 0.7]]) Chris@87: Chris@87: large = identity(10) Chris@87: t = large[1,:].copy() Chris@87: large[1,:] = large[0,:] Chris@87: large[0,:] = t Chris@87: Chris@87: def test_large_power(self): Chris@87: assert_equal(matrix_power(self.R90, 2**100+2**10+2**5+1), self.R90) Chris@87: Chris@87: def test_large_power_trailing_zero(self): Chris@87: assert_equal(matrix_power(self.R90, 2**100+2**10+2**5), identity(2)) Chris@87: Chris@87: def testip_zero(self): Chris@87: def tz(M): Chris@87: mz = matrix_power(M, 0) Chris@87: assert_equal(mz, identity(M.shape[0])) Chris@87: assert_equal(mz.dtype, M.dtype) Chris@87: for M in [self.Arb22, self.arbfloat, self.large]: Chris@87: yield tz, M Chris@87: Chris@87: def testip_one(self): Chris@87: def tz(M): Chris@87: mz = matrix_power(M, 1) Chris@87: assert_equal(mz, M) Chris@87: assert_equal(mz.dtype, M.dtype) Chris@87: for M in [self.Arb22, self.arbfloat, self.large]: Chris@87: yield tz, M Chris@87: Chris@87: def testip_two(self): Chris@87: def tz(M): Chris@87: mz = matrix_power(M, 2) Chris@87: assert_equal(mz, dot(M, M)) Chris@87: assert_equal(mz.dtype, M.dtype) Chris@87: for M in [self.Arb22, self.arbfloat, self.large]: Chris@87: yield tz, M Chris@87: Chris@87: def testip_invert(self): Chris@87: def tz(M): Chris@87: mz = matrix_power(M, -1) Chris@87: assert_almost_equal(identity(M.shape[0]), dot(mz, M)) Chris@87: for M in [self.R90, self.Arb22, self.arbfloat, self.large]: Chris@87: yield tz, M Chris@87: Chris@87: def test_invert_noninvertible(self): Chris@87: import numpy.linalg Chris@87: assert_raises(numpy.linalg.linalg.LinAlgError, Chris@87: lambda: matrix_power(self.noninv, -1)) Chris@87: Chris@87: Chris@87: class TestBoolPower(object): Chris@87: def test_square(self): Chris@87: A = array([[True, False], [True, True]]) Chris@87: assert_equal(matrix_power(A, 2), A) Chris@87: Chris@87: Chris@87: class TestEigvalsh(HermitianTestCase, HermitianGeneralizedTestCase): Chris@87: def do(self, a, b): Chris@87: # note that eigenvalue arrays must be sorted since Chris@87: # their order isn't guaranteed. Chris@87: ev = linalg.eigvalsh(a, 'L') Chris@87: evalues, evectors = linalg.eig(a) Chris@87: ev.sort(axis=-1) Chris@87: evalues.sort(axis=-1) Chris@87: assert_allclose(ev, evalues, Chris@87: rtol=get_rtol(ev.dtype)) Chris@87: Chris@87: ev2 = linalg.eigvalsh(a, 'U') Chris@87: ev2.sort(axis=-1) Chris@87: assert_allclose(ev2, evalues, Chris@87: rtol=get_rtol(ev.dtype)) Chris@87: Chris@87: def test_types(self): Chris@87: def check(dtype): Chris@87: x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype) Chris@87: w = np.linalg.eigvalsh(x) Chris@87: assert_equal(w.dtype, get_real_dtype(dtype)) Chris@87: for dtype in [single, double, csingle, cdouble]: Chris@87: yield check, dtype Chris@87: Chris@87: def test_invalid(self): Chris@87: x = np.array([[1, 0.5], [0.5, 1]], dtype=np.float32) Chris@87: assert_raises(ValueError, np.linalg.eigvalsh, x, UPLO="lrong") Chris@87: assert_raises(ValueError, np.linalg.eigvalsh, x, "lower") Chris@87: assert_raises(ValueError, np.linalg.eigvalsh, x, "upper") Chris@87: Chris@87: def test_UPLO(self): Chris@87: Klo = np.array([[0, 0],[1, 0]], dtype=np.double) Chris@87: Kup = np.array([[0, 1],[0, 0]], dtype=np.double) Chris@87: tgt = np.array([-1, 1], dtype=np.double) Chris@87: rtol = get_rtol(np.double) Chris@87: Chris@87: # Check default is 'L' Chris@87: w = np.linalg.eigvalsh(Klo) Chris@87: assert_allclose(np.sort(w), tgt, rtol=rtol) Chris@87: # Check 'L' Chris@87: w = np.linalg.eigvalsh(Klo, UPLO='L') Chris@87: assert_allclose(np.sort(w), tgt, rtol=rtol) Chris@87: # Check 'l' Chris@87: w = np.linalg.eigvalsh(Klo, UPLO='l') Chris@87: assert_allclose(np.sort(w), tgt, rtol=rtol) Chris@87: # Check 'U' Chris@87: w = np.linalg.eigvalsh(Kup, UPLO='U') Chris@87: assert_allclose(np.sort(w), tgt, rtol=rtol) Chris@87: # Check 'u' Chris@87: w = np.linalg.eigvalsh(Kup, UPLO='u') Chris@87: assert_allclose(np.sort(w), tgt, rtol=rtol) Chris@87: Chris@87: Chris@87: class TestEigh(HermitianTestCase, HermitianGeneralizedTestCase): Chris@87: def do(self, a, b): Chris@87: # note that eigenvalue arrays must be sorted since Chris@87: # their order isn't guaranteed. Chris@87: ev, evc = linalg.eigh(a) Chris@87: evalues, evectors = linalg.eig(a) Chris@87: ev.sort(axis=-1) Chris@87: evalues.sort(axis=-1) Chris@87: assert_almost_equal(ev, evalues) Chris@87: Chris@87: assert_allclose(dot_generalized(a, evc), Chris@87: np.asarray(ev)[...,None,:] * np.asarray(evc), Chris@87: rtol=get_rtol(ev.dtype)) Chris@87: Chris@87: ev2, evc2 = linalg.eigh(a, 'U') Chris@87: ev2.sort(axis=-1) Chris@87: assert_almost_equal(ev2, evalues) Chris@87: Chris@87: assert_allclose(dot_generalized(a, evc2), Chris@87: np.asarray(ev2)[...,None,:] * np.asarray(evc2), Chris@87: rtol=get_rtol(ev.dtype), err_msg=repr(a)) Chris@87: Chris@87: def test_types(self): Chris@87: def check(dtype): Chris@87: x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype) Chris@87: w, v = np.linalg.eigh(x) Chris@87: assert_equal(w.dtype, get_real_dtype(dtype)) Chris@87: assert_equal(v.dtype, dtype) Chris@87: for dtype in [single, double, csingle, cdouble]: Chris@87: yield check, dtype Chris@87: Chris@87: def test_invalid(self): Chris@87: x = np.array([[1, 0.5], [0.5, 1]], dtype=np.float32) Chris@87: assert_raises(ValueError, np.linalg.eigh, x, UPLO="lrong") Chris@87: assert_raises(ValueError, np.linalg.eigh, x, "lower") Chris@87: assert_raises(ValueError, np.linalg.eigh, x, "upper") Chris@87: Chris@87: def test_UPLO(self): Chris@87: Klo = np.array([[0, 0],[1, 0]], dtype=np.double) Chris@87: Kup = np.array([[0, 1],[0, 0]], dtype=np.double) Chris@87: tgt = np.array([-1, 1], dtype=np.double) Chris@87: rtol = get_rtol(np.double) Chris@87: Chris@87: # Check default is 'L' Chris@87: w, v = np.linalg.eigh(Klo) Chris@87: assert_allclose(np.sort(w), tgt, rtol=rtol) Chris@87: # Check 'L' Chris@87: w, v = np.linalg.eigh(Klo, UPLO='L') Chris@87: assert_allclose(np.sort(w), tgt, rtol=rtol) Chris@87: # Check 'l' Chris@87: w, v = np.linalg.eigh(Klo, UPLO='l') Chris@87: assert_allclose(np.sort(w), tgt, rtol=rtol) Chris@87: # Check 'U' Chris@87: w, v = np.linalg.eigh(Kup, UPLO='U') Chris@87: assert_allclose(np.sort(w), tgt, rtol=rtol) Chris@87: # Check 'u' Chris@87: w, v = np.linalg.eigh(Kup, UPLO='u') Chris@87: assert_allclose(np.sort(w), tgt, rtol=rtol) Chris@87: Chris@87: Chris@87: class _TestNorm(object): Chris@87: Chris@87: dt = None Chris@87: dec = None Chris@87: Chris@87: def test_empty(self): Chris@87: assert_equal(norm([]), 0.0) Chris@87: assert_equal(norm(array([], dtype=self.dt)), 0.0) Chris@87: assert_equal(norm(atleast_2d(array([], dtype=self.dt))), 0.0) Chris@87: Chris@87: def test_vector(self): Chris@87: a = [1, 2, 3, 4] Chris@87: b = [-1, -2, -3, -4] Chris@87: c = [-1, 2, -3, 4] Chris@87: Chris@87: def _test(v): Chris@87: np.testing.assert_almost_equal(norm(v), 30**0.5, Chris@87: decimal=self.dec) Chris@87: np.testing.assert_almost_equal(norm(v, inf), 4.0, Chris@87: decimal=self.dec) Chris@87: np.testing.assert_almost_equal(norm(v, -inf), 1.0, Chris@87: decimal=self.dec) Chris@87: np.testing.assert_almost_equal(norm(v, 1), 10.0, Chris@87: decimal=self.dec) Chris@87: np.testing.assert_almost_equal(norm(v, -1), 12.0/25, Chris@87: decimal=self.dec) Chris@87: np.testing.assert_almost_equal(norm(v, 2), 30**0.5, Chris@87: decimal=self.dec) Chris@87: np.testing.assert_almost_equal(norm(v, -2), ((205./144)**-0.5), Chris@87: decimal=self.dec) Chris@87: np.testing.assert_almost_equal(norm(v, 0), 4, Chris@87: decimal=self.dec) Chris@87: Chris@87: for v in (a, b, c,): Chris@87: _test(v) Chris@87: Chris@87: for v in (array(a, dtype=self.dt), array(b, dtype=self.dt), Chris@87: array(c, dtype=self.dt)): Chris@87: _test(v) Chris@87: Chris@87: def test_matrix(self): Chris@87: A = matrix([[1, 3], [5, 7]], dtype=self.dt) Chris@87: assert_almost_equal(norm(A), 84**0.5) Chris@87: assert_almost_equal(norm(A, 'fro'), 84**0.5) Chris@87: assert_almost_equal(norm(A, inf), 12.0) Chris@87: assert_almost_equal(norm(A, -inf), 4.0) Chris@87: assert_almost_equal(norm(A, 1), 10.0) Chris@87: assert_almost_equal(norm(A, -1), 6.0) Chris@87: assert_almost_equal(norm(A, 2), 9.1231056256176615) Chris@87: assert_almost_equal(norm(A, -2), 0.87689437438234041) Chris@87: Chris@87: assert_raises(ValueError, norm, A, 'nofro') Chris@87: assert_raises(ValueError, norm, A, -3) Chris@87: assert_raises(ValueError, norm, A, 0) Chris@87: Chris@87: def test_axis(self): Chris@87: # Vector norms. Chris@87: # Compare the use of `axis` with computing the norm of each row Chris@87: # or column separately. Chris@87: A = array([[1, 2, 3], [4, 5, 6]], dtype=self.dt) Chris@87: for order in [None, -1, 0, 1, 2, 3, np.Inf, -np.Inf]: Chris@87: expected0 = [norm(A[:, k], ord=order) for k in range(A.shape[1])] Chris@87: assert_almost_equal(norm(A, ord=order, axis=0), expected0) Chris@87: expected1 = [norm(A[k,:], ord=order) for k in range(A.shape[0])] Chris@87: assert_almost_equal(norm(A, ord=order, axis=1), expected1) Chris@87: Chris@87: # Matrix norms. Chris@87: B = np.arange(1, 25, dtype=self.dt).reshape(2, 3, 4) Chris@87: Chris@87: for order in [None, -2, 2, -1, 1, np.Inf, -np.Inf, 'fro']: Chris@87: assert_almost_equal(norm(A, ord=order), norm(A, ord=order, Chris@87: axis=(0, 1))) Chris@87: Chris@87: n = norm(B, ord=order, axis=(1, 2)) Chris@87: expected = [norm(B[k], ord=order) for k in range(B.shape[0])] Chris@87: assert_almost_equal(n, expected) Chris@87: Chris@87: n = norm(B, ord=order, axis=(2, 1)) Chris@87: expected = [norm(B[k].T, ord=order) for k in range(B.shape[0])] Chris@87: assert_almost_equal(n, expected) Chris@87: Chris@87: n = norm(B, ord=order, axis=(0, 2)) Chris@87: expected = [norm(B[:, k,:], ord=order) for k in range(B.shape[1])] Chris@87: assert_almost_equal(n, expected) Chris@87: Chris@87: n = norm(B, ord=order, axis=(0, 1)) Chris@87: expected = [norm(B[:,:, k], ord=order) for k in range(B.shape[2])] Chris@87: assert_almost_equal(n, expected) Chris@87: Chris@87: def test_bad_args(self): Chris@87: # Check that bad arguments raise the appropriate exceptions. Chris@87: Chris@87: A = array([[1, 2, 3], [4, 5, 6]], dtype=self.dt) Chris@87: B = np.arange(1, 25, dtype=self.dt).reshape(2, 3, 4) Chris@87: Chris@87: # Using `axis=` or passing in a 1-D array implies vector Chris@87: # norms are being computed, so also using `ord='fro'` raises a Chris@87: # ValueError. Chris@87: assert_raises(ValueError, norm, A, 'fro', 0) Chris@87: assert_raises(ValueError, norm, [3, 4], 'fro', None) Chris@87: Chris@87: # Similarly, norm should raise an exception when ord is any finite Chris@87: # number other than 1, 2, -1 or -2 when computing matrix norms. Chris@87: for order in [0, 3]: Chris@87: assert_raises(ValueError, norm, A, order, None) Chris@87: assert_raises(ValueError, norm, A, order, (0, 1)) Chris@87: assert_raises(ValueError, norm, B, order, (1, 2)) Chris@87: Chris@87: # Invalid axis Chris@87: assert_raises(ValueError, norm, B, None, 3) Chris@87: assert_raises(ValueError, norm, B, None, (2, 3)) Chris@87: assert_raises(ValueError, norm, B, None, (0, 1, 2)) Chris@87: Chris@87: def test_longdouble_norm(self): Chris@87: # Non-regression test: p-norm of longdouble would previously raise Chris@87: # UnboundLocalError. Chris@87: x = np.arange(10, dtype=np.longdouble) Chris@87: old_assert_almost_equal(norm(x, ord=3), 12.65, decimal=2) Chris@87: Chris@87: def test_intmin(self): Chris@87: # Non-regression test: p-norm of signed integer would previously do Chris@87: # float cast and abs in the wrong order. Chris@87: x = np.array([-2 ** 31], dtype=np.int32) Chris@87: old_assert_almost_equal(norm(x, ord=3), 2 ** 31, decimal=5) Chris@87: Chris@87: def test_complex_high_ord(self): Chris@87: # gh-4156 Chris@87: d = np.empty((2,), dtype=np.clongdouble) Chris@87: d[0] = 6+7j Chris@87: d[1] = -6+7j Chris@87: res = 11.615898132184 Chris@87: old_assert_almost_equal(np.linalg.norm(d, ord=3), res, decimal=10) Chris@87: d = d.astype(np.complex128) Chris@87: old_assert_almost_equal(np.linalg.norm(d, ord=3), res, decimal=9) Chris@87: d = d.astype(np.complex64) Chris@87: old_assert_almost_equal(np.linalg.norm(d, ord=3), res, decimal=5) Chris@87: Chris@87: Chris@87: class TestNormDouble(_TestNorm): Chris@87: dt = np.double Chris@87: dec = 12 Chris@87: Chris@87: Chris@87: class TestNormSingle(_TestNorm): Chris@87: dt = np.float32 Chris@87: dec = 6 Chris@87: Chris@87: Chris@87: class TestNormInt64(_TestNorm): Chris@87: dt = np.int64 Chris@87: dec = 12 Chris@87: Chris@87: Chris@87: class TestMatrixRank(object): Chris@87: def test_matrix_rank(self): Chris@87: # Full rank matrix Chris@87: yield assert_equal, 4, matrix_rank(np.eye(4)) Chris@87: # rank deficient matrix Chris@87: I=np.eye(4); I[-1, -1] = 0. Chris@87: yield assert_equal, matrix_rank(I), 3 Chris@87: # All zeros - zero rank Chris@87: yield assert_equal, matrix_rank(np.zeros((4, 4))), 0 Chris@87: # 1 dimension - rank 1 unless all 0 Chris@87: yield assert_equal, matrix_rank([1, 0, 0, 0]), 1 Chris@87: yield assert_equal, matrix_rank(np.zeros((4,))), 0 Chris@87: # accepts array-like Chris@87: yield assert_equal, matrix_rank([1]), 1 Chris@87: # greater than 2 dimensions raises error Chris@87: yield assert_raises, TypeError, matrix_rank, np.zeros((2, 2, 2)) Chris@87: # works on scalar Chris@87: yield assert_equal, matrix_rank(1), 1 Chris@87: Chris@87: Chris@87: def test_reduced_rank(): Chris@87: # Test matrices with reduced rank Chris@87: rng = np.random.RandomState(20120714) Chris@87: for i in range(100): Chris@87: # Make a rank deficient matrix Chris@87: X = rng.normal(size=(40, 10)) Chris@87: X[:, 0] = X[:, 1] + X[:, 2] Chris@87: # Assert that matrix_rank detected deficiency Chris@87: assert_equal(matrix_rank(X), 9) Chris@87: X[:, 3] = X[:, 4] + X[:, 5] Chris@87: assert_equal(matrix_rank(X), 8) Chris@87: Chris@87: Chris@87: class TestQR(object): Chris@87: Chris@87: def check_qr(self, a): Chris@87: # This test expects the argument `a` to be an ndarray or Chris@87: # a subclass of an ndarray of inexact type. Chris@87: a_type = type(a) Chris@87: a_dtype = a.dtype Chris@87: m, n = a.shape Chris@87: k = min(m, n) Chris@87: Chris@87: # mode == 'complete' Chris@87: q, r = linalg.qr(a, mode='complete') Chris@87: assert_(q.dtype == a_dtype) Chris@87: assert_(r.dtype == a_dtype) Chris@87: assert_(isinstance(q, a_type)) Chris@87: assert_(isinstance(r, a_type)) Chris@87: assert_(q.shape == (m, m)) Chris@87: assert_(r.shape == (m, n)) Chris@87: assert_almost_equal(dot(q, r), a) Chris@87: assert_almost_equal(dot(q.T.conj(), q), np.eye(m)) Chris@87: assert_almost_equal(np.triu(r), r) Chris@87: Chris@87: # mode == 'reduced' Chris@87: q1, r1 = linalg.qr(a, mode='reduced') Chris@87: assert_(q1.dtype == a_dtype) Chris@87: assert_(r1.dtype == a_dtype) Chris@87: assert_(isinstance(q1, a_type)) Chris@87: assert_(isinstance(r1, a_type)) Chris@87: assert_(q1.shape == (m, k)) Chris@87: assert_(r1.shape == (k, n)) Chris@87: assert_almost_equal(dot(q1, r1), a) Chris@87: assert_almost_equal(dot(q1.T.conj(), q1), np.eye(k)) Chris@87: assert_almost_equal(np.triu(r1), r1) Chris@87: Chris@87: # mode == 'r' Chris@87: r2 = linalg.qr(a, mode='r') Chris@87: assert_(r2.dtype == a_dtype) Chris@87: assert_(isinstance(r2, a_type)) Chris@87: assert_almost_equal(r2, r1) Chris@87: Chris@87: def test_qr_empty(self): Chris@87: a = np.zeros((0, 2)) Chris@87: assert_raises(linalg.LinAlgError, linalg.qr, a) Chris@87: Chris@87: def test_mode_raw(self): Chris@87: # The factorization is not unique and varies between libraries, Chris@87: # so it is not possible to check against known values. Functional Chris@87: # testing is a possibility, but awaits the exposure of more Chris@87: # of the functions in lapack_lite. Consequently, this test is Chris@87: # very limited in scope. Note that the results are in FORTRAN Chris@87: # order, hence the h arrays are transposed. Chris@87: a = array([[1, 2], [3, 4], [5, 6]], dtype=np.double) Chris@87: b = a.astype(np.single) Chris@87: Chris@87: # Test double Chris@87: h, tau = linalg.qr(a, mode='raw') Chris@87: assert_(h.dtype == np.double) Chris@87: assert_(tau.dtype == np.double) Chris@87: assert_(h.shape == (2, 3)) Chris@87: assert_(tau.shape == (2,)) Chris@87: Chris@87: h, tau = linalg.qr(a.T, mode='raw') Chris@87: assert_(h.dtype == np.double) Chris@87: assert_(tau.dtype == np.double) Chris@87: assert_(h.shape == (3, 2)) Chris@87: assert_(tau.shape == (2,)) Chris@87: Chris@87: def test_mode_all_but_economic(self): Chris@87: a = array([[1, 2], [3, 4]]) Chris@87: b = array([[1, 2], [3, 4], [5, 6]]) Chris@87: for dt in "fd": Chris@87: m1 = a.astype(dt) Chris@87: m2 = b.astype(dt) Chris@87: self.check_qr(m1) Chris@87: self.check_qr(m2) Chris@87: self.check_qr(m2.T) Chris@87: self.check_qr(matrix(m1)) Chris@87: for dt in "fd": Chris@87: m1 = 1 + 1j * a.astype(dt) Chris@87: m2 = 1 + 1j * b.astype(dt) Chris@87: self.check_qr(m1) Chris@87: self.check_qr(m2) Chris@87: self.check_qr(m2.T) Chris@87: self.check_qr(matrix(m1)) Chris@87: Chris@87: Chris@87: def test_byteorder_check(): Chris@87: # Byte order check should pass for native order Chris@87: if sys.byteorder == 'little': Chris@87: native = '<' Chris@87: else: Chris@87: native = '>' Chris@87: Chris@87: for dtt in (np.float32, np.float64): Chris@87: arr = np.eye(4, dtype=dtt) Chris@87: n_arr = arr.newbyteorder(native) Chris@87: sw_arr = arr.newbyteorder('S').byteswap() Chris@87: assert_equal(arr.dtype.byteorder, '=') Chris@87: for routine in (linalg.inv, linalg.det, linalg.pinv): Chris@87: # Normal call Chris@87: res = routine(arr) Chris@87: # Native but not '=' Chris@87: assert_array_equal(res, routine(n_arr)) Chris@87: # Swapped Chris@87: assert_array_equal(res, routine(sw_arr)) Chris@87: Chris@87: Chris@87: def test_generalized_raise_multiloop(): Chris@87: # It should raise an error even if the error doesn't occur in the Chris@87: # last iteration of the ufunc inner loop Chris@87: Chris@87: invertible = np.array([[1, 2], [3, 4]]) Chris@87: non_invertible = np.array([[1, 1], [1, 1]]) Chris@87: Chris@87: x = np.zeros([4, 4, 2, 2])[1::2] Chris@87: x[...] = invertible Chris@87: x[0, 0] = non_invertible Chris@87: Chris@87: assert_raises(np.linalg.LinAlgError, np.linalg.inv, x) Chris@87: Chris@87: def test_xerbla_override(): Chris@87: # Check that our xerbla has been successfully linked in. If it is not, Chris@87: # the default xerbla routine is called, which prints a message to stdout Chris@87: # and may, or may not, abort the process depending on the LAPACK package. Chris@87: from nose import SkipTest Chris@87: Chris@87: try: Chris@87: pid = os.fork() Chris@87: except (OSError, AttributeError): Chris@87: # fork failed, or not running on POSIX Chris@87: raise SkipTest("Not POSIX or fork failed.") Chris@87: Chris@87: if pid == 0: Chris@87: # child; close i/o file handles Chris@87: os.close(1) Chris@87: os.close(0) Chris@87: # Avoid producing core files. Chris@87: import resource Chris@87: resource.setrlimit(resource.RLIMIT_CORE, (0, 0)) Chris@87: # These calls may abort. Chris@87: try: Chris@87: np.linalg.lapack_lite.xerbla() Chris@87: except ValueError: Chris@87: pass Chris@87: except: Chris@87: os._exit(os.EX_CONFIG) Chris@87: Chris@87: try: Chris@87: a = np.array([[1.]]) Chris@87: np.linalg.lapack_lite.dorgqr( Chris@87: 1, 1, 1, a, Chris@87: 0, # <- invalid value Chris@87: a, a, 0, 0) Chris@87: except ValueError as e: Chris@87: if "DORGQR parameter number 5" in str(e): Chris@87: # success Chris@87: os._exit(os.EX_OK) Chris@87: Chris@87: # Did not abort, but our xerbla was not linked in. Chris@87: os._exit(os.EX_CONFIG) Chris@87: else: Chris@87: # parent Chris@87: pid, status = os.wait() Chris@87: if os.WEXITSTATUS(status) != os.EX_OK or os.WIFSIGNALED(status): Chris@87: raise SkipTest('Numpy xerbla not linked in.') Chris@87: Chris@87: Chris@87: if __name__ == "__main__": Chris@87: run_module_suite()