diff --git a/nnmnkwii/paramgen/_bandmat/misc.pyx b/nnmnkwii/paramgen/_bandmat/misc.pyx deleted file mode 100644 index b1bcd33..0000000 --- a/nnmnkwii/paramgen/_bandmat/misc.pyx +++ /dev/null @@ -1,100 +0,0 @@ -# cython: language_level=3 - -"""Assorted helpful functions.""" - -# Copyright 2013, 2014, 2015, 2016, 2017, 2018 Matt Shannon - -# This file is part of bandmat. -# See `License` for details of license and warranty. - -cimport numpy as cnp -cimport cython - -cnp.import_array() -cnp.import_ufunc() - -def fancy_plus_equals(cnp.ndarray[cnp.int_t, ndim=1] target_index_seq, - cnp.ndarray[cnp.float64_t, ndim=1] source, - cnp.ndarray[cnp.float64_t, ndim=1] target): - """Implements a += method with fancy indexing. - - Does what you might expect - target[target_index_seq] += source - to do. - """ - cdef unsigned long source_size - - source_size = source.shape[0] - assert target_index_seq.shape[0] == source_size - - cdef unsigned long source_index - cdef long target_index - - for source_index in range(source_size): - target_index = target_index_seq[source_index] - target[target_index] += source[source_index] - - return - -def fancy_plus_equals_2d(cnp.ndarray[cnp.int_t, ndim=1] target_index_seq, - cnp.ndarray[cnp.float64_t, ndim=2] source, - cnp.ndarray[cnp.float64_t, ndim=2] target): - """Implements a += method with fancy indexing. - - Does what you might expect - target[target_index_seq] += source - to do. - """ - cdef unsigned long source_size - cdef unsigned long size1 - - source_size = source.shape[0] - assert target_index_seq.shape[0] == source_size - size1 = source.shape[1] - assert target.shape[1] == size1 - - cdef unsigned long source_index - cdef long target_index - cdef unsigned long index1 - - for source_index in range(source_size): - target_index = target_index_seq[source_index] - for index1 in range(size1): - target[target_index, index1] += source[source_index, index1] - - return - -def fancy_plus_equals_3d(cnp.ndarray[cnp.int_t, ndim=1] target_index_seq, - cnp.ndarray[cnp.float64_t, ndim=3] source, - cnp.ndarray[cnp.float64_t, ndim=3] target): - """Implements a += method with fancy indexing. - - Does what you might expect - target[target_index_seq] += source - to do. - """ - cdef unsigned long source_size - cdef unsigned long size1 - cdef unsigned long size2 - - source_size = source.shape[0] - assert target_index_seq.shape[0] == source_size - size1 = source.shape[1] - assert target.shape[1] == size1 - size2 = source.shape[2] - assert target.shape[2] == size2 - - cdef unsigned long source_index - cdef long target_index - cdef unsigned long index1 - cdef unsigned long index2 - - for source_index in range(source_size): - target_index = target_index_seq[source_index] - for index1 in range(size1): - for index2 in range(size2): - target[target_index, index1, index2] += ( - source[source_index, index1, index2] - ) - - return diff --git a/nnmnkwii/paramgen/_bandmat/testhelp.py b/nnmnkwii/paramgen/_bandmat/testhelp.py new file mode 100644 index 0000000..8442cf0 --- /dev/null +++ b/nnmnkwii/paramgen/_bandmat/testhelp.py @@ -0,0 +1,79 @@ +"""Helper functions for testing.""" + +# Copyright 2013, 2014, 2015, 2016, 2017, 2018 Matt Shannon + +# This file is part of bandmat. +# See `License` for details of license and warranty. + +import numpy as np +from numpy.random import randn + + +def assert_allclose( + actual, desired, rtol=1e-7, atol=1e-14, msg="items not almost equal" +): + if np.shape(actual) != np.shape(desired): + raise AssertionError( + "%s (wrong shape)\n ACTUAL: %r\n DESIRED: %r" % (msg, actual, desired) + ) + if not np.allclose(actual, desired, rtol, atol): + abs_err = np.abs(actual - desired) + rel_err = np.abs((actual - desired) / desired) + raise AssertionError( + "%s\n ACTUAL:\n%r\n DESIRED:\n%r\n" + " ABS ERR: %r (max %s)\n REL ERR: %r (max %s)" + % (msg, actual, desired, abs_err, np.max(abs_err), rel_err, np.max(rel_err)) + ) + + +def assert_allequal(actual, desired, msg="items not equal"): + if np.shape(actual) != np.shape(desired): + raise AssertionError( + "%s (wrong shape)\n ACTUAL: %r\n DESIRED: %r" % (msg, actual, desired) + ) + if not np.all(actual == desired): + raise AssertionError("%s\n ACTUAL:\n%r\n DESIRED:\n%r" % (msg, actual, desired)) + + +def randomize_extra_entries(lo, u, mat_rect): + """Randomizes the extra entries of a rectangular matrix. + + See the docstring for `band_c` for the definition of extra entries. + + N.B. in-place, i.e. mutates `mat_rect`. + """ + assert lo >= 0 + assert u >= 0 + assert mat_rect.shape[0] == lo + u + 1 + + size = mat_rect.shape[1] + + for offset in range(1, u + 1): + mat_rect[u - offset, 0 : min(offset, size)] = randn() + for offset in range(1, lo + 1): + mat_rect[u + offset, max(size - offset, 0) : size] = randn() + + +def randomize_extra_entries_bm(mat_bm): + if mat_bm.transposed: + randomize_extra_entries(mat_bm.u, mat_bm.l, mat_bm.data) + else: + randomize_extra_entries(mat_bm.l, mat_bm.u, mat_bm.data) + + +def get_array_mem(*arrays): + """Returns a representation of the memory layout of an array. + + This is intended to be used to check whether the memory used by a given + numpy array, or how this memory is mapped into the logical indices of the + tensor it represents, changes between two points in time. + + Example usage: + >>> import numpy as np + >>> x = np.array([2.0, 3.0, 4.0]) + >>> array_mem = get_array_mem(x) + >>> # some potentially complicated operation + >>> x *= 2.0 + >>> assert get_array_mem(x) == array_mem + """ + return [array.__array_interface__ for array in arrays] diff --git a/pyproject.toml b/pyproject.toml index 7f10099..01481ca 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -24,3 +24,8 @@ mypy_ignore_packages=["nnmnkwii.paramgen._bandmat.*", "nnmnkwii.frontend.*"] [[tool.pysen.lint.mypy_targets]] paths = ["nnmnkwii", "tests"] + +[tool.pysen.lint.source] +excludes = [ + "tests/bandmat", +] diff --git a/setup.py b/setup.py index 10598ff..b74c41f 100644 --- a/setup.py +++ b/setup.py @@ -108,7 +108,6 @@ def run(self): ("nnmnkwii", "paramgen", "_bandmat", "core"), ("nnmnkwii", "paramgen", "_bandmat", "tensor"), ("nnmnkwii", "paramgen", "_bandmat", "linalg"), - ("nnmnkwii", "paramgen", "_bandmat", "misc"), ("nnmnkwii", "paramgen", "_bandmat", "overlap"), ] ext_modules.extend( diff --git a/tests/bandmat/test_core.py b/tests/bandmat/test_core.py new file mode 100644 index 0000000..f267cf3 --- /dev/null +++ b/tests/bandmat/test_core.py @@ -0,0 +1,640 @@ +"""Tests for core banded matrix definitions and functions.""" + +# Copyright 2013, 2014, 2015, 2016, 2017, 2018 Matt Shannon + +# This file is part of bandmat. +# See `License` for details of license and warranty. + +import doctest +import random +import sys +import unittest + +import numpy as np +from nnmnkwii.paramgen import _bandmat as bm +from nnmnkwii.paramgen._bandmat import full as fl +from nnmnkwii.paramgen._bandmat.testhelp import ( + assert_allclose, + assert_allequal, + get_array_mem, +) +from numpy.random import randint, randn + + +def rand_bool(): + return randint(0, 2) == 0 + + +def gen_BandMat(size, l=None, u=None, transposed=None): + """Generates a random BandMat.""" + if l is None: + l = random.choice([0, 1, randint(0, 10)]) + if u is None: + u = random.choice([0, 1, randint(0, 10)]) + data = randn(l + u + 1, size) + if transposed is None: + transposed = rand_bool() + return bm.BandMat(l, u, data, transposed=transposed) + + +def load_tests(loader, tests, ignore): + # package-level doctests (N.B. includes other modules, not just core) + tests.addTests(doctest.DocTestSuite(bm)) + tests.addTests(doctest.DocTestSuite(bm.core)) + return tests + + +class TestCore(unittest.TestCase): + def test_BandMat_basic(self, its=50): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + a_bm = gen_BandMat(size) + + assert a_bm.size == size + + def test_BandMat_full(self, its=50): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + a_bm = gen_BandMat(size) + a_full = a_bm.full() + l, u = a_bm.l, a_bm.u + + # N.B. these tests are not really testing much of anything (they + # are virtually identical to the implementation of BandMat.full), + # but this is not that surprising since the lines below are kind + # of the definition of the representation used by BandMat in the + # two cases (transposed True and transposed False). + if a_bm.transposed: + assert_allequal(a_full.T, fl.band_c(u, l, a_bm.data)) + else: + assert_allequal(a_full, fl.band_c(l, u, a_bm.data)) + + assert not np.may_share_memory(a_full, a_bm.data) + + def test_BandMat_T(self, its=50): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + a_bm = gen_BandMat(size) + + assert_allequal(a_bm.T.full(), a_bm.full().T) + + assert a_bm.T.data is a_bm.data + if size > 0: + assert np.may_share_memory(a_bm.T.data, a_bm.data) + + def test_BandMat_copy_exact(self, its=50): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + mat_bm = gen_BandMat(size) + mat_full_orig = mat_bm.full().copy() + + mat_bm_new = mat_bm.copy_exact() + assert mat_bm_new.l == mat_bm.l + assert mat_bm_new.u == mat_bm.u + assert mat_bm_new.transposed == mat_bm.transposed + + # check that copy represents the same matrix + assert_allequal(mat_bm_new.full(), mat_full_orig) + + # check that copy does not share memory with original + assert not np.may_share_memory(mat_bm_new.data, mat_bm.data) + + # check that mutating the copy does not change the original + mat_bm_new.data += 1.0 + assert_allequal(mat_bm.full(), mat_full_orig) + + def test_BandMat_copy(self, its=50): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + mat_bm = gen_BandMat(size) + mat_full_orig = mat_bm.full().copy() + + mat_bm_new = mat_bm.copy() + assert mat_bm_new.l == mat_bm.l + assert mat_bm_new.u == mat_bm.u + assert not mat_bm_new.transposed + + # check that copy represents the same matrix + assert_allequal(mat_bm_new.full(), mat_full_orig) + + # check that copy does not share memory with original + assert not np.may_share_memory(mat_bm_new.data, mat_bm.data) + + # check that mutating the copy does not change the original + mat_bm_new.data += 1.0 + assert_allequal(mat_bm.full(), mat_full_orig) + + def test_BandMat_equiv(self, its=50): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + mat_bm = gen_BandMat(size) + l_new = random.choice([None, 0, 1, randint(0, 10)]) + u_new = random.choice([None, 0, 1, randint(0, 10)]) + transposed_new = random.choice([None, True, False]) + zero_extra = rand_bool() + + l_new_value = mat_bm.l if l_new is None else l_new + u_new_value = mat_bm.u if u_new is None else u_new + transposed_new_value = ( + mat_bm.transposed if transposed_new is None else transposed_new + ) + + if l_new_value < mat_bm.l or u_new_value < mat_bm.u: + self.assertRaises( + AssertionError, + mat_bm.equiv, + l_new=l_new, + u_new=u_new, + transposed_new=transposed_new, + zero_extra=zero_extra, + ) + else: + mat_bm_new = mat_bm.equiv( + l_new=l_new, + u_new=u_new, + transposed_new=transposed_new, + zero_extra=zero_extra, + ) + assert mat_bm_new.l == l_new_value + assert mat_bm_new.u == u_new_value + assert mat_bm_new.transposed == transposed_new_value + assert_allequal(mat_bm_new.full(), mat_bm.full()) + assert not np.may_share_memory(mat_bm_new.data, mat_bm.data) + + if zero_extra: + mat_new_data_good = ( + (fl.band_e(u_new_value, l_new_value, mat_bm.full().T)) + if mat_bm_new.transposed + else (fl.band_e(l_new_value, u_new_value, mat_bm.full())) + ) + assert_allequal(mat_bm_new.data, mat_new_data_good) + + def test_BandMat_plus_equals_band_of(self, its=100): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + mult = randn() + target_bm = gen_BandMat(size) + mat_bm = gen_BandMat(size) + target_full = target_bm.full() + mat_full = mat_bm.full() + array_mem = get_array_mem(target_bm.data, mat_bm.data) + + target_bm.plus_equals_band_of(mat_bm, mult) + target_full += fl.band_ec(target_bm.l, target_bm.u, mat_full) * mult + assert_allclose(target_bm.full(), target_full) + assert get_array_mem(target_bm.data, mat_bm.data) == array_mem + + def test_BandMat_add(self, its=100): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + a_bm = gen_BandMat(size) + b_bm = gen_BandMat(size) + a_full = a_bm.full() + b_full = b_bm.full() + + c_bm = a_bm + b_bm + c_full = a_full + b_full + assert_allclose(c_bm.full(), c_full) + assert not np.may_share_memory(c_bm.data, a_bm.data) + assert not np.may_share_memory(c_bm.data, b_bm.data) + + c_bm = a_bm + 0 + c_full = a_full + 0 + assert_allclose(c_bm.full(), c_full) + assert not np.may_share_memory(c_bm.data, a_bm.data) + + c_bm = 0 + a_bm + c_full = 0 + a_full + assert_allclose(c_bm.full(), c_full) + assert not np.may_share_memory(c_bm.data, a_bm.data) + + with self.assertRaises(TypeError): + a_bm + 1.0 + with self.assertRaises(TypeError): + 1.0 + a_bm + + def test_BandMat_sum(self, its=50): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + num_terms = randint(10) + a_bms = [gen_BandMat(size) for _ in range(num_terms)] + a_fulls = [a_bm.full() for a_bm in a_bms] + + b_bm = sum(a_bms) + b_full = sum(a_fulls) + if num_terms > 0: + assert_allclose(b_bm.full(), b_full) + for a_bm in a_bms: + assert not np.may_share_memory(b_bm.data, a_bm.data) + + def test_BandMat_sub(self, its=100): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + a_bm = gen_BandMat(size) + b_bm = gen_BandMat(size) + a_full = a_bm.full() + b_full = b_bm.full() + + c_bm = a_bm - b_bm + c_full = a_full - b_full + assert_allclose(c_bm.full(), c_full) + assert not np.may_share_memory(c_bm.data, a_bm.data) + assert not np.may_share_memory(c_bm.data, b_bm.data) + + c_bm = a_bm - 0 + c_full = a_full - 0 + assert_allclose(c_bm.full(), c_full) + assert not np.may_share_memory(c_bm.data, a_bm.data) + + c_bm = 0 - a_bm + c_full = 0 - a_full + assert_allclose(c_bm.full(), c_full) + assert not np.may_share_memory(c_bm.data, a_bm.data) + + with self.assertRaises(TypeError): + a_bm - 1.0 + with self.assertRaises(TypeError): + 1.0 - a_bm + + def test_BandMat_iadd(self, its=100): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + a_bm = gen_BandMat(size) + b_bm = gen_BandMat(size) + a_full = a_bm.full() + b_full = b_bm.full() + array_mem = get_array_mem(a_bm.data, b_bm.data) + + if a_bm.l >= b_bm.l and a_bm.u >= b_bm.u: + a_bm += b_bm + a_full += b_full + assert_allclose(a_bm.full(), a_full) + assert get_array_mem(a_bm.data, b_bm.data) == array_mem + else: + with self.assertRaises(AssertionError): + a_bm += b_bm + + a_bm += 0 + a_full += 0 + assert_allclose(a_bm.full(), a_full) + assert get_array_mem(a_bm.data, b_bm.data) == array_mem + + with self.assertRaises(TypeError): + a_bm += 1.0 + + def test_BandMat_isub(self, its=100): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + a_bm = gen_BandMat(size) + b_bm = gen_BandMat(size) + a_full = a_bm.full() + b_full = b_bm.full() + array_mem = get_array_mem(a_bm.data, b_bm.data) + + if a_bm.l >= b_bm.l and a_bm.u >= b_bm.u: + a_bm -= b_bm + a_full -= b_full + assert_allclose(a_bm.full(), a_full) + assert get_array_mem(a_bm.data, b_bm.data) == array_mem + else: + with self.assertRaises(AssertionError): + a_bm -= b_bm + + a_bm -= 0 + a_full -= 0 + assert_allclose(a_bm.full(), a_full) + assert get_array_mem(a_bm.data, b_bm.data) == array_mem + + with self.assertRaises(TypeError): + a_bm -= 1.0 + + def test_BandMat_pos(self, its=100): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + a_bm = gen_BandMat(size) + a_full = a_bm.full() + + b_bm = +a_bm + b_full = +a_full + assert_allclose(b_bm.full(), b_full) + assert not np.may_share_memory(b_bm.data, a_bm.data) + + def test_BandMat_neg(self, its=100): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + a_bm = gen_BandMat(size) + a_full = a_bm.full() + + b_bm = -a_bm + b_full = -a_full + assert_allclose(b_bm.full(), b_full) + assert not np.may_share_memory(b_bm.data, a_bm.data) + + def test_BandMat_mul_and_rmul(self, its=100): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + mult = randn() + a_bm = gen_BandMat(size) + a_full = a_bm.full() + + b_bm = a_bm * mult + b_full = a_full * mult + assert b_bm.l == a_bm.l + assert b_bm.u == a_bm.u + assert_allclose(b_bm.full(), b_full) + assert not np.may_share_memory(b_bm.data, a_bm.data) + + c_bm = mult * a_bm + c_full = mult * a_full + assert c_bm.l == a_bm.l + assert c_bm.u == a_bm.u + assert_allclose(c_bm.full(), c_full) + assert not np.may_share_memory(c_bm.data, a_bm.data) + + with self.assertRaises(TypeError): + a_bm * a_bm + + def test_BandMat_various_divs(self, its=100): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + mult = randn() + a_bm = gen_BandMat(size) + a_full = a_bm.full() + + b_bm = a_bm // mult + b_full = a_full // mult + assert b_bm.l == a_bm.l + assert b_bm.u == a_bm.u + assert_allclose(b_bm.full(), b_full) + assert not np.may_share_memory(b_bm.data, a_bm.data) + + b_bm = a_bm / mult + b_full = a_full / mult + assert b_bm.l == a_bm.l + assert b_bm.u == a_bm.u + assert_allclose(b_bm.full(), b_full) + assert not np.may_share_memory(b_bm.data, a_bm.data) + + b_bm = a_bm.__floordiv__(mult) + b_full = a_full.__floordiv__(mult) + assert b_bm.l == a_bm.l + assert b_bm.u == a_bm.u + assert_allclose(b_bm.full(), b_full) + assert not np.may_share_memory(b_bm.data, a_bm.data) + + # __div__ does not exist in python3 + if sys.version_info[0] < 3: + b_bm = a_bm.__div__(mult) + b_full = a_full.__div__(mult) + assert b_bm.l == a_bm.l + assert b_bm.u == a_bm.u + assert_allclose(b_bm.full(), b_full) + assert not np.may_share_memory(b_bm.data, a_bm.data) + + b_bm = a_bm.__truediv__(mult) + b_full = a_full.__truediv__(mult) + assert b_bm.l == a_bm.l + assert b_bm.u == a_bm.u + assert_allclose(b_bm.full(), b_full) + assert not np.may_share_memory(b_bm.data, a_bm.data) + + with self.assertRaises(TypeError): + a_bm // a_bm + with self.assertRaises(TypeError): + a_bm / a_bm + with self.assertRaises(TypeError): + 1.0 // a_bm + with self.assertRaises(TypeError): + 1.0 / a_bm + + def test_BandMat_imul(self, its=100): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + mult = randn() + a_bm = gen_BandMat(size) + a_full = a_bm.full() + array_mem = get_array_mem(a_bm.data) + + a_bm *= mult + a_full *= mult + assert_allclose(a_bm.full(), a_full) + assert get_array_mem(a_bm.data) == array_mem + + with self.assertRaises(TypeError): + a_bm *= a_bm + + def test_BandMat_various_idivs(self, its=100): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + mult = randn() + + a_bm = gen_BandMat(size) + a_full = a_bm.full() + array_mem = get_array_mem(a_bm.data) + a_bm //= mult + a_full //= mult + assert_allclose(a_bm.full(), a_full) + assert get_array_mem(a_bm.data) == array_mem + + a_bm = gen_BandMat(size) + a_full = a_bm.full() + array_mem = get_array_mem(a_bm.data) + a_bm /= mult + a_full /= mult + assert_allclose(a_bm.full(), a_full) + assert get_array_mem(a_bm.data) == array_mem + + a_bm = gen_BandMat(size) + a_full = a_bm.full() + array_mem = get_array_mem(a_bm.data) + a_bm.__ifloordiv__(mult) + a_full.__ifloordiv__(mult) + assert_allclose(a_bm.full(), a_full) + assert get_array_mem(a_bm.data) == array_mem + + # __idiv__ does not exist in python3 + if sys.version_info[0] < 3: + a_bm = gen_BandMat(size) + a_full = a_bm.full() + array_mem = get_array_mem(a_bm.data) + a_bm.__idiv__(mult) + a_full.__idiv__(mult) + assert_allclose(a_bm.full(), a_full) + assert get_array_mem(a_bm.data) == array_mem + + a_bm = gen_BandMat(size) + a_full = a_bm.full() + array_mem = get_array_mem(a_bm.data) + a_bm.__itruediv__(mult) + a_full.__itruediv__(mult) + assert_allclose(a_bm.full(), a_full) + assert get_array_mem(a_bm.data) == array_mem + + with self.assertRaises(TypeError): + a_bm //= a_bm + with self.assertRaises(TypeError): + a_bm /= a_bm + + def test_BandMat_reverse_view(self, its=100): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + a_bm = gen_BandMat(size) + a_full = a_bm.full() + + b_bm = a_bm.reverse_view() + b_full = a_full[::-1, ::-1] + assert_allclose(b_bm.full(), b_full) + assert b_bm.data.base is a_bm.data + if size > 0: + assert np.may_share_memory(b_bm.data, a_bm.data) + + def test_BandMat_sub_matrix_view(self, its=100): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + start = randint(size + 1) + end = randint(size + 1) + if start > end: + start, end = end, start + a_bm = gen_BandMat(size) + a_full = a_bm.full() + + b_bm = a_bm.sub_matrix_view(start, end) + b_full = a_full[start:end, start:end] + assert_allclose(b_bm.full(), b_full) + assert b_bm.data.base is a_bm.data + if end > start: + assert np.may_share_memory(b_bm.data, a_bm.data) + + def test_BandMat_embed_as_sub_matrix(self, its=100): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + start = randint(size + 1) + end = randint(size + 1) + if start > end: + start, end = end, start + a_bm = gen_BandMat(end - start) + a_full = a_bm.full() + + b_bm = a_bm.embed_as_sub_matrix(start, size) + b_full = np.zeros((size, size)) + b_full[start:end, start:end] = a_full + assert_allclose(b_bm.full(), b_full) + assert not np.may_share_memory(b_bm.data, a_bm.data) + + def test_zeros(self, its=50): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + l = random.choice([0, 1, randint(0, 10)]) + u = random.choice([0, 1, randint(0, 10)]) + + mat_bm = bm.zeros(l, u, size) + assert mat_bm.l == l + assert mat_bm.u == u + assert_allequal(mat_bm.full(), np.zeros((size, size))) + + def test_from_full(self, its=50): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + l = random.choice([0, 1, randint(0, 10)]) + u = random.choice([0, 1, randint(0, 10)]) + mat_full = gen_BandMat(size).full() + zero_outside_band = np.all(fl.band_ec(l, u, mat_full) == mat_full) + + if zero_outside_band: + mat_bm = bm.from_full(l, u, mat_full) + assert mat_bm.l == l + assert mat_bm.u == u + assert_allequal(mat_bm.full(), mat_full) + assert not np.may_share_memory(mat_bm.data, mat_full) + else: + self.assertRaises(AssertionError, bm.from_full, l, u, mat_full) + + def test_band_c_bm(self, its=50): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + l = random.choice([0, 1, randint(0, 10)]) + u = random.choice([0, 1, randint(0, 10)]) + mat_rect = randn(l + u + 1, size) + + mat_bm = bm.band_c_bm(l, u, mat_rect) + + mat_full_good = fl.band_c(l, u, mat_rect) + assert_allequal(mat_bm.full(), mat_full_good) + assert not np.may_share_memory(mat_bm.data, mat_rect) + + def test_band_e_bm(self, its=50): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + mat_bm = gen_BandMat(size) + l = random.choice([0, 1, randint(0, 10)]) + u = random.choice([0, 1, randint(0, 10)]) + + mat_rect = bm.band_e_bm(l, u, mat_bm) + + mat_rect_good = fl.band_e(l, u, mat_bm.full()) + assert_allequal(mat_rect, mat_rect_good) + assert not np.may_share_memory(mat_rect, mat_bm.data) + + def test_band_ec_bm_view(self, its=50): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + a_bm = gen_BandMat(size) + l = random.choice([0, 1, randint(0, 10)]) + u = random.choice([0, 1, randint(0, 10)]) + + b_bm = bm.band_ec_bm_view(l, u, a_bm) + + b_full_good = fl.band_ec(l, u, a_bm.full()) + assert_allequal(b_bm.full(), b_full_good) + assert b_bm.data.base is a_bm.data + if size > 0: + assert np.may_share_memory(b_bm.data, a_bm.data) + + def test_band_ec_bm(self, its=50): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + a_bm = gen_BandMat(size) + l = random.choice([0, 1, randint(0, 10)]) + u = random.choice([0, 1, randint(0, 10)]) + + b_bm = bm.band_ec_bm(l, u, a_bm) + + b_full_good = fl.band_ec(l, u, a_bm.full()) + assert_allequal(b_bm.full(), b_full_good) + assert not np.may_share_memory(b_bm.data, a_bm.data) + + def test_band_e_bm_common(self, its=50): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + num_bms = randint(5) + mat_bms = [gen_BandMat(size) for _ in range(num_bms)] + if num_bms > 0: + l = max([mat_bm.l for mat_bm in mat_bms]) + u = max([mat_bm.u for mat_bm in mat_bms]) + + mat_rects = bm.band_e_bm_common(*mat_bms) + for mat_bm, mat_rect in zip(mat_bms, mat_rects): + assert_allclose(mat_rect, bm.band_e_bm(l, u, mat_bm)) + + def test_diag(self, its=50): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + + vec = randn(size) + mat_bm = bm.diag(vec) + assert isinstance(mat_bm, bm.BandMat) + assert_allequal(mat_bm.full(), np.diag(vec)) + assert mat_bm.data.base is vec + if size > 0: + assert np.may_share_memory(mat_bm.data, vec) + + mat_bm = gen_BandMat(size) + vec = bm.diag(mat_bm) + assert_allequal(vec, np.diag(mat_bm.full())) + assert vec.base is mat_bm.data + if size > 0: + assert np.may_share_memory(vec, mat_bm.data) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/bandmat/test_full.py b/tests/bandmat/test_full.py new file mode 100644 index 0000000..24d4d93 --- /dev/null +++ b/tests/bandmat/test_full.py @@ -0,0 +1,193 @@ +"""Tests for operations involving the bands of square matrices.""" + +# Copyright 2013, 2014, 2015, 2016, 2017, 2018 Matt Shannon + +# This file is part of bandmat. +# See `License` for details of license and warranty. + +import random +import unittest + +import numpy as np +from nnmnkwii.paramgen._bandmat import full as fl +from nnmnkwii.paramgen._bandmat.testhelp import ( + assert_allclose, + assert_allequal, + get_array_mem, +) +from numpy.random import randint, randn + + +class TestFull(unittest.TestCase): + def test_band_c_basis(self, its=100): + """Checks band_c behaves correctly on canonical basis matrices.""" + for it in range(its): + l = random.choice([0, 1, randint(0, 10)]) + u = random.choice([0, 1, randint(0, 10)]) + # size >= 1 (there are no canonical basis matrices if size == 0) + size = random.choice([1, 2, randint(1, 10), randint(1, 100)]) + + # pick a random canonical basis matrix + i = randint(-u, l + 1) + j = randint(size) + mat_rect = np.zeros((l + u + 1, size)) + mat_rect[u + i, j] = 1.0 + + mat_full = fl.band_c(l, u, mat_rect) + + k = i + j + mat_full_good = np.zeros((size, size)) + if 0 <= k < size: + mat_full_good[k, j] = 1.0 + + assert_allequal(mat_full, mat_full_good) + + def test_band_c_linear(self, its=100): + """Checks band_c is linear.""" + for it in range(its): + l = random.choice([0, 1, randint(0, 10)]) + u = random.choice([0, 1, randint(0, 10)]) + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + + # check additive + mat_rect1 = randn(l + u + 1, size) + mat_rect2 = randn(l + u + 1, size) + assert_allclose( + fl.band_c(l, u, mat_rect1 + mat_rect2), + fl.band_c(l, u, mat_rect1) + fl.band_c(l, u, mat_rect2), + ) + + # check homogeneous + mat_rect = randn(l + u + 1, size) + mult = random.choice([0.0, randn(), randn(), randn()]) + assert_allclose( + fl.band_c(l, u, mat_rect * mult), fl.band_c(l, u, mat_rect) * mult + ) + + # check output is a newly-created array + mat_full = fl.band_c(l, u, mat_rect) + assert not np.may_share_memory(mat_full, mat_rect) + + def test_band_e_basis(self, its=100): + """Checks band_e behaves correctly on canonical basis matrices.""" + for it in range(its): + l = random.choice([0, 1, randint(0, 10)]) + u = random.choice([0, 1, randint(0, 10)]) + # size >= 1 (there are no canonical basis matrices if size == 0) + size = random.choice([1, 2, randint(1, 10), randint(1, 100)]) + + # pick a random canonical basis matrix + k = randint(size) + j = randint(size) + mat_full = np.zeros((size, size)) + mat_full[k, j] = 1.0 + + mat_rect = fl.band_e(l, u, mat_full) + + i = k - j + mat_rect_good = np.zeros((l + u + 1, size)) + if -u <= i <= l: + mat_rect_good[u + i, j] = 1.0 + + assert_allequal(mat_rect, mat_rect_good) + + def test_band_e_linear(self, its=100): + """Checks band_e is linear.""" + for it in range(its): + l = random.choice([0, 1, randint(0, 10)]) + u = random.choice([0, 1, randint(0, 10)]) + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + + # check additive + mat_full1 = randn(size, size) + mat_full2 = randn(size, size) + assert_allclose( + fl.band_e(l, u, mat_full1 + mat_full2), + fl.band_e(l, u, mat_full1) + fl.band_e(l, u, mat_full2), + ) + + # check homogeneous + mat_full = randn(size, size) + mult = random.choice([0.0, randn(), randn(), randn()]) + assert_allclose( + fl.band_e(l, u, mat_full * mult), fl.band_e(l, u, mat_full) * mult + ) + + # check output is a newly-created array + mat_rect = fl.band_e(l, u, mat_full) + assert not np.may_share_memory(mat_rect, mat_full) + + def test_zero_extra_entries(self, its=100): + """Checks zero_extra_entries against its equivalent definition.""" + for it in range(its): + l = random.choice([0, 1, randint(0, 10)]) + u = random.choice([0, 1, randint(0, 10)]) + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + mat_rect = randn(l + u + 1, size) + mat_rect_good = mat_rect.copy() + array_mem = get_array_mem(mat_rect) + + fl.zero_extra_entries(l, u, mat_rect) + mat_rect_good[:] = fl.band_e(l, u, fl.band_c(l, u, mat_rect_good)) + assert_allequal(mat_rect, mat_rect_good) + assert get_array_mem(mat_rect) == array_mem + + def test_band_ce(self, its=100): + """Checks band_ce against its definition and required properties.""" + for it in range(its): + l = random.choice([0, 1, randint(0, 10)]) + u = random.choice([0, 1, randint(0, 10)]) + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + mat_rect = randn(l + u + 1, size) + + mat_rect_new = fl.band_ce(l, u, mat_rect) + mat_rect_new_good = fl.band_e(l, u, fl.band_c(l, u, mat_rect)) + assert_allequal(mat_rect_new, mat_rect_new_good) + assert not np.may_share_memory(mat_rect_new, mat_rect) + + # check idempotent + assert_allequal(fl.band_ce(l, u, mat_rect_new), mat_rect_new) + + def test_band_ec(self, its=100): + """Checks band_ec against its definition and required properties.""" + for it in range(its): + l = random.choice([0, 1, randint(0, 10)]) + u = random.choice([0, 1, randint(0, 10)]) + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + mat_full = randn(size, size) + + mat_full_new = fl.band_ec(l, u, mat_full) + mat_full_new_good = fl.band_c(l, u, fl.band_e(l, u, mat_full)) + assert_allequal(mat_full_new, mat_full_new_good) + assert not np.may_share_memory(mat_full_new, mat_full) + + # check idempotent + assert_allequal(fl.band_ec(l, u, mat_full_new), mat_full_new) + + def test_band_cTe(self, its=100): + """Checks band_cTe against its definition and required properties.""" + for it in range(its): + l = random.choice([0, 1, randint(0, 10)]) + u = random.choice([0, 1, randint(0, 10)]) + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + mat_rect = randn(l + u + 1, size) + + mat_rect_new = fl.band_cTe(l, u, mat_rect) + mat_rect_new_good = fl.band_e(u, l, fl.band_c(l, u, mat_rect).T) + assert_allequal(mat_rect_new, mat_rect_new_good) + assert not np.may_share_memory(mat_rect_new, mat_rect) + + # check a property to do with doing band_cTe twice + assert_allequal(fl.band_cTe(u, l, mat_rect_new), fl.band_ce(l, u, mat_rect)) + + # check version that uses pre-specified target + mat_rect_new2 = np.empty((l + u + 1, size)) + array_mem = get_array_mem(mat_rect, mat_rect_new2) + ret = fl.band_cTe(l, u, mat_rect, target_rect=mat_rect_new2) + self.assertIsNone(ret) + assert_allequal(mat_rect_new2, mat_rect_new) + assert get_array_mem(mat_rect, mat_rect_new2) == array_mem + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/bandmat/test_linalg.py b/tests/bandmat/test_linalg.py new file mode 100644 index 0000000..30ab8d7 --- /dev/null +++ b/tests/bandmat/test_linalg.py @@ -0,0 +1,349 @@ +"""Tests for linear algebra operations for banded matrices.""" + +# Copyright 2013, 2014, 2015, 2016, 2017, 2018 Matt Shannon + +# This file is part of bandmat. +# See `License` for details of license and warranty. + +import random +import re +import unittest + +import numpy as np +import numpy.linalg as la +import scipy.linalg as sla +from nnmnkwii.paramgen import _bandmat as bm +from nnmnkwii.paramgen._bandmat import full as fl +from nnmnkwii.paramgen._bandmat import linalg as bla +from nnmnkwii.paramgen._bandmat.testhelp import ( + assert_allclose, + randomize_extra_entries_bm, +) +from numpy.random import randint, randn +from test_core import gen_BandMat + + +def rand_bool(): + return randint(0, 2) == 0 + + +def gen_symmetric_BandMat(size, depth=None, transposed=None): + if depth is None: + depth = random.choice([0, 1, randint(0, 10)]) + if transposed is None: + transposed = rand_bool() + a_bm = gen_BandMat(size, l=depth, u=depth, transposed=transposed) + b_bm = a_bm + a_bm.T + randomize_extra_entries_bm(b_bm) + return b_bm + + +def gen_pos_def_BandMat(size, depth=None, contrib_rank=2, transposed=None): + """Generates a random positive definite BandMat.""" + assert contrib_rank >= 0 + if depth is None: + depth = random.choice([0, 1, randint(0, 10)]) + if transposed is None: + transposed = rand_bool() + mat_bm = bm.zeros(depth, depth, size) + for _ in range(contrib_rank): + diff = randint(0, depth + 1) + chol_bm = gen_BandMat(size, l=depth - diff, u=diff) + bm.dot_mm_plus_equals(chol_bm, chol_bm.T, mat_bm) + if transposed: + mat_bm = mat_bm.T + randomize_extra_entries_bm(mat_bm) + return mat_bm + + +def gen_chol_factor_BandMat(size, depth=None, contrib_rank=2, transposed=None): + """Generates a random Cholesky factor BandMat. + + This works by generating a random positive definite matrix and then + computing its Cholesky factor, since using a random matrix as a Cholesky + factor seems to often lead to ill-conditioned matrices. + """ + if transposed is None: + transposed = rand_bool() + mat_bm = gen_pos_def_BandMat(size, depth=depth, contrib_rank=contrib_rank) + chol_bm = bla.cholesky(mat_bm, lower=rand_bool()) + if transposed: + chol_bm = chol_bm.T + assert chol_bm.l == 0 or chol_bm.u == 0 + assert chol_bm.l + chol_bm.u == mat_bm.l + randomize_extra_entries_bm(chol_bm) + return chol_bm + + +class TestLinAlg(unittest.TestCase): + def test_cholesky_banded_upper_scipy_test(self): + """Basic test copied from scipy.linalg.tests.test_decomp_cholesky.""" + # Symmetric positive definite banded matrix `a` + a = np.array( + [ + [4.0, 1.0, 0.0, 0.0], + [1.0, 4.0, 0.5, 0.0], + [0.0, 0.5, 4.0, 0.2], + [0.0, 0.0, 0.2, 4.0], + ] + ) + # Banded storage form of `a`. + ab = np.array([[-1.0, 1.0, 0.5, 0.2], [4.0, 4.0, 4.0, 4.0]]) + c = bla._cholesky_banded(ab, lower=False) + ufac = np.zeros_like(a) + ufac[range(4), range(4)] = c[-1] + ufac[(0, 1, 2), (1, 2, 3)] = c[0, 1:] + assert_allclose(a, np.dot(ufac.T, ufac)) + + def test_cholesky_banded_lower_scipy_test(self): + """Basic test copied from scipy.linalg.tests.test_decomp_cholesky.""" + # Symmetric positive definite banded matrix `a` + a = np.array( + [ + [4.0, 1.0, 0.0, 0.0], + [1.0, 4.0, 0.5, 0.0], + [0.0, 0.5, 4.0, 0.2], + [0.0, 0.0, 0.2, 4.0], + ] + ) + # Banded storage form of `a`. + ab = np.array([[4.0, 4.0, 4.0, 4.0], [1.0, 0.5, 0.2, -1.0]]) + c = bla._cholesky_banded(ab, lower=True) + lfac = np.zeros_like(a) + lfac[range(4), range(4)] = c[0] + lfac[(1, 2, 3), (0, 1, 2)] = c[1, :3] + assert_allclose(a, np.dot(lfac, lfac.T)) + + def test__cholesky_banded(self, its=100): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + if rand_bool(): + mat_bm = gen_pos_def_BandMat(size, transposed=False) + else: + mat_bm = gen_symmetric_BandMat(size, transposed=False) + # make it a bit more likely to be pos def + bm.diag(mat_bm)[:] = np.abs(bm.diag(mat_bm)) + 0.1 + depth = mat_bm.l + lower = rand_bool() + if lower: + mat_half_data = mat_bm.data[depth:] + else: + mat_half_data = mat_bm.data[: (depth + 1)] + overwrite = rand_bool() + + mat_half_data_arg = mat_half_data.copy() + try: + chol_data = bla._cholesky_banded( + mat_half_data_arg, overwrite_ab=overwrite, lower=lower + ) + except la.LinAlgError as e: + # First part of the message is e.g. "2-th leading minor". + msgRe = r"^" + re.escape(str(e)[:15]) + r".*not positive definite$" + with self.assertRaisesRegexp(la.LinAlgError, msgRe): + sla.cholesky(mat_bm.full(), lower=lower) + else: + assert np.shape(chol_data) == (depth + 1, size) + if lower: + chol_bm = bm.BandMat(depth, 0, chol_data) + mat_bm_again = bm.dot_mm(chol_bm, chol_bm.T) + else: + chol_bm = bm.BandMat(0, depth, chol_data) + mat_bm_again = bm.dot_mm(chol_bm.T, chol_bm) + assert_allclose(mat_bm_again.full(), mat_bm.full()) + + if size > 0: + self.assertEqual( + np.may_share_memory(chol_data, mat_half_data_arg), overwrite + ) + + if not overwrite: + assert np.all(mat_half_data_arg == mat_half_data) + + def test__solve_triangular_banded(self, its=100): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + b = randn(size) + chol_bm = gen_chol_factor_BandMat(size, transposed=False) + chol_data = chol_bm.data + depth = chol_bm.l + chol_bm.u + lower = chol_bm.u == 0 + if size > 0 and rand_bool() and rand_bool(): + badFrame = randint(size) + chol_data[0 if lower else depth, badFrame] = 0.0 + else: + badFrame = None + transposed = rand_bool() + overwrite_b = rand_bool() + chol_full = chol_bm.full() + + b_arg = b.copy() + if badFrame is not None: + msg = "singular matrix: resolution failed at diagonal %d" % badFrame + msgRe = "^" + re.escape(msg) + "$" + with self.assertRaisesRegexp(la.LinAlgError, msgRe): + bla._solve_triangular_banded( + chol_data, + b_arg, + transposed=transposed, + lower=lower, + overwrite_b=overwrite_b, + ) + with self.assertRaisesRegexp(la.LinAlgError, msgRe): + sla.solve_triangular(chol_full, b, trans=transposed, lower=lower) + else: + x = bla._solve_triangular_banded( + chol_data, + b_arg, + transposed=transposed, + lower=lower, + overwrite_b=overwrite_b, + ) + if transposed: + assert_allclose(bm.dot_mv(chol_bm.T, x), b) + else: + assert_allclose(bm.dot_mv(chol_bm, x), b) + if size == 0: + x_good = np.zeros((size,)) + else: + x_good = sla.solve_triangular( + chol_full, b, trans=transposed, lower=lower + ) + assert_allclose(x, x_good) + assert not np.may_share_memory(x, chol_data) + if size > 0: + self.assertEqual(np.may_share_memory(x, b_arg), overwrite_b) + + if not overwrite_b: + assert np.all(b_arg == b) + + def test_cholesky(self, its=50): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + mat_bm = gen_pos_def_BandMat(size) + depth = mat_bm.l + lower = rand_bool() + alternative = rand_bool() + + chol_bm = bla.cholesky(mat_bm, lower=lower, alternative=alternative) + assert chol_bm.l == (depth if lower else 0) + assert chol_bm.u == (0 if lower else depth) + assert not np.may_share_memory(chol_bm.data, mat_bm.data) + + if lower != alternative: + mat_bm_again = bm.dot_mm(chol_bm, chol_bm.T) + else: + mat_bm_again = bm.dot_mm(chol_bm.T, chol_bm) + assert_allclose(mat_bm_again.full(), mat_bm.full()) + + def test_solve_triangular(self, its=50): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + b = randn(size) + chol_bm = gen_chol_factor_BandMat(size) + depth = chol_bm.l + chol_bm.u + lower = chol_bm.u == 0 + chol_lower_bm = chol_bm if lower else chol_bm.T + chol_full = chol_bm.full() + + x = bla.solve_triangular(chol_bm, b) + assert_allclose(bm.dot_mv(chol_bm, x), b) + if size == 0: + x_good = np.zeros((size,)) + else: + x_good = sla.solve_triangular(chol_full, b, lower=lower) + assert_allclose(x, x_good) + assert not np.may_share_memory(x, chol_bm.data) + assert not np.may_share_memory(x, b) + + def test_cho_solve(self, its=50): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + b = randn(size) + chol_bm = gen_chol_factor_BandMat(size) + depth = chol_bm.l + chol_bm.u + lower = chol_bm.u == 0 + chol_lower_bm = chol_bm if lower else chol_bm.T + chol_full = chol_bm.full() + + x = bla.cho_solve(chol_bm, b) + assert_allclose(bm.dot_mv(chol_lower_bm, bm.dot_mv(chol_lower_bm.T, x)), b) + if size == 0: + x_good = np.zeros((size,)) + else: + x_good = sla.cho_solve((chol_full, lower), b) + assert_allclose(x, x_good) + assert not np.may_share_memory(x, chol_bm.data) + assert not np.may_share_memory(x, b) + + def test_solve(self, its=50): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + b = randn(size) + # the below tries to ensure the matrix is well-conditioned + a_bm = gen_BandMat(size) + bm.diag(np.ones((size,)) * 10.0) + a_full = a_bm.full() + + x = bla.solve(a_bm, b) + assert_allclose(bm.dot_mv(a_bm, x), b) + if size == 0: + x_good = np.zeros((size,)) + else: + x_good = sla.solve(a_full, b) + assert_allclose(x, x_good) + assert not np.may_share_memory(x, a_bm.data) + assert not np.may_share_memory(x, b) + + def test_solveh(self, its=50): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + b = randn(size) + a_bm = gen_pos_def_BandMat(size) + a_full = a_bm.full() + + x = bla.solveh(a_bm, b) + assert_allclose(bm.dot_mv(a_bm, x), b) + if size == 0: + x_good = np.zeros((size,)) + else: + # x_good = sla.solve(a_full, b, sym_pos=True) + x_good = sla.solve(a_full, b, assume_a="pos") + assert_allclose(x, x_good) + assert not np.may_share_memory(x, a_bm.data) + assert not np.may_share_memory(x, b) + + def test_band_of_inverse_from_chol(self, its=50): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + chol_bm = gen_chol_factor_BandMat(size) + depth = chol_bm.l + chol_bm.u + + band_of_inv_bm = bla.band_of_inverse_from_chol(chol_bm) + assert not np.may_share_memory(band_of_inv_bm.data, chol_bm.data) + + mat_bm = ( + bm.dot_mm(chol_bm, chol_bm.T) + if chol_bm.u == 0 + else bm.dot_mm(chol_bm.T, chol_bm) + ) + band_of_inv_full_good = fl.band_ec( + depth, depth, np.eye(0, 0) if size == 0 else la.inv(mat_bm.full()) + ) + assert_allclose(band_of_inv_bm.full(), band_of_inv_full_good) + + def test_band_of_inverse(self, its=50): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + mat_bm = gen_pos_def_BandMat(size) + depth = mat_bm.l + + band_of_inv_bm = bla.band_of_inverse(mat_bm) + assert not np.may_share_memory(band_of_inv_bm.data, mat_bm.data) + + band_of_inv_full_good = fl.band_ec( + depth, depth, np.eye(0, 0) if size == 0 else la.inv(mat_bm.full()) + ) + assert_allclose(band_of_inv_bm.full(), band_of_inv_full_good) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/bandmat/test_overlap.py b/tests/bandmat/test_overlap.py new file mode 100644 index 0000000..652e786 --- /dev/null +++ b/tests/bandmat/test_overlap.py @@ -0,0 +1,262 @@ +"""Tests for functions to do with overlapping subtensors.""" + +# Copyright 2013, 2014, 2015, 2016, 2017, 2018 Matt Shannon + +# This file is part of bandmat. +# See `License` for details of license and warranty. + +import random +import unittest + +import numpy as np +from nnmnkwii.paramgen import _bandmat as bm +from nnmnkwii.paramgen._bandmat import overlap as bmo +from nnmnkwii.paramgen._bandmat.testhelp import assert_allclose, assert_allequal +from numpy.random import randint, randn +from test_core import gen_BandMat + +cc = bm.band_e_bm_common + + +def rand_bool(): + return randint(0, 2) == 0 + + +def chunk_randomly(xs): + size = len(xs) + + num_divs = random.choice([0, randint(size // 2 + 1), randint(size + 3)]) + divs = [0] + sorted([randint(size + 1) for _ in range(num_divs)]) + [size] + + for start, end in zip(divs, divs[1:]): + yield start, end, xs[start:end] + + +class TestOverlap(unittest.TestCase): + def test_sum_overlapping_v(self, its=50): + for it in range(its): + num_contribs = random.choice([0, 1, randint(10), randint(100)]) + step = random.choice([1, randint(2), randint(10)]) + width = step + random.choice([0, 1, randint(10)]) + overlap = width - step + vec_size = num_contribs * step + overlap + + contribs = randn(num_contribs, width) + target = randn(vec_size) + target_orig = target.copy() + + vec = bmo.sum_overlapping_v(contribs, step=step) + assert vec.shape == (vec_size,) + + # check target-based version adds to target correctly + bmo.sum_overlapping_v(contribs, step=step, target=target) + assert_allclose(target, target_orig + vec) + + if num_contribs == 0: + # check action for no contributions + assert_allequal(vec, np.zeros((overlap,))) + elif num_contribs == 1: + # check action for a single contribution + assert_allequal(vec, contribs[0]) + else: + # check action under splitting list of contributions in two + split_pos = randint(num_contribs + 1) + vec_again = np.zeros((vec_size,)) + bmo.sum_overlapping_v( + contribs[:split_pos], + step=step, + target=vec_again[0 : (split_pos * step + overlap)], + ) + bmo.sum_overlapping_v( + contribs[split_pos:], + step=step, + target=vec_again[(split_pos * step) : vec_size], + ) + assert_allclose(vec, vec_again) + + def test_sum_overlapping_m(self, its=50): + for it in range(its): + num_contribs = random.choice([0, 1, randint(10), randint(100)]) + step = random.choice([1, randint(2), randint(10)]) + width = max(step, 1) + random.choice([0, 1, randint(10)]) + depth = width - 1 + assert depth >= 0 + overlap = width - step + mat_size = num_contribs * step + overlap + + contribs = randn(num_contribs, width, width) + target_bm = gen_BandMat(mat_size, l=depth, u=depth) + target_bm_orig = target_bm.copy() + + mat_bm = bmo.sum_overlapping_m(contribs, step=step) + assert mat_bm.size == mat_size + assert mat_bm.l == mat_bm.u == depth + + # check target-based version adds to target_bm correctly + bmo.sum_overlapping_m(contribs, step=step, target_bm=target_bm) + assert_allclose(*cc(target_bm, target_bm_orig + mat_bm)) + + if num_contribs == 0: + # check action for no contributions + assert_allequal(mat_bm.full(), np.zeros((overlap, overlap))) + elif num_contribs == 1: + # check action for a single contribution + assert_allequal(mat_bm.full(), contribs[0]) + else: + # check action under splitting list of contributions in two + split_pos = randint(num_contribs + 1) + mat_bm_again = bm.zeros(depth, depth, mat_size) + bmo.sum_overlapping_m( + contribs[:split_pos], + step=step, + target_bm=mat_bm_again.sub_matrix_view( + 0, split_pos * step + overlap + ), + ) + bmo.sum_overlapping_m( + contribs[split_pos:], + step=step, + target_bm=mat_bm_again.sub_matrix_view(split_pos * step, mat_size), + ) + assert_allclose(*cc(mat_bm, mat_bm_again)) + + def test_extract_overlapping_v(self, its=50): + for it in range(its): + num_subs = random.choice([0, 1, randint(10), randint(100)]) + step = random.choice([1, randint(1, 10)]) + width = step + random.choice([0, 1, randint(10)]) + overlap = width - step + vec_size = num_subs * step + overlap + + vec = randn(vec_size) + target = None if rand_bool() else randn(num_subs, width) + + if target is None: + subvectors = bmo.extract_overlapping_v(vec, width, step=step) + assert subvectors.shape == (num_subs, width) + else: + bmo.extract_overlapping_v(vec, width, step=step, target=target) + subvectors = target + + for index in range(num_subs): + assert_allequal( + subvectors[index], vec[(index * step) : (index * step + width)] + ) + + def test_extract_overlapping_m(self, its=50): + for it in range(its): + num_subs = random.choice([0, 1, randint(10), randint(100)]) + step = random.choice([1, randint(1, 10)]) + width = step + random.choice([0, 1, randint(10)]) + depth = width - 1 + assert depth >= 0 + overlap = width - step + mat_size = num_subs * step + overlap + + mat_bm = gen_BandMat(mat_size, l=depth, u=depth) + target = None if rand_bool() else randn(num_subs, width, width) + + if target is None: + submats = bmo.extract_overlapping_m(mat_bm, step=step) + assert submats.shape == (num_subs, width, width) + else: + bmo.extract_overlapping_m(mat_bm, step=step, target=target) + submats = target + + for index in range(num_subs): + assert_allequal( + submats[index], + mat_bm.sub_matrix_view(index * step, index * step + width).full(), + ) + + def test_sum_overlapping_v_chunked(self, its=50): + for it in range(its): + num_contribs = random.choice([0, 1, randint(10), randint(100)]) + step = random.choice([1, randint(2), randint(10)]) + width = step + random.choice([0, 1, randint(10)]) + overlap = width - step + vec_size = num_contribs * step + overlap + + contribs = randn(num_contribs, width) + contribs_chunks = chunk_randomly(contribs) + target = randn(vec_size) + target_orig = target.copy() + + bmo.sum_overlapping_v_chunked(contribs_chunks, width, target, step=step) + vec_good = bmo.sum_overlapping_v(contribs, step=step) + assert_allclose(target, target_orig + vec_good) + + def test_sum_overlapping_m_chunked(self, its=50): + for it in range(its): + num_contribs = random.choice([0, 1, randint(10), randint(100)]) + step = random.choice([1, randint(2), randint(10)]) + width = max(step, 1) + random.choice([0, 1, randint(10)]) + depth = width - 1 + assert depth >= 0 + overlap = width - step + mat_size = num_contribs * step + overlap + + contribs = randn(num_contribs, width, width) + contribs_chunks = chunk_randomly(contribs) + target_bm = gen_BandMat(mat_size, l=depth, u=depth) + target_bm_orig = target_bm.copy() + + bmo.sum_overlapping_m_chunked(contribs_chunks, target_bm, step=step) + mat_bm_good = bmo.sum_overlapping_m(contribs, step=step) + assert_allclose(*cc(target_bm, target_bm_orig + mat_bm_good)) + + def test_extract_overlapping_v_chunked(self, its=50): + for it in range(its): + num_subs = random.choice([0, 1, randint(10), randint(100)]) + step = random.choice([1, randint(1, 10)]) + width = step + random.choice([0, 1, randint(10)]) + overlap = width - step + vec_size = num_subs * step + overlap + chunk_size = random.choice([1, randint(1, 10), randint(1, 10)]) + + vec = randn(vec_size) + + indices_remaining = set(range(num_subs)) + subvectors_all = np.empty((num_subs, width)) + for start, end, subvectors in bmo.extract_overlapping_v_chunked( + vec, width, chunk_size, step=step + ): + assert end >= start + 1 + for index in range(start, end): + assert index in indices_remaining + indices_remaining.remove(index) + subvectors_all[start:end] = subvectors + + subvectors_good = bmo.extract_overlapping_v(vec, width, step=step) + assert_allclose(subvectors_all, subvectors_good) + + def test_extract_overlapping_m_chunked(self, its=50): + for it in range(its): + num_subs = random.choice([0, 1, randint(10), randint(100)]) + step = random.choice([1, randint(1, 10)]) + width = step + random.choice([0, 1, randint(10)]) + depth = width - 1 + assert depth >= 0 + overlap = width - step + mat_size = num_subs * step + overlap + chunk_size = random.choice([1, randint(1, 10), randint(1, 10)]) + + mat_bm = gen_BandMat(mat_size, l=depth, u=depth) + + indices_remaining = set(range(num_subs)) + submats_all = np.empty((num_subs, width, width)) + for start, end, submats in bmo.extract_overlapping_m_chunked( + mat_bm, chunk_size, step=step + ): + assert end >= start + 1 + for index in range(start, end): + assert index in indices_remaining + indices_remaining.remove(index) + submats_all[start:end] = submats + + submats_good = bmo.extract_overlapping_m(mat_bm, step=step) + assert_allclose(submats_all, submats_good) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/bandmat/test_tensor.py b/tests/bandmat/test_tensor.py new file mode 100644 index 0000000..a759f13 --- /dev/null +++ b/tests/bandmat/test_tensor.py @@ -0,0 +1,177 @@ +"""Tests for multiplication, etc using banded matrices.""" + +# Copyright 2013, 2014, 2015, 2016, 2017, 2018 Matt Shannon + +# This file is part of bandmat. +# See `License` for details of license and warranty. + +import random +import unittest + +import numpy as np +from nnmnkwii.paramgen import _bandmat as bm +from nnmnkwii.paramgen._bandmat import full as fl +from nnmnkwii.paramgen._bandmat.testhelp import assert_allclose, get_array_mem +from numpy.random import randint, randn +from test_core import gen_BandMat + + +def rand_bool(): + return randint(0, 2) == 0 + + +class TestTensor(unittest.TestCase): + def test_dot_mv_plus_equals(self, its=100): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + a_bm = gen_BandMat(size) + b = randn(size) + c = randn(size) + a_full = a_bm.full() + c_good = c.copy() + array_mem = get_array_mem(a_bm.data, b, c) + + bm.dot_mv_plus_equals(a_bm, b, c) + c_good += np.dot(a_full, b) + assert_allclose(c, c_good) + assert get_array_mem(a_bm.data, b, c) == array_mem + + def test_dot_mv(self, its=100): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + a_bm = gen_BandMat(size) + b = randn(size) + a_full = a_bm.full() + + c = bm.dot_mv(a_bm, b) + c_good = np.dot(a_full, b) + assert_allclose(c, c_good) + assert not np.may_share_memory(c, a_bm.data) + assert not np.may_share_memory(c, b) + + def test_dot_mm_plus_equals(self, its=50): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + a_bm = gen_BandMat(size) + b_bm = gen_BandMat(size) + c_bm = gen_BandMat(size) + diag = None if rand_bool() else randn(size) + diag_value = np.ones((size,)) if diag is None else diag + a_full = a_bm.full() + b_full = b_bm.full() + c_full = c_bm.full() + l = c_bm.l + u = c_bm.u + array_mem = get_array_mem(a_bm.data, b_bm.data, c_bm.data) + if diag is not None: + diag_mem = get_array_mem(diag) + + bm.dot_mm_plus_equals(a_bm, b_bm, c_bm, diag=diag) + c_full += fl.band_ec( + l, u, np.dot(np.dot(a_full, np.diag(diag_value)), b_full) + ) + assert_allclose(c_bm.full(), c_full) + assert get_array_mem(a_bm.data, b_bm.data, c_bm.data) == array_mem + if diag is not None: + assert get_array_mem(diag) == diag_mem + + def test_dot_mm(self, its=50): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + a_bm = gen_BandMat(size) + b_bm = gen_BandMat(size) + diag = None if rand_bool() else randn(size) + diag_value = np.ones((size,)) if diag is None else diag + a_full = a_bm.full() + b_full = b_bm.full() + + c_bm = bm.dot_mm(a_bm, b_bm, diag=diag) + c_full = np.dot(np.dot(a_full, np.diag(diag_value)), b_full) + assert c_bm.l == a_bm.l + b_bm.l + assert c_bm.u == a_bm.u + b_bm.u + assert c_bm.size == size + assert_allclose(c_bm.full(), c_full) + assert not np.may_share_memory(c_bm.data, a_bm.data) + assert not np.may_share_memory(c_bm.data, b_bm.data) + if diag is not None: + assert not np.may_share_memory(c_bm.data, diag) + + def test_dot_mm_partial(self, its=50): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + a_bm = gen_BandMat(size) + b_bm = gen_BandMat(size) + l = random.choice([0, 1, randint(0, 10)]) + u = random.choice([0, 1, randint(0, 10)]) + diag = None if rand_bool() else randn(size) + diag_value = np.ones((size,)) if diag is None else diag + a_full = a_bm.full() + b_full = b_bm.full() + + c_bm = bm.dot_mm_partial(l, u, a_bm, b_bm, diag=diag) + c_full = fl.band_ec( + l, u, np.dot(np.dot(a_full, np.diag(diag_value)), b_full) + ) + assert c_bm.l == l + assert c_bm.u == u + assert c_bm.size == size + assert_allclose(c_bm.full(), c_full) + assert not np.may_share_memory(c_bm.data, a_bm.data) + assert not np.may_share_memory(c_bm.data, b_bm.data) + if diag is not None: + assert not np.may_share_memory(c_bm.data, diag) + + def test_dot_mmm_partial(self, its=50): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + a_bm = gen_BandMat(size) + b_bm = gen_BandMat(size) + c_bm = gen_BandMat(size) + l = random.choice([0, 1, randint(0, 10)]) + u = random.choice([0, 1, randint(0, 10)]) + a_full = a_bm.full() + b_full = b_bm.full() + c_full = c_bm.full() + + d_bm = bm.dot_mmm_partial(l, u, a_bm, b_bm, c_bm) + d_full = fl.band_ec(l, u, np.dot(a_full, np.dot(b_full, c_full))) + assert d_bm.l == l + assert d_bm.u == u + assert d_bm.size == size + assert_allclose(d_bm.full(), d_full) + assert not np.may_share_memory(d_bm.data, a_bm.data) + assert not np.may_share_memory(d_bm.data, b_bm.data) + assert not np.may_share_memory(d_bm.data, c_bm.data) + + def test_band_of_outer_plus_equals(self, its=50): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + a_vec = randn(size) + b_vec = randn(size) + mult = randn() + mat_bm = gen_BandMat(size) + mat_full = mat_bm.full() + l = mat_bm.l + u = mat_bm.u + array_mem = get_array_mem(a_vec, b_vec, mat_bm.data) + + bm.band_of_outer_plus_equals(a_vec, b_vec, mat_bm, mult=mult) + mat_full += fl.band_ec(l, u, np.outer(a_vec, b_vec) * mult) + assert_allclose(mat_bm.full(), mat_full) + assert get_array_mem(a_vec, b_vec, mat_bm.data) == array_mem + + def test_trace_dot(self, its=50): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + a_bm = gen_BandMat(size) + b_bm = gen_BandMat(size) + a_full = a_bm.full() + b_full = b_bm.full() + + c = bm.trace_dot(a_bm, b_bm) + c_good = np.trace(np.dot(a_full.T, b_full)) + assert_allclose(c, c_good) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/bandmat/test_testhelp.py b/tests/bandmat/test_testhelp.py new file mode 100644 index 0000000..e4e2c00 --- /dev/null +++ b/tests/bandmat/test_testhelp.py @@ -0,0 +1,127 @@ +"""Tests for helper functions for testing.""" + +# Copyright 2013, 2014, 2015, 2016, 2017, 2018 Matt Shannon + +# This file is part of bandmat. +# See `License` for details of license and warranty. + +import doctest +import random +import unittest + +import numpy as np +from nnmnkwii.paramgen import _bandmat as bm +from nnmnkwii.paramgen._bandmat import full as fl +from nnmnkwii.paramgen._bandmat import testhelp as th +from numpy.random import randint, randn + + +def rand_bool(): + return randint(0, 2) == 0 + + +def gen_array(ranks=[0, 1, 2, 3]): + rank = random.choice(ranks) + shape = tuple([randint(5) for _ in range(rank)]) + return np.asarray(randn(*shape)) + + +def gen_BandMat_simple(size): + """Generates a random BandMat.""" + l = random.choice([0, 1, randint(0, 10)]) + u = random.choice([0, 1, randint(0, 10)]) + data = randn(l + u + 1, size) + transposed = rand_bool() + return bm.BandMat(l, u, data, transposed=transposed) + + +def load_tests(loader, tests, ignore): + tests.addTests(doctest.DocTestSuite(th)) + return tests + + +class TestTestHelp(unittest.TestCase): + def test_assert_allclose(self): + a0 = np.array([2.0, 3.0, 4.0]) + a1 = np.array([2.0, 3.0, 4.0]) + a2 = np.array([2.0, 3.0]) + a3 = np.array([2.0, 3.0, 5.0]) + a4 = np.array([[2.0, 3.0, 4.0]]) + th.assert_allclose(a0, a0) + th.assert_allclose(a0, a1) + self.assertRaises(AssertionError, th.assert_allclose, a0, a2) + self.assertRaises(AssertionError, th.assert_allclose, a0, a3) + self.assertRaises(AssertionError, th.assert_allclose, a0, a4) + + def test_assert_allequal(self): + a0 = np.array([2.0, 3.0, 4.0]) + a1 = np.array([2.0, 3.0, 4.0]) + a2 = np.array([2.0, 3.0]) + a3 = np.array([2.0, 3.0, 5.0]) + a4 = np.array([[2.0, 3.0, 4.0]]) + th.assert_allequal(a0, a0) + th.assert_allequal(a0, a1) + self.assertRaises(AssertionError, th.assert_allequal, a0, a2) + self.assertRaises(AssertionError, th.assert_allequal, a0, a3) + self.assertRaises(AssertionError, th.assert_allequal, a0, a4) + + def test_randomize_extra_entries(self, its=50): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + l = random.choice([0, 1, randint(0, 10)]) + u = random.choice([0, 1, randint(0, 10)]) + + mat_rect = randn(l + u + 1, size) + assert np.all(mat_rect != 0.0) + fl.zero_extra_entries(l, u, mat_rect) + th.randomize_extra_entries(l, u, mat_rect) + assert np.all(mat_rect != 0.0) + + mat_rect = np.zeros((l + u + 1, size)) + assert np.all(mat_rect == 0.0) + th.randomize_extra_entries(l, u, mat_rect) + fl.zero_extra_entries(l, u, mat_rect) + assert np.all(mat_rect == 0.0) + + def test_randomize_extra_entries_bm(self, its=50): + for it in range(its): + size = random.choice([0, 1, randint(0, 10), randint(0, 100)]) + mat_bm = gen_BandMat_simple(size) + + mat_full = mat_bm.full() + th.randomize_extra_entries_bm(mat_bm) + th.assert_allequal(mat_bm.full(), mat_full) + + def test_get_array_mem(self, its=50): + # (FIXME : these are not great tests, since not familiar enough with + # numpy internals to know what sorts of changes in memory layout are + # possible and how they might arise in a realistic program) + for it in range(its): + x = gen_array() + array_mem = th.get_array_mem(x) + x *= 2.0 + assert th.get_array_mem(x) == array_mem + + x = gen_array() + array_mem = th.get_array_mem(x) + x.shape = x.shape + (1,) + assert th.get_array_mem(x) != array_mem + + x = gen_array(ranks=[1, 2, 3]) + array_mem = th.get_array_mem(x) + shape = x.shape + strides = x.strides + shape_new = x.T.shape + strides_new = x.T.strides + if np.prod(shape_new) != 0: + x.shape = shape_new + x.strides = strides_new + if shape_new != shape or strides_new != strides: + # FIXME : re-enable once I understand better when this may + # fail (i.e. when memory may be unexpectedly shared). + # assert th.get_array_mem(x) != array_mem + pass + + +if __name__ == "__main__": + unittest.main()