2697 lines
94 KiB
Python
2697 lines
94 KiB
Python
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, EnableNRTStatsMixin)
|
|
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(EnableNRTStatsMixin, 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()
|
|
super(TestProduct, self).setUp()
|
|
|
|
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(EnableNRTStatsMixin, 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()
|
|
super(TestLinalgBase, self).setUp()
|
|
|
|
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.
|
|
"""
|
|
|
|
# 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 uses type-specific
|
|
# routines and therefore sometimes comes out with a different answer to
|
|
# NumPy (orthonormal bases are not unique, etc.).
|
|
|
|
def check_reconstruction(self, a, got, expected):
|
|
u, sv, vt = got
|
|
|
|
# Check they are dimensionally correct
|
|
for k in range(len(expected)):
|
|
self.assertEqual(got[k].shape, expected[k].shape)
|
|
|
|
# Columns 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
|
|
)
|
|
|
|
@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 use_reconstruction:
|
|
self.check_reconstruction(a, got, expected)
|
|
|
|
# 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)
|
|
|
|
try:
|
|
for e_a, g_a in zip(expected, got):
|
|
np.testing.assert_allclose(e_a, g_a)
|
|
except AssertionError:
|
|
self.check_reconstruction(X, got, expected)
|
|
|
|
|
|
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()
|