import contextlib import gc from itertools import product, cycle import sys import warnings from numbers import Number, Integral import platform import numpy as np from numba import jit, njit, typeof from numba.core import errors from numba.tests.support import (TestCase, tag, needs_lapack, needs_blas, _is_armv7l) from .matmul_usecase import matmul_usecase import unittest def dot2(a, b): return np.dot(a, b) def dot3(a, b, out): return np.dot(a, b, out=out) def vdot(a, b): return np.vdot(a, b) class TestProduct(TestCase): """ Tests for dot products. """ dtypes = (np.float64, np.float32, np.complex128, np.complex64) def setUp(self): # Collect leftovers from previous test cases before checking for leaks gc.collect() def sample_vector(self, n, dtype): # Be careful to generate only exactly representable float values, # to avoid rounding discrepancies between Numpy and Numba base = np.arange(n) if issubclass(dtype, np.complexfloating): return (base * (1 - 0.5j) + 2j).astype(dtype) else: return (base * 0.5 - 1).astype(dtype) def sample_matrix(self, m, n, dtype): return self.sample_vector(m * n, dtype).reshape((m, n)) @contextlib.contextmanager def check_contiguity_warning(self, pyfunc): """ Check performance warning(s) for non-contiguity. """ with warnings.catch_warnings(record=True) as w: warnings.simplefilter('always', errors.NumbaPerformanceWarning) yield self.assertGreaterEqual(len(w), 1) self.assertIs(w[0].category, errors.NumbaPerformanceWarning) self.assertIn("faster on contiguous arrays", str(w[0].message)) self.assertEqual(w[0].filename, pyfunc.__code__.co_filename) # This works because our functions are one-liners self.assertEqual(w[0].lineno, pyfunc.__code__.co_firstlineno + 1) def check_func(self, pyfunc, cfunc, args): with self.assertNoNRTLeak(): expected = pyfunc(*args) got = cfunc(*args) self.assertPreciseEqual(got, expected, ignore_sign_on_zero=True) del got, expected def _aligned_copy(self, arr): # This exists for armv7l because NumPy wants aligned arrays for the # `out` arg of functions, but np.empty/np.copy doesn't seem to always # produce them, in particular for complex dtypes size = (arr.size + 1) * arr.itemsize + 1 datasize = arr.size * arr.itemsize tmp = np.empty(size, dtype=np.uint8) for i in range(arr.itemsize + 1): new = tmp[i : i + datasize].view(dtype=arr.dtype) if new.flags.aligned: break else: raise Exception("Could not obtain aligned array") if arr.flags.c_contiguous: new = np.reshape(new, arr.shape, order='C') else: new = np.reshape(new, arr.shape, order='F') new[:] = arr[:] assert new.flags.aligned return new def check_func_out(self, pyfunc, cfunc, args, out): copier = self._aligned_copy if _is_armv7l else np.copy with self.assertNoNRTLeak(): expected = copier(out) got = copier(out) self.assertIs(pyfunc(*args, out=expected), expected) self.assertIs(cfunc(*args, out=got), got) self.assertPreciseEqual(got, expected, ignore_sign_on_zero=True) del got, expected def assert_mismatching_sizes(self, cfunc, args, is_out=False): with self.assertRaises(ValueError) as raises: cfunc(*args) msg = ("incompatible output array size" if is_out else "incompatible array sizes") self.assertIn(msg, str(raises.exception)) def assert_mismatching_dtypes(self, cfunc, args, func_name="np.dot()"): with self.assertRaises(errors.TypingError) as raises: cfunc(*args) self.assertIn("%s arguments must all have the same dtype" % (func_name,), str(raises.exception)) def check_dot_vv(self, pyfunc, func_name): n = 3 cfunc = jit(nopython=True)(pyfunc) for dtype in self.dtypes: a = self.sample_vector(n, dtype) b = self.sample_vector(n, dtype) self.check_func(pyfunc, cfunc, (a, b)) # Non-contiguous self.check_func(pyfunc, cfunc, (a[::-1], b[::-1])) # Mismatching sizes a = self.sample_vector(n - 1, np.float64) b = self.sample_vector(n, np.float64) self.assert_mismatching_sizes(cfunc, (a, b)) # Mismatching dtypes a = self.sample_vector(n, np.float32) b = self.sample_vector(n, np.float64) self.assert_mismatching_dtypes(cfunc, (a, b), func_name=func_name) @needs_blas def test_dot_vv(self): """ Test vector * vector np.dot() """ self.check_dot_vv(dot2, "np.dot()") @needs_blas def test_vdot(self): """ Test np.vdot() """ self.check_dot_vv(vdot, "np.vdot()") def check_dot_vm(self, pyfunc2, pyfunc3, func_name): def samples(m, n): for order in 'CF': a = self.sample_matrix(m, n, np.float64).copy(order=order) b = self.sample_vector(n, np.float64) yield a, b for dtype in self.dtypes: a = self.sample_matrix(m, n, dtype) b = self.sample_vector(n, dtype) yield a, b # Non-contiguous yield a[::-1], b[::-1] cfunc2 = jit(nopython=True)(pyfunc2) if pyfunc3 is not None: cfunc3 = jit(nopython=True)(pyfunc3) for m, n in [(2, 3), (3, 0), (0, 3) ]: for a, b in samples(m, n): self.check_func(pyfunc2, cfunc2, (a, b)) self.check_func(pyfunc2, cfunc2, (b, a.T)) if pyfunc3 is not None: for a, b in samples(m, n): out = np.empty(m, dtype=a.dtype) self.check_func_out(pyfunc3, cfunc3, (a, b), out) self.check_func_out(pyfunc3, cfunc3, (b, a.T), out) # Mismatching sizes m, n = 2, 3 a = self.sample_matrix(m, n - 1, np.float64) b = self.sample_vector(n, np.float64) self.assert_mismatching_sizes(cfunc2, (a, b)) self.assert_mismatching_sizes(cfunc2, (b, a.T)) if pyfunc3 is not None: out = np.empty(m, np.float64) self.assert_mismatching_sizes(cfunc3, (a, b, out)) self.assert_mismatching_sizes(cfunc3, (b, a.T, out)) a = self.sample_matrix(m, m, np.float64) b = self.sample_vector(m, np.float64) out = np.empty(m - 1, np.float64) self.assert_mismatching_sizes(cfunc3, (a, b, out), is_out=True) self.assert_mismatching_sizes(cfunc3, (b, a.T, out), is_out=True) # Mismatching dtypes a = self.sample_matrix(m, n, np.float32) b = self.sample_vector(n, np.float64) self.assert_mismatching_dtypes(cfunc2, (a, b), func_name) if pyfunc3 is not None: a = self.sample_matrix(m, n, np.float64) b = self.sample_vector(n, np.float64) out = np.empty(m, np.float32) self.assert_mismatching_dtypes(cfunc3, (a, b, out), func_name) @needs_blas def test_dot_vm(self): """ Test vector * matrix and matrix * vector np.dot() """ self.check_dot_vm(dot2, dot3, "np.dot()") def check_dot_mm(self, pyfunc2, pyfunc3, func_name): def samples(m, n, k): for order_a, order_b in product('CF', 'CF'): a = self.sample_matrix(m, k, np.float64).copy(order=order_a) b = self.sample_matrix(k, n, np.float64).copy(order=order_b) yield a, b for dtype in self.dtypes: a = self.sample_matrix(m, k, dtype) b = self.sample_matrix(k, n, dtype) yield a, b # Non-contiguous yield a[::-1], b[::-1] cfunc2 = jit(nopython=True)(pyfunc2) if pyfunc3 is not None: cfunc3 = jit(nopython=True)(pyfunc3) # Test generic matrix * matrix as well as "degenerate" cases # where one of the outer dimensions is 1 (i.e. really represents # a vector, which may select a different implementation), # one of the matrices is empty, or both matrices are empty. for m, n, k in [(2, 3, 4), # Generic matrix * matrix (1, 3, 4), # 2d vector * matrix (1, 1, 4), # 2d vector * 2d vector (0, 3, 2), # Empty matrix * matrix, empty output (3, 0, 2), # Matrix * empty matrix, empty output (0, 0, 3), # Both arguments empty, empty output (3, 2, 0), # Both arguments empty, nonempty output ]: for a, b in samples(m, n, k): self.check_func(pyfunc2, cfunc2, (a, b)) self.check_func(pyfunc2, cfunc2, (b.T, a.T)) if pyfunc3 is not None: for a, b in samples(m, n, k): out = np.empty((m, n), dtype=a.dtype) self.check_func_out(pyfunc3, cfunc3, (a, b), out) out = np.empty((n, m), dtype=a.dtype) self.check_func_out(pyfunc3, cfunc3, (b.T, a.T), out) # Mismatching sizes m, n, k = 2, 3, 4 a = self.sample_matrix(m, k - 1, np.float64) b = self.sample_matrix(k, n, np.float64) self.assert_mismatching_sizes(cfunc2, (a, b)) if pyfunc3 is not None: out = np.empty((m, n), np.float64) self.assert_mismatching_sizes(cfunc3, (a, b, out)) a = self.sample_matrix(m, k, np.float64) b = self.sample_matrix(k, n, np.float64) out = np.empty((m, n - 1), np.float64) self.assert_mismatching_sizes(cfunc3, (a, b, out), is_out=True) # Mismatching dtypes a = self.sample_matrix(m, k, np.float32) b = self.sample_matrix(k, n, np.float64) self.assert_mismatching_dtypes(cfunc2, (a, b), func_name) if pyfunc3 is not None: a = self.sample_matrix(m, k, np.float64) b = self.sample_matrix(k, n, np.float64) out = np.empty((m, n), np.float32) self.assert_mismatching_dtypes(cfunc3, (a, b, out), func_name) @needs_blas def test_dot_mm(self): """ Test matrix * matrix np.dot() """ self.check_dot_mm(dot2, dot3, "np.dot()") @needs_blas def test_matmul_vv(self): """ Test vector @ vector """ self.check_dot_vv(matmul_usecase, "'@'") @needs_blas def test_matmul_vm(self): """ Test vector @ matrix and matrix @ vector """ self.check_dot_vm(matmul_usecase, None, "'@'") @needs_blas def test_matmul_mm(self): """ Test matrix @ matrix """ self.check_dot_mm(matmul_usecase, None, "'@'") @needs_blas def test_contiguity_warnings(self): m, k, n = 2, 3, 4 dtype = np.float64 a = self.sample_matrix(m, k, dtype)[::-1] b = self.sample_matrix(k, n, dtype)[::-1] out = np.empty((m, n), dtype) cfunc = jit(nopython=True)(dot2) with self.check_contiguity_warning(cfunc.py_func): cfunc(a, b) cfunc = jit(nopython=True)(dot3) with self.check_contiguity_warning(cfunc.py_func): cfunc(a, b, out) a = self.sample_vector(n, dtype)[::-1] b = self.sample_vector(n, dtype)[::-1] cfunc = jit(nopython=True)(vdot) with self.check_contiguity_warning(cfunc.py_func): cfunc(a, b) # Implementation definitions for the purpose of jitting. def invert_matrix(a): return np.linalg.inv(a) def cholesky_matrix(a): return np.linalg.cholesky(a) def eig_matrix(a): return np.linalg.eig(a) def eigvals_matrix(a): return np.linalg.eigvals(a) def eigh_matrix(a): return np.linalg.eigh(a) def eigvalsh_matrix(a): return np.linalg.eigvalsh(a) def svd_matrix(a, full_matrices=1): return np.linalg.svd(a, full_matrices) def qr_matrix(a): return np.linalg.qr(a) def lstsq_system(A, B, rcond=-1): return np.linalg.lstsq(A, B, rcond) def solve_system(A, B): return np.linalg.solve(A, B) def pinv_matrix(A, rcond=1e-15): # 1e-15 from numpy impl return np.linalg.pinv(A) def slogdet_matrix(a): return np.linalg.slogdet(a) def det_matrix(a): return np.linalg.det(a) def norm_matrix(a, ord=None): return np.linalg.norm(a, ord) def cond_matrix(a, p=None): return np.linalg.cond(a, p) def matrix_rank_matrix(a, tol=None): return np.linalg.matrix_rank(a, tol) def matrix_power_matrix(a, n): return np.linalg.matrix_power(a, n) def trace_matrix(a, offset=0): return np.trace(a, offset) def trace_matrix_no_offset(a): return np.trace(a) def outer_matrix(a, b, out=None): return np.outer(a, b, out=out) def kron_matrix(a, b): return np.kron(a, b) class TestLinalgBase(TestCase): """ Provides setUp and common data/error modes for testing np.linalg functions. """ # supported dtypes dtypes = (np.float64, np.float32, np.complex128, np.complex64) def setUp(self): # Collect leftovers from previous test cases before checking for leaks gc.collect() def sample_vector(self, n, dtype): # Be careful to generate only exactly representable float values, # to avoid rounding discrepancies between Numpy and Numba base = np.arange(n) if issubclass(dtype, np.complexfloating): return (base * (1 - 0.5j) + 2j).astype(dtype) else: return (base * 0.5 + 1).astype(dtype) def specific_sample_matrix( self, size, dtype, order, rank=None, condition=None): """ Provides a sample matrix with an optionally specified rank or condition number. size: (rows, columns), the dimensions of the returned matrix. dtype: the dtype for the returned matrix. order: the memory layout for the returned matrix, 'F' or 'C'. rank: the rank of the matrix, an integer value, defaults to full rank. condition: the condition number of the matrix (defaults to 1.) NOTE: Only one of rank or condition may be set. """ # default condition d_cond = 1. if len(size) != 2: raise ValueError("size must be a length 2 tuple.") if order not in ['F', 'C']: raise ValueError("order must be one of 'F' or 'C'.") if dtype not in [np.float32, np.float64, np.complex64, np.complex128]: raise ValueError("dtype must be a numpy floating point type.") if rank is not None and condition is not None: raise ValueError("Only one of rank or condition can be specified.") if condition is None: condition = d_cond if condition < 1: raise ValueError("Condition number must be >=1.") np.random.seed(0) # repeatable seed m, n = size if m < 0 or n < 0: raise ValueError("Negative dimensions given for matrix shape.") minmn = min(m, n) if rank is None: rv = minmn else: if rank <= 0: raise ValueError("Rank must be greater than zero.") if not isinstance(rank, Integral): raise ValueError("Rank must an integer.") rv = rank if rank > minmn: raise ValueError("Rank given greater than full rank.") if m == 1 or n == 1: # vector, must be rank 1 (enforced above) # condition of vector is also 1 if condition != d_cond: raise ValueError( "Condition number was specified for a vector (always 1.).") maxmn = max(m, n) Q = self.sample_vector(maxmn, dtype).reshape(m, n) else: # Build a sample matrix via combining SVD like inputs. # Create matrices of left and right singular vectors. # This could use Modified Gram-Schmidt and perhaps be quicker, # at present it uses QR decompositions to obtain orthonormal # matrices. tmp = self.sample_vector(m * m, dtype).reshape(m, m) U, _ = np.linalg.qr(tmp) # flip the second array, else for m==n the identity matrix appears tmp = self.sample_vector(n * n, dtype)[::-1].reshape(n, n) V, _ = np.linalg.qr(tmp) # create singular values. sv = np.linspace(d_cond, condition, rv) S = np.zeros((m, n)) idx = np.nonzero(np.eye(m, n)) S[idx[0][:rv], idx[1][:rv]] = sv Q = np.dot(np.dot(U, S), V.T) # construct Q = np.array(Q, dtype=dtype, order=order) # sort out order/type return Q def assert_error(self, cfunc, args, msg, err=ValueError): with self.assertRaises(err) as raises: cfunc(*args) self.assertIn(msg, str(raises.exception)) def assert_non_square(self, cfunc, args): msg = "Last 2 dimensions of the array must be square." self.assert_error(cfunc, args, msg, np.linalg.LinAlgError) def assert_wrong_dtype(self, name, cfunc, args): msg = "np.linalg.%s() only supported on float and complex arrays" % name self.assert_error(cfunc, args, msg, errors.TypingError) def assert_wrong_dimensions(self, name, cfunc, args, la_prefix=True): prefix = "np.linalg" if la_prefix else "np" msg = "%s.%s() only supported on 2-D arrays" % (prefix, name) self.assert_error(cfunc, args, msg, errors.TypingError) def assert_no_nan_or_inf(self, cfunc, args): msg = "Array must not contain infs or NaNs." self.assert_error(cfunc, args, msg, np.linalg.LinAlgError) def assert_contig_sanity(self, got, expected_contig): """ This checks that in a computed result from numba (array, possibly tuple of arrays) all the arrays are contiguous in memory and that they are all at least one of "C_CONTIGUOUS" or "F_CONTIGUOUS". The computed result of the contiguousness is then compared against a hardcoded expected result. got: is the computed results from numba expected_contig: is "C" or "F" and is the expected type of contiguousness across all input values (and therefore tests). """ if isinstance(got, tuple): # tuple present, check all results for a in got: self.assert_contig_sanity(a, expected_contig) else: if not isinstance(got, Number): # else a single array is present c_contig = got.flags.c_contiguous f_contig = got.flags.f_contiguous # check that the result (possible set of) is at least one of # C or F contiguous. msg = "Results are not at least one of all C or F contiguous." self.assertTrue(c_contig | f_contig, msg) msg = "Computed contiguousness does not match expected." if expected_contig == "C": self.assertTrue(c_contig, msg) elif expected_contig == "F": self.assertTrue(f_contig, msg) else: raise ValueError("Unknown contig") def assert_raise_on_singular(self, cfunc, args): msg = "Matrix is singular to machine precision." self.assert_error(cfunc, args, msg, err=np.linalg.LinAlgError) def assert_is_identity_matrix(self, got, rtol=None, atol=None): """ Checks if a matrix is equal to the identity matrix. """ # check it is square self.assertEqual(got.shape[-1], got.shape[-2]) # create identity matrix eye = np.eye(got.shape[-1], dtype=got.dtype) resolution = 5 * np.finfo(got.dtype).resolution if rtol is None: rtol = 10 * resolution if atol is None: atol = 100 * resolution # zeros tend to be fuzzy # check it matches np.testing.assert_allclose(got, eye, rtol, atol) def assert_invalid_norm_kind(self, cfunc, args): """ For use in norm() and cond() tests. """ msg = "Invalid norm order for matrices." self.assert_error(cfunc, args, msg, ValueError) def assert_raise_on_empty(self, cfunc, args): msg = 'Arrays cannot be empty' self.assert_error(cfunc, args, msg, np.linalg.LinAlgError) class TestTestLinalgBase(TestCase): """ The sample matrix code TestLinalgBase.specific_sample_matrix() is a bit involved, this class tests it works as intended. """ def test_specific_sample_matrix(self): # add a default test to the ctor, it never runs so doesn't matter inst = TestLinalgBase('specific_sample_matrix') sizes = [(7, 1), (11, 5), (5, 11), (3, 3), (1, 7)] # test loop for size, dtype, order in product(sizes, inst.dtypes, 'FC'): m, n = size minmn = min(m, n) # test default full rank A = inst.specific_sample_matrix(size, dtype, order) self.assertEqual(A.shape, size) self.assertEqual(np.linalg.matrix_rank(A), minmn) # test reduced rank if a reduction is possible if minmn > 1: rank = minmn - 1 A = inst.specific_sample_matrix(size, dtype, order, rank=rank) self.assertEqual(A.shape, size) self.assertEqual(np.linalg.matrix_rank(A), rank) resolution = 5 * np.finfo(dtype).resolution # test default condition A = inst.specific_sample_matrix(size, dtype, order) self.assertEqual(A.shape, size) np.testing.assert_allclose(np.linalg.cond(A), 1., rtol=resolution, atol=resolution) # test specified condition if matrix is > 1D if minmn > 1: condition = 10. A = inst.specific_sample_matrix( size, dtype, order, condition=condition) self.assertEqual(A.shape, size) np.testing.assert_allclose(np.linalg.cond(A), 10., rtol=resolution, atol=resolution) # check errors are raised appropriately def check_error(args, msg, err=ValueError): with self.assertRaises(err) as raises: inst.specific_sample_matrix(*args) self.assertIn(msg, str(raises.exception)) # check the checker runs ok with self.assertRaises(AssertionError) as raises: msg = "blank" check_error(((2, 3), np.float64, 'F'), msg, err=ValueError) # check invalid inputs... # bad size msg = "size must be a length 2 tuple." check_error(((1,), np.float64, 'F'), msg, err=ValueError) # bad order msg = "order must be one of 'F' or 'C'." check_error(((2, 3), np.float64, 'z'), msg, err=ValueError) # bad type msg = "dtype must be a numpy floating point type." check_error(((2, 3), np.int32, 'F'), msg, err=ValueError) # specifying both rank and condition msg = "Only one of rank or condition can be specified." check_error(((2, 3), np.float64, 'F', 1, 1), msg, err=ValueError) # specifying negative condition msg = "Condition number must be >=1." check_error(((2, 3), np.float64, 'F', None, -1), msg, err=ValueError) # specifying negative matrix dimension msg = "Negative dimensions given for matrix shape." check_error(((2, -3), np.float64, 'F'), msg, err=ValueError) # specifying negative rank msg = "Rank must be greater than zero." check_error(((2, 3), np.float64, 'F', -1), msg, err=ValueError) # specifying a rank greater than maximum rank msg = "Rank given greater than full rank." check_error(((2, 3), np.float64, 'F', 4), msg, err=ValueError) # specifying a condition number for a vector msg = "Condition number was specified for a vector (always 1.)." check_error(((1, 3), np.float64, 'F', None, 10), msg, err=ValueError) # specifying a non integer rank msg = "Rank must an integer." check_error(((2, 3), np.float64, 'F', 1.5), msg, err=ValueError) class TestLinalgInv(TestLinalgBase): """ Tests for np.linalg.inv. """ @needs_lapack def test_linalg_inv(self): """ Test np.linalg.inv """ n = 10 cfunc = jit(nopython=True)(invert_matrix) def check(a, **kwargs): expected = invert_matrix(a) got = cfunc(a) self.assert_contig_sanity(got, "F") use_reconstruction = False # try strict try: np.testing.assert_array_almost_equal_nulp(got, expected, nulp=10) except AssertionError: # fall back to reconstruction use_reconstruction = True if use_reconstruction: rec = np.dot(got, a) self.assert_is_identity_matrix(rec) # Ensure proper resource management with self.assertNoNRTLeak(): cfunc(a) for dtype, order in product(self.dtypes, 'CF'): a = self.specific_sample_matrix((n, n), dtype, order) check(a) # 0 dimensioned matrix check(np.empty((0, 0))) # Non square matrix self.assert_non_square(cfunc, (np.ones((2, 3)),)) # Wrong dtype self.assert_wrong_dtype("inv", cfunc, (np.ones((2, 2), dtype=np.int32),)) # Dimension issue self.assert_wrong_dimensions("inv", cfunc, (np.ones(10),)) # Singular matrix self.assert_raise_on_singular(cfunc, (np.zeros((2, 2)),)) @needs_lapack def test_no_input_mutation(self): X = np.array([[1., 3, 2, 7,], [-5, 4, 2, 3,], [9, -3, 1, 1,], [2, -2, 2, 8,]], order='F') X_orig = np.copy(X) @jit(nopython=True) def ainv(X, test): if test: # not executed, but necessary to trigger A ordering in X X = X[1:2, :] return np.linalg.inv(X) expected = ainv.py_func(X, False) np.testing.assert_allclose(X, X_orig) got = ainv(X, False) np.testing.assert_allclose(X, X_orig) np.testing.assert_allclose(expected, got) class TestLinalgCholesky(TestLinalgBase): """ Tests for np.linalg.cholesky. """ def sample_matrix(self, m, dtype, order): # pd. (positive definite) matrix has eigenvalues in Z+ np.random.seed(0) # repeatable seed A = np.random.rand(m, m) # orthonormal q needed to form up q^{-1}*D*q # no "orth()" in numpy q, _ = np.linalg.qr(A) L = np.arange(1, m + 1) # some positive eigenvalues Q = np.dot(np.dot(q.T, np.diag(L)), q) # construct Q = np.array(Q, dtype=dtype, order=order) # sort out order/type return Q def assert_not_pd(self, cfunc, args): msg = "Matrix is not positive definite." self.assert_error(cfunc, args, msg, np.linalg.LinAlgError) @needs_lapack def test_linalg_cholesky(self): """ Test np.linalg.cholesky """ n = 10 cfunc = jit(nopython=True)(cholesky_matrix) def check(a): expected = cholesky_matrix(a) got = cfunc(a) use_reconstruction = False # check that the computed results are contig and in the same way self.assert_contig_sanity(got, "C") # try strict try: np.testing.assert_array_almost_equal_nulp(got, expected, nulp=10) except AssertionError: # fall back to reconstruction use_reconstruction = True # try via reconstruction if use_reconstruction: rec = np.dot(got, np.conj(got.T)) resolution = 5 * np.finfo(a.dtype).resolution np.testing.assert_allclose( a, rec, rtol=resolution, atol=resolution ) # Ensure proper resource management with self.assertNoNRTLeak(): cfunc(a) for dtype, order in product(self.dtypes, 'FC'): a = self.sample_matrix(n, dtype, order) check(a) # 0 dimensioned matrix check(np.empty((0, 0))) rn = "cholesky" # Non square matrices self.assert_non_square(cfunc, (np.ones((2, 3), dtype=np.float64),)) # Wrong dtype self.assert_wrong_dtype(rn, cfunc, (np.ones((2, 2), dtype=np.int32),)) # Dimension issue self.assert_wrong_dimensions(rn, cfunc, (np.ones(10, dtype=np.float64),)) # not pd self.assert_not_pd(cfunc, (np.ones(4, dtype=np.float64).reshape(2, 2),)) class TestLinalgEigenSystems(TestLinalgBase): """ Tests for np.linalg.eig/eigvals. """ def sample_matrix(self, m, dtype, order): # This is a tridiag with the same but skewed values on the diagonals v = self.sample_vector(m, dtype) Q = np.diag(v) idx = np.nonzero(np.eye(Q.shape[0], Q.shape[1], 1)) Q[idx] = v[1:] idx = np.nonzero(np.eye(Q.shape[0], Q.shape[1], -1)) Q[idx] = v[:-1] Q = np.array(Q, dtype=dtype, order=order) return Q def assert_no_domain_change(self, name, cfunc, args): msg = name + "() argument must not cause a domain change." self.assert_error(cfunc, args, msg) def _check_worker(self, cfunc, name, expected_res_len, check_for_domain_change): def check(*args): expected = cfunc.py_func(*args) got = cfunc(*args) a = args[0] # check that the returned tuple is same length self.assertEqual(len(expected), len(got)) # and that dimension is correct res_is_tuple = False if isinstance(got, tuple): res_is_tuple = True self.assertEqual(len(got), expected_res_len) else: # its an array self.assertEqual(got.ndim, expected_res_len) # and that the computed results are contig and in the same way self.assert_contig_sanity(got, "F") use_reconstruction = False # try plain match of each array to np first for k in range(len(expected)): try: np.testing.assert_array_almost_equal_nulp( got[k], expected[k], nulp=10) except AssertionError: # plain match failed, test by reconstruction use_reconstruction = True # If plain match fails then reconstruction is used. # this checks that A*V ~== V*diag(W) # i.e. eigensystem ties out # this is required as numpy uses only double precision lapack # routines and computation of eigenvectors is numerically # sensitive, numba uses the type specific routines therefore # sometimes comes out with a different (but entirely # valid) answer (eigenvectors are not unique etc.). # This is only applicable if eigenvectors are computed # along with eigenvalues i.e. result is a tuple. resolution = 5 * np.finfo(a.dtype).resolution if use_reconstruction: if res_is_tuple: w, v = got # modify 'a' if hermitian eigensystem functionality is # being tested. 'L' for use lower part is default and # the only thing used at present so we conjugate transpose # the lower part into the upper for use in the # reconstruction. By construction the sample matrix is # tridiag so this is just a question of copying the lower # diagonal into the upper and conjugating on the way. if name[-1] == 'h': idxl = np.nonzero(np.eye(a.shape[0], a.shape[1], -1)) idxu = np.nonzero(np.eye(a.shape[0], a.shape[1], 1)) cfunc(*args) # upper idx must match lower for default uplo="L" # if complex, conjugate a[idxu] = np.conj(a[idxl]) # also, only the real part of the diagonals is # considered in the calculation so the imag is zeroed # out for the purposes of use in reconstruction. a[np.diag_indices(a.shape[0])] = np.real(np.diag(a)) lhs = np.dot(a, v) rhs = np.dot(v, np.diag(w)) np.testing.assert_allclose( lhs.real, rhs.real, rtol=resolution, atol=resolution ) if np.iscomplexobj(v): np.testing.assert_allclose( lhs.imag, rhs.imag, rtol=resolution, atol=resolution ) else: # This isn't technically reconstruction but is here to # deal with that the order of the returned eigenvalues # may differ in the case of routines just returning # eigenvalues and there's no true reconstruction # available with which to perform a check. np.testing.assert_allclose( np.sort(expected), np.sort(got), rtol=resolution, atol=resolution ) # Ensure proper resource management with self.assertNoNRTLeak(): cfunc(*args) return check def checker_for_linalg_eig( self, name, func, expected_res_len, check_for_domain_change=None): """ Test np.linalg.eig """ n = 10 cfunc = jit(nopython=True)(func) check = self._check_worker(cfunc, name, expected_res_len, check_for_domain_change) # The main test loop for dtype, order in product(self.dtypes, 'FC'): a = self.sample_matrix(n, dtype, order) check(a) # Test both a real and complex type as the impls are different for ty in [np.float32, np.complex64]: # 0 dimensioned matrix check(np.empty((0, 0), dtype=ty)) # Non square matrices self.assert_non_square(cfunc, (np.ones((2, 3), dtype=ty),)) # Wrong dtype self.assert_wrong_dtype(name, cfunc, (np.ones((2, 2), dtype=np.int32),)) # Dimension issue self.assert_wrong_dimensions(name, cfunc, (np.ones(10, dtype=ty),)) # no nans or infs self.assert_no_nan_or_inf(cfunc, (np.array([[1., 2., ], [np.inf, np.nan]], dtype=ty),)) if check_for_domain_change: # By design numba does not support dynamic return types, numpy does # and uses this in the case of returning eigenvalues/vectors of # a real matrix. The return type of np.linalg.eig(), when # operating on a matrix in real space depends on the values present # in the matrix itself (recalling that eigenvalues are the roots of the # characteristic polynomial of the system matrix, which will by # construction depend on the values present in the system matrix). # This test asserts that if a domain change is required on the return # type, i.e. complex eigenvalues from a real input, an error is raised. # For complex types, regardless of the value of the imaginary part of # the returned eigenvalues, a complex type will be returned, this # follows numpy and fits in with numba. # First check that the computation is valid (i.e. in complex space) A = np.array([[1, -2], [2, 1]]) check(A.astype(np.complex128)) # and that the imaginary part is nonzero l, _ = func(A) self.assertTrue(np.any(l.imag)) # Now check that the computation fails in real space for ty in [np.float32, np.float64]: self.assert_no_domain_change(name, cfunc, (A.astype(ty),)) @needs_lapack def test_linalg_eig(self): self.checker_for_linalg_eig("eig", eig_matrix, 2, True) @needs_lapack def test_linalg_eigvals(self): self.checker_for_linalg_eig("eigvals", eigvals_matrix, 1, True) @needs_lapack def test_linalg_eigh(self): self.checker_for_linalg_eig("eigh", eigh_matrix, 2, False) @needs_lapack def test_linalg_eigvalsh(self): self.checker_for_linalg_eig("eigvalsh", eigvalsh_matrix, 1, False) @needs_lapack def test_no_input_mutation(self): # checks inputs are not mutated for c in (('eig', 2, True), ('eigvals', 1, True), ('eigh', 2, False), ('eigvalsh', 1, False)): m, nout, domain_change = c meth = getattr(np.linalg, m) @jit(nopython=True) def func(X, test): if test: # not executed, but necessary to trigger A ordering in X X = X[1:2, :] return meth(X) check = self._check_worker(func, m, nout, domain_change) for dtype in (np.float64, np.complex128): with self.subTest(meth=meth, dtype=dtype): # trivial system, doesn't matter, just checking if it gets # mutated X = np.array([[10., 1, 0, 1], [1, 9, 0, 0], [0, 0, 8, 0], [1, 0, 0, 7], ], order='F', dtype=dtype) X_orig = np.copy(X) expected = func.py_func(X, False) np.testing.assert_allclose(X, X_orig) got = func(X, False) np.testing.assert_allclose(X, X_orig) check(X, False) class TestLinalgSvd(TestLinalgBase): """ Tests for np.linalg.svd. """ @needs_lapack def test_linalg_svd(self): """ Test np.linalg.svd """ cfunc = jit(nopython=True)(svd_matrix) def check(a, **kwargs): expected = svd_matrix(a, **kwargs) got = cfunc(a, **kwargs) # check that the returned tuple is same length self.assertEqual(len(expected), len(got)) # and that length is 3 self.assertEqual(len(got), 3) # and that the computed results are contig and in the same way self.assert_contig_sanity(got, "F") use_reconstruction = False # try plain match of each array to np first for k in range(len(expected)): try: np.testing.assert_array_almost_equal_nulp( got[k], expected[k], nulp=10) except AssertionError: # plain match failed, test by reconstruction use_reconstruction = True # if plain match fails then reconstruction is used. # this checks that A ~= U*S*V**H # i.e. SV decomposition ties out # this is required as numpy uses only double precision lapack # routines and computation of svd is numerically # sensitive, numba using the type specific routines therefore # sometimes comes out with a different answer (orthonormal bases # are not unique etc.). if use_reconstruction: u, sv, vt = got # check they are dimensionally correct for k in range(len(expected)): self.assertEqual(got[k].shape, expected[k].shape) # regardless of full_matrices cols in u and rows in vt # dictates the working size of s s = np.zeros((u.shape[1], vt.shape[0])) np.fill_diagonal(s, sv) rec = np.dot(np.dot(u, s), vt) resolution = np.finfo(a.dtype).resolution np.testing.assert_allclose( a, rec, rtol=10 * resolution, atol=100 * resolution # zeros tend to be fuzzy ) # Ensure proper resource management with self.assertNoNRTLeak(): cfunc(a, **kwargs) # test: column vector, tall, wide, square, row vector # prime sizes sizes = [(7, 1), (7, 5), (5, 7), (3, 3), (1, 7)] # flip on reduced or full matrices full_matrices = (True, False) # test loop for size, dtype, fmat, order in \ product(sizes, self.dtypes, full_matrices, 'FC'): a = self.specific_sample_matrix(size, dtype, order) check(a, full_matrices=fmat) rn = "svd" # Wrong dtype self.assert_wrong_dtype(rn, cfunc, (np.ones((2, 2), dtype=np.int32),)) # Dimension issue self.assert_wrong_dimensions(rn, cfunc, (np.ones(10, dtype=np.float64),)) # no nans or infs self.assert_no_nan_or_inf(cfunc, (np.array([[1., 2., ], [np.inf, np.nan]], dtype=np.float64),)) # empty for sz in [(0, 1), (1, 0), (0, 0)]: args = (np.empty(sz), True) self.assert_raise_on_empty(cfunc, args) @needs_lapack def test_no_input_mutation(self): X = np.array([[1., 3, 2, 7,], [-5, 4, 2, 3,], [9, -3, 1, 1,], [2, -2, 2, 8,]], order='F') X_orig = np.copy(X) @jit(nopython=True) def func(X, test): if test: # not executed, but necessary to trigger A ordering in X X = X[1:2, :] return np.linalg.svd(X) expected = func.py_func(X, False) np.testing.assert_allclose(X, X_orig) got = func(X, False) np.testing.assert_allclose(X, X_orig) for e_a, g_a in zip(expected, got): np.testing.assert_allclose(e_a, g_a) class TestLinalgQr(TestLinalgBase): """ Tests for np.linalg.qr. """ @needs_lapack def test_linalg_qr(self): """ Test np.linalg.qr """ cfunc = jit(nopython=True)(qr_matrix) def check(a, **kwargs): expected = qr_matrix(a, **kwargs) got = cfunc(a, **kwargs) # check that the returned tuple is same length self.assertEqual(len(expected), len(got)) # and that length is 2 self.assertEqual(len(got), 2) # and that the computed results are contig and in the same way self.assert_contig_sanity(got, "F") use_reconstruction = False # try plain match of each array to np first for k in range(len(expected)): try: np.testing.assert_array_almost_equal_nulp( got[k], expected[k], nulp=10) except AssertionError: # plain match failed, test by reconstruction use_reconstruction = True # if plain match fails then reconstruction is used. # this checks that A ~= Q*R and that (Q^H)*Q = I # i.e. QR decomposition ties out # this is required as numpy uses only double precision lapack # routines and computation of qr is numerically # sensitive, numba using the type specific routines therefore # sometimes comes out with a different answer (orthonormal bases # are not unique etc.). if use_reconstruction: q, r = got # check they are dimensionally correct for k in range(len(expected)): self.assertEqual(got[k].shape, expected[k].shape) # check A=q*r rec = np.dot(q, r) resolution = np.finfo(a.dtype).resolution np.testing.assert_allclose( a, rec, rtol=10 * resolution, atol=100 * resolution # zeros tend to be fuzzy ) # check q is orthonormal self.assert_is_identity_matrix(np.dot(np.conjugate(q.T), q)) # Ensure proper resource management with self.assertNoNRTLeak(): cfunc(a, **kwargs) # test: column vector, tall, wide, square, row vector # prime sizes sizes = [(7, 1), (11, 5), (5, 11), (3, 3), (1, 7)] # test loop for size, dtype, order in \ product(sizes, self.dtypes, 'FC'): a = self.specific_sample_matrix(size, dtype, order) check(a) rn = "qr" # Wrong dtype self.assert_wrong_dtype(rn, cfunc, (np.ones((2, 2), dtype=np.int32),)) # Dimension issue self.assert_wrong_dimensions(rn, cfunc, (np.ones(10, dtype=np.float64),)) # no nans or infs self.assert_no_nan_or_inf(cfunc, (np.array([[1., 2., ], [np.inf, np.nan]], dtype=np.float64),)) # empty for sz in [(0, 1), (1, 0), (0, 0)]: self.assert_raise_on_empty(cfunc, (np.empty(sz),)) @needs_lapack def test_no_input_mutation(self): X = np.array([[1., 3, 2, 7,], [-5, 4, 2, 3,], [9, -3, 1, 1,], [2, -2, 2, 8,]], order='F') X_orig = np.copy(X) @jit(nopython=True) def func(X, test): if test: # not executed, but necessary to trigger A ordering in X X = X[1:2, :] return np.linalg.qr(X) expected = func.py_func(X, False) np.testing.assert_allclose(X, X_orig) got = func(X, False) np.testing.assert_allclose(X, X_orig) for e_a, g_a in zip(expected, got): np.testing.assert_allclose(e_a, g_a) class TestLinalgSystems(TestLinalgBase): """ Base class for testing "system" solvers from np.linalg. Namely np.linalg.solve() and np.linalg.lstsq(). """ # check for RHS with dimension > 2 raises def assert_wrong_dimensions_1D(self, name, cfunc, args, la_prefix=True): prefix = "np.linalg" if la_prefix else "np" msg = "%s.%s() only supported on 1 and 2-D arrays" % (prefix, name) self.assert_error(cfunc, args, msg, errors.TypingError) # check that a dimensionally invalid system raises def assert_dimensionally_invalid(self, cfunc, args): msg = "Incompatible array sizes, system is not dimensionally valid." self.assert_error(cfunc, args, msg, np.linalg.LinAlgError) # check that args with differing dtypes raise def assert_homogeneous_dtypes(self, name, cfunc, args): msg = "np.linalg.%s() only supports inputs that have homogeneous dtypes." % name self.assert_error(cfunc, args, msg, errors.TypingError) class TestLinalgLstsq(TestLinalgSystems): """ Tests for np.linalg.lstsq. """ # NOTE: The testing of this routine is hard as it has to handle numpy # using double precision routines on single precision input, this has # a knock on effect especially in rank deficient cases and cases where # conditioning is generally poor. As a result computed ranks can differ # and consequently the calculated residual can differ. # The tests try and deal with this as best as they can through the use # of reconstruction and measures like residual norms. # Suggestions for improvements are welcomed! @needs_lapack def test_linalg_lstsq(self): """ Test np.linalg.lstsq """ cfunc = jit(nopython=True)(lstsq_system) def check(A, B, **kwargs): expected = lstsq_system(A, B, **kwargs) got = cfunc(A, B, **kwargs) # check that the returned tuple is same length self.assertEqual(len(expected), len(got)) # and that length is 4 self.assertEqual(len(got), 4) # and that the computed results are contig and in the same way self.assert_contig_sanity(got, "C") use_reconstruction = False # check the ranks are the same and continue to a standard # match if that is the case (if ranks differ, then output # in e.g. residual array is of different size!). try: self.assertEqual(got[2], expected[2]) # try plain match of each array to np first for k in range(len(expected)): try: np.testing.assert_array_almost_equal_nulp( got[k], expected[k], nulp=10) except AssertionError: # plain match failed, test by reconstruction use_reconstruction = True except AssertionError: use_reconstruction = True if use_reconstruction: x, res, rank, s = got # indicies in the output which are ndarrays out_array_idx = [0, 1, 3] try: # check the ranks are the same self.assertEqual(rank, expected[2]) # check they are dimensionally correct, skip [2] = rank. for k in out_array_idx: if isinstance(expected[k], np.ndarray): self.assertEqual(got[k].shape, expected[k].shape) except AssertionError: # check the rank differs by 1. (numerical fuzz) self.assertTrue(abs(rank - expected[2]) < 2) # check if A*X = B resolution = np.finfo(A.dtype).resolution try: # this will work so long as the conditioning is # ok and the rank is full rec = np.dot(A, x) np.testing.assert_allclose( B, rec, rtol=10 * resolution, atol=10 * resolution ) except AssertionError: # system is probably under/over determined and/or # poorly conditioned. Check slackened equality # and that the residual norm is the same. for k in out_array_idx: try: np.testing.assert_allclose( expected[k], got[k], rtol=100 * resolution, atol=100 * resolution ) except AssertionError: # check the fail is likely due to bad conditioning c = np.linalg.cond(A) self.assertGreater(10 * c, (1. / resolution)) # make sure the residual 2-norm is ok # if this fails its probably due to numpy using double # precision LAPACK routines for singles. res_expected = np.linalg.norm( B - np.dot(A, expected[0])) res_got = np.linalg.norm(B - np.dot(A, x)) # rtol = 10. as all the systems are products of orthonormals # and on the small side (rows, cols) < 100. np.testing.assert_allclose( res_expected, res_got, rtol=10.) # Ensure proper resource management with self.assertNoNRTLeak(): cfunc(A, B, **kwargs) # test: column vector, tall, wide, square, row vector # prime sizes, the A's sizes = [(7, 1), (11, 5), (5, 11), (3, 3), (1, 7)] # compatible B's for Ax=B must have same number of rows and 1 or more # columns # This test takes ages! So combinations are trimmed via cycling # gets a dtype cycle_dt = cycle(self.dtypes) orders = ['F', 'C'] # gets a memory order flag cycle_order = cycle(orders) # a specific condition number to use in the following tests # there is nothing special about it other than it is not magic specific_cond = 10. # inner test loop, extracted as there's additional logic etc required # that'd end up with this being repeated a lot def inner_test_loop_fn(A, dt, **kwargs): # test solve Ax=B for (column, matrix) B, same dtype as A b_sizes = (1, 13) for b_size in b_sizes: # check 2D B b_order = next(cycle_order) B = self.specific_sample_matrix( (A.shape[0], b_size), dt, b_order) check(A, B, **kwargs) # check 1D B b_order = next(cycle_order) tmp = B[:, 0].copy(order=b_order) check(A, tmp, **kwargs) # test loop for a_size in sizes: dt = next(cycle_dt) a_order = next(cycle_order) # A full rank, well conditioned system A = self.specific_sample_matrix(a_size, dt, a_order) # run the test loop inner_test_loop_fn(A, dt) m, n = a_size minmn = min(m, n) # operations that only make sense with a 2D matrix system if m != 1 and n != 1: # Test a rank deficient system r = minmn - 1 A = self.specific_sample_matrix( a_size, dt, a_order, rank=r) # run the test loop inner_test_loop_fn(A, dt) # Test a system with a given condition number for use in # testing the rcond parameter. # This works because the singular values in the # specific_sample_matrix code are linspace (1, cond, [0... if # rank deficient]) A = self.specific_sample_matrix( a_size, dt, a_order, condition=specific_cond) # run the test loop rcond = 1. / specific_cond approx_half_rank_rcond = minmn * rcond inner_test_loop_fn(A, dt, rcond=approx_half_rank_rcond) # check empty arrays empties = [ [(0, 1), (1,)], # empty A, valid b [(1, 0), (1,)], # empty A, valid b [(1, 1), (0,)], # valid A, empty 1D b [(1, 1), (1, 0)], # valid A, empty 2D b ] for A, b in empties: args = (np.empty(A), np.empty(b)) self.assert_raise_on_empty(cfunc, args) # Test input validation ok = np.array([[1., 2.], [3., 4.]], dtype=np.float64) # check ok input is ok cfunc, (ok, ok) # check bad inputs rn = "lstsq" # Wrong dtype bad = np.array([[1, 2], [3, 4]], dtype=np.int32) self.assert_wrong_dtype(rn, cfunc, (ok, bad)) self.assert_wrong_dtype(rn, cfunc, (bad, ok)) # different dtypes bad = np.array([[1, 2], [3, 4]], dtype=np.float32) self.assert_homogeneous_dtypes(rn, cfunc, (ok, bad)) self.assert_homogeneous_dtypes(rn, cfunc, (bad, ok)) # Dimension issue bad = np.array([1, 2], dtype=np.float64) self.assert_wrong_dimensions(rn, cfunc, (bad, ok)) # no nans or infs bad = np.array([[1., 2., ], [np.inf, np.nan]], dtype=np.float64) self.assert_no_nan_or_inf(cfunc, (ok, bad)) self.assert_no_nan_or_inf(cfunc, (bad, ok)) # check 1D is accepted for B (2D is done previously) # and then that anything of higher dimension raises oneD = np.array([1., 2.], dtype=np.float64) cfunc, (ok, oneD) bad = np.array([[[1, 2], [3, 4]], [[5, 6], [7, 8]]], dtype=np.float64) self.assert_wrong_dimensions_1D(rn, cfunc, (ok, bad)) # check a dimensionally invalid system raises (1D and 2D cases # checked) bad1D = np.array([1.], dtype=np.float64) bad2D = np.array([[1.], [2.], [3.]], dtype=np.float64) self.assert_dimensionally_invalid(cfunc, (ok, bad1D)) self.assert_dimensionally_invalid(cfunc, (ok, bad2D)) @needs_lapack def test_issue3368(self): X = np.array([[1., 7.54, 6.52], [1., 2.70, 4.00], [1., 2.50, 3.80], [1., 1.15, 5.64], [1., 4.22, 3.27], [1., 1.41, 5.70],], order='F') X_orig = np.copy(X) y = np.array([1., 2., 3., 4., 5., 6.]) @jit(nopython=True) def f2(X, y, test): if test: # never executed, but necessary to trigger the bug X = X[1:2, :] return np.linalg.lstsq(X, y) f2(X, y, False) np.testing.assert_allclose(X, X_orig) class TestLinalgSolve(TestLinalgSystems): """ Tests for np.linalg.solve. """ @needs_lapack def test_linalg_solve(self): """ Test np.linalg.solve """ cfunc = jit(nopython=True)(solve_system) def check(a, b, **kwargs): expected = solve_system(a, b, **kwargs) got = cfunc(a, b, **kwargs) # check that the computed results are contig and in the same way self.assert_contig_sanity(got, "F") use_reconstruction = False # try plain match of the result first try: np.testing.assert_array_almost_equal_nulp( got, expected, nulp=10) except AssertionError: # plain match failed, test by reconstruction use_reconstruction = True # If plain match fails then reconstruction is used, # this checks that AX ~= B. # Plain match can fail due to numerical fuzziness associated # with system size and conditioning, or more simply from # numpy using double precision routines for computation that # could be done in single precision (which is what numba does). # Therefore minor differences in results can appear due to # e.g. numerical roundoff being different between two precisions. if use_reconstruction: # check they are dimensionally correct self.assertEqual(got.shape, expected.shape) # check AX=B rec = np.dot(a, got) resolution = np.finfo(a.dtype).resolution np.testing.assert_allclose( b, rec, rtol=10 * resolution, atol=100 * resolution # zeros tend to be fuzzy ) # Ensure proper resource management with self.assertNoNRTLeak(): cfunc(a, b, **kwargs) # test: prime size squares sizes = [(1, 1), (3, 3), (7, 7)] # test loop for size, dtype, order in \ product(sizes, self.dtypes, 'FC'): A = self.specific_sample_matrix(size, dtype, order) b_sizes = (1, 13) for b_size, b_order in product(b_sizes, 'FC'): # check 2D B B = self.specific_sample_matrix( (A.shape[0], b_size), dtype, b_order) check(A, B) # check 1D B tmp = B[:, 0].copy(order=b_order) check(A, tmp) # check empty cfunc(np.empty((0, 0)), np.empty((0,))) # Test input validation ok = np.array([[1., 0.], [0., 1.]], dtype=np.float64) # check ok input is ok cfunc(ok, ok) # check bad inputs rn = "solve" # Wrong dtype bad = np.array([[1, 0], [0, 1]], dtype=np.int32) self.assert_wrong_dtype(rn, cfunc, (ok, bad)) self.assert_wrong_dtype(rn, cfunc, (bad, ok)) # different dtypes bad = np.array([[1, 2], [3, 4]], dtype=np.float32) self.assert_homogeneous_dtypes(rn, cfunc, (ok, bad)) self.assert_homogeneous_dtypes(rn, cfunc, (bad, ok)) # Dimension issue bad = np.array([1, 0], dtype=np.float64) self.assert_wrong_dimensions(rn, cfunc, (bad, ok)) # no nans or infs bad = np.array([[1., 0., ], [np.inf, np.nan]], dtype=np.float64) self.assert_no_nan_or_inf(cfunc, (ok, bad)) self.assert_no_nan_or_inf(cfunc, (bad, ok)) # check 1D is accepted for B (2D is done previously) # and then that anything of higher dimension raises ok_oneD = np.array([1., 2.], dtype=np.float64) cfunc(ok, ok_oneD) bad = np.array([[[1, 2], [3, 4]], [[5, 6], [7, 8]]], dtype=np.float64) self.assert_wrong_dimensions_1D(rn, cfunc, (ok, bad)) # check an invalid system raises (1D and 2D cases checked) bad1D = np.array([1.], dtype=np.float64) bad2D = np.array([[1.], [2.], [3.]], dtype=np.float64) self.assert_dimensionally_invalid(cfunc, (ok, bad1D)) self.assert_dimensionally_invalid(cfunc, (ok, bad2D)) # check that a singular system raises bad2D = self.specific_sample_matrix((2, 2), np.float64, 'C', rank=1) self.assert_raise_on_singular(cfunc, (bad2D, ok)) @needs_lapack def test_no_input_mutation(self): X = np.array([[1., 1, 1, 1], [0., 1, 1, 1], [0., 0, 1, 1], [1., 0, 0, 1],], order='F') X_orig = np.copy(X) y = np.array([1., 2., 3., 4]) y_orig = np.copy(y) @jit(nopython=True) def func(X, y, test): if test: # not executed, triggers A order in X X = X[1:2, :] return np.linalg.solve(X, y) expected = func.py_func(X, y, False) np.testing.assert_allclose(X, X_orig) np.testing.assert_allclose(y, y_orig) got = func(X, y, False) np.testing.assert_allclose(X, X_orig) np.testing.assert_allclose(y, y_orig) np.testing.assert_allclose(expected, got) class TestLinalgPinv(TestLinalgBase): """ Tests for np.linalg.pinv. """ @needs_lapack def test_linalg_pinv(self): """ Test np.linalg.pinv """ cfunc = jit(nopython=True)(pinv_matrix) def check(a, **kwargs): expected = pinv_matrix(a, **kwargs) got = cfunc(a, **kwargs) # check that the computed results are contig and in the same way self.assert_contig_sanity(got, "F") use_reconstruction = False # try plain match of each array to np first try: np.testing.assert_array_almost_equal_nulp( got, expected, nulp=10) except AssertionError: # plain match failed, test by reconstruction use_reconstruction = True # If plain match fails then reconstruction is used. # This can occur due to numpy using double precision # LAPACK when single can be used, this creates round off # problems. Also, if the matrix has machine precision level # zeros in its singular values then the singular vectors are # likely to vary depending on round off. if use_reconstruction: # check they are dimensionally correct self.assertEqual(got.shape, expected.shape) # check pinv(A)*A~=eye # if the problem is numerical fuzz then this will probably # work, if the problem is rank deficiency then it won't! rec = np.dot(got, a) try: self.assert_is_identity_matrix(rec) except AssertionError: # check A=pinv(pinv(A)) resolution = 5 * np.finfo(a.dtype).resolution rec = cfunc(got) np.testing.assert_allclose( rec, a, rtol=10 * resolution, atol=100 * resolution # zeros tend to be fuzzy ) if a.shape[0] >= a.shape[1]: # if it is overdetermined or fully determined # use numba lstsq function (which is type specific) to # compute the inverse and check against that. lstsq = jit(nopython=True)(lstsq_system) lstsq_pinv = lstsq( a, np.eye( a.shape[0]).astype( a.dtype), **kwargs)[0] np.testing.assert_allclose( got, lstsq_pinv, rtol=10 * resolution, atol=100 * resolution # zeros tend to be fuzzy ) # check the 2 norm of the difference is small self.assertLess(np.linalg.norm(got - expected), resolution) # Ensure proper resource management with self.assertNoNRTLeak(): cfunc(a, **kwargs) # test: column vector, tall, wide, square, row vector # prime sizes sizes = [(7, 1), (11, 5), (5, 11), (3, 3), (1, 7)] # When required, a specified condition number specific_cond = 10. # test loop for size, dtype, order in \ product(sizes, self.dtypes, 'FC'): # check a full rank matrix a = self.specific_sample_matrix(size, dtype, order) check(a) m, n = size if m != 1 and n != 1: # check a rank deficient matrix minmn = min(m, n) a = self.specific_sample_matrix(size, dtype, order, condition=specific_cond) rcond = 1. / specific_cond approx_half_rank_rcond = minmn * rcond check(a, rcond=approx_half_rank_rcond) # check empty for sz in [(0, 1), (1, 0)]: check(np.empty(sz)) rn = "pinv" # Wrong dtype self.assert_wrong_dtype(rn, cfunc, (np.ones((2, 2), dtype=np.int32),)) # Dimension issue self.assert_wrong_dimensions(rn, cfunc, (np.ones(10, dtype=np.float64),)) # no nans or infs self.assert_no_nan_or_inf(cfunc, (np.array([[1., 2., ], [np.inf, np.nan]], dtype=np.float64),)) @needs_lapack def test_issue5870(self): # testing for mutation of input matrix @jit(nopython=True) def some_fn(v): return np.linalg.pinv(v[0]) v_data = np.array([[1., 3, 2, 7,], [-5, 4, 2, 3,], [9, -3, 1, 1,], [2, -2, 2, 8,]], order='F') v_orig = np.copy(v_data) reshaped_v = v_data.reshape((1, 4, 4)) expected = some_fn.py_func(reshaped_v) np.testing.assert_allclose(v_data, v_orig) got = some_fn(reshaped_v) np.testing.assert_allclose(v_data, v_orig) np.testing.assert_allclose(expected, got) class TestLinalgDetAndSlogdet(TestLinalgBase): """ Tests for np.linalg.det. and np.linalg.slogdet. Exactly the same inputs are used for both tests as det() is a trivial function of slogdet(), the tests are therefore combined. """ def check_det(self, cfunc, a, **kwargs): expected = det_matrix(a, **kwargs) got = cfunc(a, **kwargs) resolution = 5 * np.finfo(a.dtype).resolution # check the determinants are the same np.testing.assert_allclose(got, expected, rtol=resolution) # Ensure proper resource management with self.assertNoNRTLeak(): cfunc(a, **kwargs) def check_slogdet(self, cfunc, a, **kwargs): expected = slogdet_matrix(a, **kwargs) got = cfunc(a, **kwargs) # As numba returns python floats types and numpy returns # numpy float types, some more adjustment and different # types of comparison than those used with array based # results is required. # check that the returned tuple is same length self.assertEqual(len(expected), len(got)) # and that length is 2 self.assertEqual(len(got), 2) # check that the domain of the results match for k in range(2): self.assertEqual( np.iscomplexobj(got[k]), np.iscomplexobj(expected[k])) # turn got[0] into the same dtype as `a` # this is so checking with nulp will work got_conv = a.dtype.type(got[0]) np.testing.assert_array_almost_equal_nulp( got_conv, expected[0], nulp=10) # compare log determinant magnitude with a more fuzzy value # as numpy values come from higher precision lapack routines resolution = 5 * np.finfo(a.dtype).resolution np.testing.assert_allclose( got[1], expected[1], rtol=resolution, atol=resolution) # Ensure proper resource management with self.assertNoNRTLeak(): cfunc(a, **kwargs) def do_test(self, rn, check, cfunc): # test: 1x1 as it is unusual, 4x4 as it is even and 7x7 as is it odd! sizes = [(1, 1), (4, 4), (7, 7)] # test loop for size, dtype, order in \ product(sizes, self.dtypes, 'FC'): # check a full rank matrix a = self.specific_sample_matrix(size, dtype, order) check(cfunc, a) # use a matrix of zeros to trip xgetrf U(i,i)=0 singular test for dtype, order in product(self.dtypes, 'FC'): a = np.zeros((3, 3), dtype=dtype) check(cfunc, a) # check empty check(cfunc, np.empty((0, 0))) # Wrong dtype self.assert_wrong_dtype(rn, cfunc, (np.ones((2, 2), dtype=np.int32),)) # Dimension issue self.assert_wrong_dimensions(rn, cfunc, (np.ones(10, dtype=np.float64),)) # no nans or infs self.assert_no_nan_or_inf(cfunc, (np.array([[1., 2., ], [np.inf, np.nan]], dtype=np.float64),)) @needs_lapack def test_linalg_det(self): cfunc = jit(nopython=True)(det_matrix) self.do_test("det", self.check_det, cfunc) @needs_lapack def test_linalg_slogdet(self): cfunc = jit(nopython=True)(slogdet_matrix) self.do_test("slogdet", self.check_slogdet, cfunc) @needs_lapack def test_no_input_mutation(self): X = np.array([[1., 3, 2, 7,], [-5, 4, 2, 3,], [9, -3, 1, 1,], [2, -2, 2, 8,]], order='F') X_orig = np.copy(X) @jit(nopython=True) def func(X, test): if test: # not executed, but necessary to trigger A ordering in X X = X[1:2, :] return np.linalg.slogdet(X) expected = func.py_func(X, False) np.testing.assert_allclose(X, X_orig) got = func(X, False) np.testing.assert_allclose(X, X_orig) np.testing.assert_allclose(expected, got) # Use TestLinalgSystems as a base to get access to additional # testing for 1 and 2D inputs. class TestLinalgNorm(TestLinalgSystems): """ Tests for np.linalg.norm. """ @needs_lapack def test_linalg_norm(self): """ Test np.linalg.norm """ cfunc = jit(nopython=True)(norm_matrix) def check(a, **kwargs): expected = norm_matrix(a, **kwargs) got = cfunc(a, **kwargs) # All results should be in the real domain self.assertTrue(not np.iscomplexobj(got)) resolution = 5 * np.finfo(a.dtype).resolution # check the norms are the same to the arg `a` precision np.testing.assert_allclose(got, expected, rtol=resolution) # Ensure proper resource management with self.assertNoNRTLeak(): cfunc(a, **kwargs) # Check 1D inputs sizes = [1, 4, 7] nrm_types = [None, np.inf, -np.inf, 0, 1, -1, 2, -2, 5, 6.7, -4.3] # standard 1D input for size, dtype, nrm_type in \ product(sizes, self.dtypes, nrm_types): a = self.sample_vector(size, dtype) check(a, ord=nrm_type) # sliced 1D input for dtype, nrm_type in \ product(self.dtypes, nrm_types): a = self.sample_vector(10, dtype)[::3] check(a, ord=nrm_type) # Check 2D inputs: # test: column vector, tall, wide, square, row vector # prime sizes sizes = [(7, 1), (11, 5), (5, 11), (3, 3), (1, 7)] nrm_types = [None, np.inf, -np.inf, 1, -1, 2, -2] # standard 2D input for size, dtype, order, nrm_type in \ product(sizes, self.dtypes, 'FC', nrm_types): # check a full rank matrix a = self.specific_sample_matrix(size, dtype, order) check(a, ord=nrm_type) # check 2D slices work for the case where xnrm2 is called from # BLAS (ord=None) to make sure it is working ok. nrm_types = [None] for dtype, nrm_type, order in \ product(self.dtypes, nrm_types, 'FC'): a = self.specific_sample_matrix((17, 13), dtype, order) # contig for C order check(a[:3], ord=nrm_type) # contig for Fortran order check(a[:, 3:], ord=nrm_type) # contig for neither order check(a[1, 4::3], ord=nrm_type) # check that numba returns zero for empty arrays. Numpy returns zero # for most norm types and raises ValueError for +/-np.inf. # there is not a great deal of consistency in Numpy's response so # it is not being emulated in Numba for dtype, nrm_type, order in \ product(self.dtypes, nrm_types, 'FC'): a = np.empty((0,), dtype=dtype, order=order) self.assertEqual(cfunc(a, nrm_type), 0.0) a = np.empty((0, 0), dtype=dtype, order=order) self.assertEqual(cfunc(a, nrm_type), 0.0) rn = "norm" # Wrong dtype self.assert_wrong_dtype(rn, cfunc, (np.ones((2, 2), dtype=np.int32),)) # Dimension issue, reuse the test from the TestLinalgSystems class self.assert_wrong_dimensions_1D( rn, cfunc, (np.ones( 12, dtype=np.float64).reshape( 2, 2, 3),)) # no nans or infs for 2d case when SVD is used (e.g 2-norm) self.assert_no_nan_or_inf(cfunc, (np.array([[1., 2.], [np.inf, np.nan]], dtype=np.float64), 2)) # assert 2D input raises for an invalid norm kind kwarg self.assert_invalid_norm_kind(cfunc, (np.array([[1., 2.], [3., 4.]], dtype=np.float64), 6)) class TestLinalgCond(TestLinalgBase): """ Tests for np.linalg.cond. """ @needs_lapack def test_linalg_cond(self): """ Test np.linalg.cond """ cfunc = jit(nopython=True)(cond_matrix) def check(a, **kwargs): expected = cond_matrix(a, **kwargs) got = cfunc(a, **kwargs) # All results should be in the real domain self.assertTrue(not np.iscomplexobj(got)) resolution = 5 * np.finfo(a.dtype).resolution # check the cond is the same to the arg `a` precision np.testing.assert_allclose(got, expected, rtol=resolution) # Ensure proper resource management with self.assertNoNRTLeak(): cfunc(a, **kwargs) # valid p values (used to indicate norm type) ps = [None, np.inf, -np.inf, 1, -1, 2, -2] sizes = [(3, 3), (7, 7)] for size, dtype, order, p in \ product(sizes, self.dtypes, 'FC', ps): a = self.specific_sample_matrix(size, dtype, order) check(a, p=p) # When p=None non-square matrices are accepted. sizes = [(7, 1), (11, 5), (5, 11), (1, 7)] for size, dtype, order in \ product(sizes, self.dtypes, 'FC'): a = self.specific_sample_matrix(size, dtype, order) check(a) # empty for sz in [(0, 1), (1, 0), (0, 0)]: self.assert_raise_on_empty(cfunc, (np.empty(sz),)) # singular systems to trip divide-by-zero x = np.array([[1, 0], [0, 0]], dtype=np.float64) check(x) check(x, p=2) x = np.array([[0, 0], [0, 0]], dtype=np.float64) check(x, p=-2) # try an ill-conditioned system with 2-norm, make sure np raises an # overflow warning as the result is `+inf` and that the result from # numba matches. with warnings.catch_warnings(): a = np.array([[1.e308, 0], [0, 0.1]], dtype=np.float64) warnings.simplefilter("ignore", RuntimeWarning) check(a) rn = "cond" # Wrong dtype self.assert_wrong_dtype(rn, cfunc, (np.ones((2, 2), dtype=np.int32),)) # Dimension issue self.assert_wrong_dimensions(rn, cfunc, (np.ones(10, dtype=np.float64),)) # no nans or infs when p="None" (default for kwarg). self.assert_no_nan_or_inf(cfunc, (np.array([[1., 2., ], [np.inf, np.nan]], dtype=np.float64),)) # assert raises for an invalid norm kind kwarg self.assert_invalid_norm_kind(cfunc, (np.array([[1., 2.], [3., 4.]], dtype=np.float64), 6)) class TestLinalgMatrixRank(TestLinalgSystems): """ Tests for np.linalg.matrix_rank. """ @needs_lapack def test_linalg_matrix_rank(self): """ Test np.linalg.matrix_rank """ cfunc = jit(nopython=True)(matrix_rank_matrix) def check(a, **kwargs): expected = matrix_rank_matrix(a, **kwargs) got = cfunc(a, **kwargs) # Ranks are integral so comparison should be trivial. # check the rank is the same np.testing.assert_allclose(got, expected) # Ensure proper resource management with self.assertNoNRTLeak(): cfunc(a, **kwargs) sizes = [(7, 1), (11, 5), (5, 11), (3, 3), (1, 7)] for size, dtype, order in \ product(sizes, self.dtypes, 'FC'): # check full rank system a = self.specific_sample_matrix(size, dtype, order) check(a) # If the system is a matrix, check rank deficiency is reported # correctly. Check all ranks from 0 to (full rank - 1). tol = 1e-13 # first check 1 to (full rank - 1) for k in range(1, min(size) - 1): # check rank k a = self.specific_sample_matrix(size, dtype, order, rank=k) self.assertEqual(cfunc(a), k) check(a) # check provision of a tolerance works as expected # create a (m x n) diagonal matrix with a singular value # guaranteed below the tolerance 1e-13 m, n = a.shape a[:, :] = 0. # reuse `a`'s memory idx = np.nonzero(np.eye(m, n)) if np.iscomplexobj(a): b = 1. + np.random.rand(k) + 1.j +\ 1.j * np.random.rand(k) # min singular value is sqrt(2)*1e-14 b[0] = 1e-14 + 1e-14j else: b = 1. + np.random.rand(k) b[0] = 1e-14 # min singular value is 1e-14 a[idx[0][:k], idx[1][:k]] = b.astype(dtype) # rank should be k-1 (as tol is present) self.assertEqual(cfunc(a, tol), k - 1) check(a, tol=tol) # then check zero rank a[:, :] = 0. self.assertEqual(cfunc(a), 0) check(a) # add in a singular value that is small if np.iscomplexobj(a): a[-1, -1] = 1e-14 + 1e-14j else: a[-1, -1] = 1e-14 # check the system has zero rank to a given tolerance self.assertEqual(cfunc(a, tol), 0) check(a, tol=tol) # check the zero vector returns rank 0 and a nonzero vector # returns rank 1. for dt in self.dtypes: a = np.zeros((5), dtype=dt) self.assertEqual(cfunc(a), 0) check(a) # make it a nonzero vector a[0] = 1. self.assertEqual(cfunc(a), 1) check(a) # empty for sz in [(0, 1), (1, 0), (0, 0)]: for tol in [None, 1e-13]: self.assert_raise_on_empty(cfunc, (np.empty(sz), tol)) rn = "matrix_rank" # Wrong dtype self.assert_wrong_dtype(rn, cfunc, (np.ones((2, 2), dtype=np.int32),)) # Dimension issue self.assert_wrong_dimensions_1D( rn, cfunc, (np.ones( 12, dtype=np.float64).reshape( 2, 2, 3),)) # no nans or infs for 2D case self.assert_no_nan_or_inf(cfunc, (np.array([[1., 2., ], [np.inf, np.nan]], dtype=np.float64),)) @needs_lapack def test_no_input_mutation(self): # this is here to test no input mutation by # numba.np.linalg._compute_singular_values # which is the workhorse for norm with 2d input, rank and cond. X = np.array([[1., 3, 2, 7,], [-5, 4, 2, 3,], [9, -3, 1, 1,], [2, -2, 2, 8,]], order='F') X_orig = np.copy(X) @jit(nopython=True) def func(X, test): if test: # not executed, but necessary to trigger A ordering in X X = X[1:2, :] return np.linalg.matrix_rank(X) expected = func.py_func(X, False) np.testing.assert_allclose(X, X_orig) got = func(X, False) np.testing.assert_allclose(X, X_orig) np.testing.assert_allclose(expected, got) class TestLinalgMatrixPower(TestLinalgBase): """ Tests for np.linalg.matrix_power. """ def assert_int_exponenent(self, cfunc, args): # validate first arg is ok cfunc(args[0], 1) # pass in both args and assert fail with self.assertRaises(errors.TypingError): cfunc(*args) @needs_lapack def test_linalg_matrix_power(self): cfunc = jit(nopython=True)(matrix_power_matrix) def check(a, pwr): expected = matrix_power_matrix(a, pwr) got = cfunc(a, pwr) # check that the computed results are contig and in the same way self.assert_contig_sanity(got, "C") res = 5 * np.finfo(a.dtype).resolution np.testing.assert_allclose(got, expected, rtol=res, atol=res) # Ensure proper resource management with self.assertNoNRTLeak(): cfunc(a, pwr) sizes = [(1, 1), (5, 5), (7, 7)] powers = [-33, -17] + list(range(-10, 10)) + [17, 33] for size, pwr, dtype, order in \ product(sizes, powers, self.dtypes, 'FC'): a = self.specific_sample_matrix(size, dtype, order) check(a, pwr) a = np.empty((0, 0), dtype=dtype, order=order) check(a, pwr) rn = "matrix_power" # Wrong dtype self.assert_wrong_dtype(rn, cfunc, (np.ones((2, 2), dtype=np.int32), 1)) # not an int power self.assert_wrong_dtype(rn, cfunc, (np.ones((2, 2), dtype=np.int32), 1)) # non square system args = (np.ones((3, 5)), 1) msg = 'input must be a square array' self.assert_error(cfunc, args, msg) # Dimension issue self.assert_wrong_dimensions(rn, cfunc, (np.ones(10, dtype=np.float64), 1)) # non-integer supplied as exponent self.assert_int_exponenent(cfunc, (np.ones((2, 2)), 1.2)) # singular matrix is not invertible self.assert_raise_on_singular(cfunc, (np.array([[0., 0], [1, 1]]), -1)) class TestTrace(TestLinalgBase): """ Tests for np.trace. """ def setUp(self): super(TestTrace, self).setUp() # compile two versions, one with and one without the offset kwarg self.cfunc_w_offset = jit(nopython=True)(trace_matrix) self.cfunc_no_offset = jit(nopython=True)(trace_matrix_no_offset) def assert_int_offset(self, cfunc, a, **kwargs): # validate first arg is ok cfunc(a) # pass in kwarg and assert fail with self.assertRaises(errors.TypingError): cfunc(a, **kwargs) def test_trace(self): def check(a, **kwargs): if 'offset' in kwargs: expected = trace_matrix(a, **kwargs) cfunc = self.cfunc_w_offset else: expected = trace_matrix_no_offset(a, **kwargs) cfunc = self.cfunc_no_offset got = cfunc(a, **kwargs) res = 5 * np.finfo(a.dtype).resolution np.testing.assert_allclose(got, expected, rtol=res, atol=res) # Ensure proper resource management with self.assertNoNRTLeak(): cfunc(a, **kwargs) # test: column vector, tall, wide, square, row vector # prime sizes sizes = [(7, 1), (11, 5), (5, 11), (3, 3), (1, 7)] # offsets to cover the range of the matrix sizes above offsets = [-13, -12, -11] + list(range(-10, 10)) + [11, 12, 13] for size, offset, dtype, order in \ product(sizes, offsets, self.dtypes, 'FC'): a = self.specific_sample_matrix(size, dtype, order) check(a, offset=offset) if offset == 0: check(a) a = np.empty((0, 0), dtype=dtype, order=order) check(a, offset=offset) if offset == 0: check(a) rn = "trace" # Dimension issue self.assert_wrong_dimensions(rn, self.cfunc_w_offset, (np.ones(10, dtype=np.float64), 1), False) self.assert_wrong_dimensions(rn, self.cfunc_no_offset, (np.ones(10, dtype=np.float64),), False) # non-integer supplied as exponent self.assert_int_offset( self.cfunc_w_offset, np.ones( (2, 2)), offset=1.2) def test_trace_w_optional_input(self): "Issue 2314" @jit("(optional(float64[:,:]),)", nopython=True) def tested(a): return np.trace(a) a = np.ones((5, 5), dtype=np.float64) tested(a) with self.assertRaises(TypeError) as raises: tested(None) errmsg = str(raises.exception) self.assertEqual('expected array(float64, 2d, A), got None', errmsg) class TestBasics(TestLinalgSystems): # TestLinalgSystems for 1d test order1 = cycle(['F', 'C', 'C', 'F']) order2 = cycle(['C', 'F', 'C', 'F']) # test: column vector, matrix, row vector, 1d sizes # (7, 1, 3) and two scalars sizes = [(7, 1), (3, 3), (1, 7), (7,), (1,), (3,), 3., 5.] def _assert_wrong_dim(self, rn, cfunc): # Dimension issue self.assert_wrong_dimensions_1D( rn, cfunc, (np.array([[[1]]], dtype=np.float64), np.ones(1)), False) self.assert_wrong_dimensions_1D( rn, cfunc, (np.ones(1), np.array([[[1]]], dtype=np.float64)), False) def _gen_input(self, size, dtype, order): if not isinstance(size, tuple): return size else: if len(size) == 1: return self.sample_vector(size[0], dtype) else: return self.sample_vector( size[0] * size[1], dtype).reshape( size, order=order) def _get_input(self, size1, size2, dtype): a = self._gen_input(size1, dtype, next(self.order1)) b = self._gen_input(size2, dtype, next(self.order2)) # force domain consistency as underlying ufuncs require it if np.iscomplexobj(a): b = b + 1j if np.iscomplexobj(b): a = a + 1j return (a, b) def test_outer(self): cfunc = jit(nopython=True)(outer_matrix) def check(a, b, **kwargs): # check without kwargs expected = outer_matrix(a, b) got = cfunc(a, b) res = 5 * np.finfo(np.asarray(a).dtype).resolution np.testing.assert_allclose(got, expected, rtol=res, atol=res) # if kwargs present check with them too if 'out' in kwargs: got = cfunc(a, b, **kwargs) np.testing.assert_allclose(got, expected, rtol=res, atol=res) np.testing.assert_allclose(kwargs['out'], expected, rtol=res, atol=res) # Ensure proper resource management with self.assertNoNRTLeak(): cfunc(a, b, **kwargs) dts = cycle(self.dtypes) for size1, size2 in product(self.sizes, self.sizes): dtype = next(dts) (a, b) = self._get_input(size1, size2, dtype) check(a, b) c = np.empty((np.asarray(a).size, np.asarray(b).size), dtype=np.asarray(a).dtype) check(a, b, out=c) self._assert_wrong_dim("outer", cfunc) def test_kron(self): cfunc = jit(nopython=True)(kron_matrix) def check(a, b, **kwargs): expected = kron_matrix(a, b) got = cfunc(a, b) res = 5 * np.finfo(np.asarray(a).dtype).resolution np.testing.assert_allclose(got, expected, rtol=res, atol=res) # Ensure proper resource management with self.assertNoNRTLeak(): cfunc(a, b) for size1, size2, dtype in \ product(self.sizes, self.sizes, self.dtypes): (a, b) = self._get_input(size1, size2, dtype) check(a, b) self._assert_wrong_dim("kron", cfunc) args = (np.empty(10)[::2], np.empty(10)[::2]) msg = "only supports 'C' or 'F' layout" self.assert_error(cfunc, args, msg, err=errors.TypingError) class TestHelpers(TestCase): def test_copy_to_fortran_order(self): from numba.np.linalg import _copy_to_fortran_order def check(udt, expectfn, shapes, dtypes, orders): for shape, dtype, order in product(shapes, dtypes, orders): a = np.arange(np.prod(shape)).reshape(shape, order=order) r = udt(a) # check correct operation self.assertPreciseEqual(expectfn(a), r) # check new copy has made self.assertNotEqual(a.ctypes.data, r.ctypes.data) @njit def direct_call(a): return _copy_to_fortran_order(a) shapes = [(3, 4), (3, 2, 5)] dtypes = [np.intp] orders = ['C', 'F'] check(direct_call, np.asfortranarray, shapes, dtypes, orders) @njit def slice_to_any(a): # make a 'any' layout slice sliced = a[::2][0] return _copy_to_fortran_order(sliced) shapes = [(3, 3, 4), (3, 3, 2, 5)] dtypes = [np.intp] orders = ['C', 'F'] def expected_slice_to_any(a): # make a 'any' layout slice sliced = a[::2][0] return np.asfortranarray(sliced) check(slice_to_any, expected_slice_to_any, shapes, dtypes, orders) if __name__ == '__main__': unittest.main()