From 41206ea8f66939973127da18844264d2b04d29c9 Mon Sep 17 00:00:00 2001 From: "LAPTOP-6UUF43EH\\mathi" Date: Fri, 16 Jul 2021 23:38:55 +0200 Subject: [PATCH 1/4] Add identity_matrix. This is to be used as an efficient implementation of dh_dv --- src/pylogit/scipy_utils.py | 142 +++++++++++++++++++ tests/test_scipy_utils.py | 282 +++++++++++++++++++++++++++++++++++++ 2 files changed, 424 insertions(+) create mode 100644 src/pylogit/scipy_utils.py create mode 100644 tests/test_scipy_utils.py diff --git a/src/pylogit/scipy_utils.py b/src/pylogit/scipy_utils.py new file mode 100644 index 0000000..f8a63d6 --- /dev/null +++ b/src/pylogit/scipy_utils.py @@ -0,0 +1,142 @@ +# -*- coding: utf-8 -*- +""" +Created on Thu Jul 7 20:05:51 2021 + +@name: Scipy Utilities +@author: Mathijs van der Vlies +@summary: Contains an efficient implementation of an identity matrix through the + identity_matrix class, which simply performs the identity operation + when multiplied with another matrix without multiplying the individual + elements. +""" +from __future__ import absolute_import +from scipy.sparse import spmatrix +from scipy.sparse.csr import csr_matrix +import numpy as np +from scipy.sparse.sputils import check_shape, isshape + + +class identity_matrix(spmatrix): + """Efficient implementation of an identity matrix. + + When multiplied with another matrix `A`, simply passes through and returns `A` + without multiplying the individual elements. + + Parameters + ---------- + arg1 : int, 2-tuple or identity_matrix + If int `n`, then the matrix is n-by-n. + If a 2-tuple, then it must be of the form `(n, n)` (square matrix), so that + the matrix becomes n-by-n. + If of type `identity_matrix`, then its size will be copied over. + copy : bool, default = False + Inactive parameter because no data is stored by the object, kept for + consistency with `spmatrix` interface. + """ + + def __init__(self, arg1, copy=False): + super().__init__() + if isinstance(arg1, int): + self.n = arg1 + elif isinstance(arg1, type(self)): + self.n = arg1.n + elif isinstance(arg1, tuple): + if isshape(arg1): + m, n = check_shape(arg1) + if m != n: + raise ValueError("Only a square shape is valid.") + self.n = n + else: + raise ValueError("Only a square shape is valid.") + else: + raise TypeError(f"Invalid input to constructor. Expected an object of " + f"type `int` or {type(self)}") + + def set_shape(self, shape): + return super().set_shape(shape) + + def get_shape(self): + return (self.n, self.n) + + shape = property(fget=get_shape, fset=set_shape) + + def getnnz(self, axis=None): + if axis is None: + return self.n + + return 1 + + def __repr__(self): + return f"<{self.n}x{self.n} identity matrix>" + + def _mul_vector(self, other): + return other + + def _mul_multivector(self, other): + return other + + def _mul_sparse_matrix(self, other): + return other + + def transpose(self, axes=None, copy=False): + """ + Reverses the dimensions of the sparse matrix. + + Parameters + ---------- + axes : None, optional + This argument is in the signature *solely* for NumPy + compatibility reasons. Do not pass in anything except + for the default value. + copy : bool, optional + Indicates whether or not attributes of `self` should be + copied whenever possible. The degree to which attributes + are copied varies depending on the type of sparse matrix + being used. + + Returns + ------- + p : `self` with the dimensions reversed. + + See Also + -------- + numpy.matrix.transpose : NumPy's implementation of 'transpose' + for matrices + """ + if copy: + return self.copy() + + return self + + def conj(self, copy=True): + """Element-wise complex conjugation. + + If the matrix is of non-complex data type and `copy` is False, + this method does nothing and the data is not copied. + + Parameters + ---------- + copy : bool, optional + If True, the result is guaranteed to not share data with self. + + Returns + ------- + A : The element-wise complex conjugate. + + """ + return self.transpose(copy=copy) + + def tocsr(self, copy=False): + """Convert this matrix to Compressed Sparse Row format. + + For this identity matrix, `copy` doesn't affect anything + because it doesn't store any underlying data. + """ + row = np.arange(self.n) + col = np.arange(self.n) + data = np.ones(self.n) + + return csr_matrix((data, (row, col)), shape=self.shape, copy=False) + + def __pow__(self, other): + return self.copy() diff --git a/tests/test_scipy_utils.py b/tests/test_scipy_utils.py new file mode 100644 index 0000000..3332eec --- /dev/null +++ b/tests/test_scipy_utils.py @@ -0,0 +1,282 @@ +import unittest + +import numpy as np +import numpy.testing as npt +from scipy.sparse.csr import csr_matrix + +from pylogit.scipy_utils import identity_matrix + + + +def sparse_assert_equal(a1, a2): + """Assert equality of two sparse matrices""" + assert type(a1) is type(a2) + npt.assert_array_equal(a1.data, a2.data) + npt.assert_array_equal(a1.indices, a2.indices) + npt.assert_array_equal(a1.indptr, a2.indptr) + + +def identity_matrix_assert_equal(a1, a2): + """Assert equality of two identity matrices""" + assert type(a1) is type(a2) + assert a1.n == a2.n + + +class IdentityMatrixTests(unittest.TestCase): + """ + Contains the tests of the `identity_matrix` class + """ + + def setUp(self): + """Create the input data needed to test the `identity_matrix`""" + + # Set number number of rows of the identity matrix + self.n = 3 + + # Create identity matrix under test + self.I = identity_matrix(self.n) + + # Create vector to multiply with `I` + self.v = np.arange(self.n) + + # Create vector to multiply with `I` which will create an error because + # it contains too many elements + self.v_too_long = np.arange(4) + + # Create matrix to multiply with `I` + self.A = np.arange(12).reshape((self.n, 4)) + + # Create sparse matrix to multiply with `I` + row = np.array([0, 0, 1, 2, 2, 2]) + col = np.array([0, 2, 2, 0, 1, 3]) + data = np.array([1, 2, 3, 4, 5, 6]) + self.C = csr_matrix((data, (row, col)), shape=(self.n, 4)) + + def test_shape(self): + """Test that the shape of the identity matrix is as expected.""" + self.assertEqual(self.I.shape, (self.n, self.n)) + + def test_nnz(self): + """Test that the number of non-zero elements is as expected""" + self.assertEqual(self.I.getnnz(), self.n) + self.assertEqual(self.I.nnz, self.n) + + def test_mul_vector(self): + """Test that matrix-vector multiplication `I * v` returns `v` again""" + npt.assert_array_equal(self.I * self.v, self.v) + + def test_mul_vector_raises(self): + """Test that matrix-vector multiplication `I * v` with wrong-sized `v` + raises an error""" + with self.assertRaises(ValueError): + self.I * self.v_too_long + + def test_dot_vector(self): + """Test that matrix-vector multiplication `I.dot(v)` returns `v` again""" + npt.assert_array_equal(self.I.dot(self.v), self.v) + + def test_dot_vector_raises(self): + """Test that matrix-vector multiplication `I.dot(v)` with wrong-sized `v` + raises an error""" + with self.assertRaises(ValueError): + self.I.dot(self.v_too_long) + + def test_matmul_vector(self): + """Test that matrix-vector multiplication `I @ v` returns `v` again""" + npt.assert_array_equal(self.I @ self.v, self.v) + + def test_matmul_vector_raises(self): + """Test that matrix-vector multiplication `I @ v` with wrong-sized `v` + raises an error""" + with self.assertRaises(ValueError): + self.I @ self.v_too_long + + def test_rmul_vector(self): + """Test that matrix-vector multiplication `v * I` returns `v` again""" + npt.assert_array_equal(self.v * self.I, self.v) + + def test_rmul_vector_raises(self): + """Test that matrix-vector multiplication `v * I` with wrong-sized `v` + raises an error""" + with self.assertRaises(ValueError): + self.v_too_long * self.I + + def test_rmatmul_vector(self): + """Test that matrix-vector multiplication `v @ I` returns `v` again""" + npt.assert_array_equal(self.v @ self.I, self.v) + + def test_rmatmul_vector_raises(self): + """Test that matrix-vector multiplication `v @ I` with wrong-sized `v` + raises an error""" + with self.assertRaises(ValueError): + self.v_too_long @ self.I + + def test_mul_multivector(self): + """Test that matrix multiplication `I * A` returns `A` again""" + npt.assert_array_equal(self.I * self.A, self.A) + + def test_mul_multivector_raises(self): + """Test that matrix multiplication `I * A.T` with wrong-sized `A.T` + raises an error""" + with self.assertRaises(ValueError): + self.I * self.A.T + + def test_dot_multivector(self): + """Test that matrix multiplication `I * A` returns `A` again""" + npt.assert_array_equal(self.I.dot(self.A), self.A) + + def test_dot_multivector_raises(self): + """Test that matrix multiplication `I * A.T` with wrong-sized `A.T` + raises an error""" + with self.assertRaises(ValueError): + self.I.dot(self.A.T) + + def test_matmul_multivector(self): + """Test that matrix multiplication `I @ A` returns `A` again""" + npt.assert_array_equal(self.I @ self.A, self.A) + + def test_matmul_multivector_raises(self): + """Test that matrix multiplication `I @ A.T` with wrong-sized `A.T` + raises an error""" + with self.assertRaises(ValueError): + self.I @ self.A.T + + def test_rmul_multivector(self): + """Test that matrix multiplication `A.T * I` returns `A.T` again""" + At = self.A.T + npt.assert_array_equal(At * self.I, At) + + def test_rmul_multivector_raises(self): + """Test that matrix multiplication `A * I` with wrong-sized `A` + raises an error""" + with self.assertRaises(ValueError): + self.A * self.I + + def test_rmatmul_multivector(self): + """Test that matrix multiplication `A.T @ I` returns `A.T` again""" + At = self.A.T + npt.assert_array_equal(At @ self.I, At) + + def test_rmatmul_multivector_raises(self): + """Test that matrix multiplication `A @ I` with wrong-sized `A` + raises an error""" + with self.assertRaises(ValueError): + self.A @ self.I + + def test_mul_sparse_matrix(self): + """Test that sparse matrix multiplication `I * C` returns `C` again""" + sparse_assert_equal(self.I * self.C, self.C) + + def test_mul_sparse_matrix_raises(self): + """Test that sparse matrix multiplication `I * C.T` with wrong-sized `C.T` + raises an error""" + with self.assertRaises(ValueError): + self.I * self.C.T + + def test_dot_sparse_matrix(self): + """Test that sparse matrix multiplication `I.dot(C)` returns `C` again""" + sparse_assert_equal(self.I.dot(self.C), self.C) + + def test_dot_sparse_matrix_raises(self): + """Test that sparse matrix multiplication `I.dot(C.T)` with wrong-sized `C.T` + raises an error""" + with self.assertRaises(ValueError): + self.I.dot(self.C.T) + + def test_matmul_sparse_matrix(self): + """Test that sparse matrix multiplication `I @ C` returns `C` again""" + sparse_assert_equal(self.I @ self.C, self.C) + + def test_matmul_sparse_matrix_raises(self): + """Test that sparse matrix multiplication `I @ C.T` with wrong-sized `C.T` + raises an error""" + with self.assertRaises(ValueError): + self.I @ self.C.T + + # FIXME Pre-multiplying a sparse matrix with `I` currently gives an error. + # For our use-case `dh_dv.dot(design)` this doesn't matter, because this is + # of the form `I * A` with `A` being a numpy array. + # def test_rmul_sparse_matrix(self): + # """Test that sparse matrix multiplication `C.T * I` returns `C.T` again""" + # Ct = self.C.T + # npt.assert_array_equal(Ct * self.I, Ct) + + # def test_rmul_sparse_matrix_raises(self): + # """Test that sparse matrix multiplication `C * I` with wrong-sized `C` + # raises an error""" + # with self.assertRaises(ValueError): + # self.C * self.I + + # def test_rmatmul_sparse_matrix(self): + # """Test that sparse matrix multiplication `C.T @ I` returns `C.T` again""" + # Ct = self.C.T + # npt.assert_array_equal(Ct @ self.I, Ct) + + # def test_rmatmul_sparse_matrix_raises(self): + # """Test that sparse matrix multiplication `C @ I` with wrong-sized `C` + # raises an error""" + # with self.assertRaises(ValueError): + # self.C @ self.I + + def test_toarray(self): + """Test that toarray() gives an identity matrix in numpy array format""" + npt.assert_array_equal(self.I.toarray(), np.identity(3)) + + def test_tocsr(self): + """Test that tocsr() returns an identity matrix in CSR format""" + row = np.arange(self.n) + col = np.arange(self.n) + data = np.ones(self.n) + csr_expected = csr_matrix((data, (row, col)), shape=(self.n, self.n), copy=False) + + sparse_assert_equal(self.I.tocsr(), csr_expected) + + def test___pow__(self): + """Test that __pow__ method always returns the identity regardless of input""" + identity_matrix_assert_equal(self.I, self.I ** 2) + identity_matrix_assert_equal(self.I, self.I ** -1) + identity_matrix_assert_equal(self.I, self.I ** 0) + identity_matrix_assert_equal(self.I, self.I ** 1) + identity_matrix_assert_equal(self.I, self.I ** 0.5) + + def test_copy(self): + """Test that `copy` returns a copy that's not the same object, but equal to the original""" + Ic = self.I.copy() + self.assertIsNot(self.I, Ic) + identity_matrix_assert_equal(self.I, Ic) + + def test_transpose(self): + """Test that `transpose` simply returns itself again""" + self.assertIs(self.I.transpose(), self.I) + + def test_transpose_copy(self): + """Test that `transpose(copy=True)` simply returns a copy of itself""" + It = self.I.transpose(copy=True) + self.assertIsNot(self.I, It) + identity_matrix_assert_equal(self.I, It) + + def test_conj(self): + """Test that `conj` simply returns itself again""" + self.assertIs(self.I.conj(copy=False), self.I) + + def test_conj_copy(self): + """Test that `conj()` simply returns a copy of itself""" + It = self.I.conj() + self.assertIsNot(self.I, It) + identity_matrix_assert_equal(self.I, It) + + def test_conjugate(self): + """Test that `conj` simply returns itself again""" + self.assertIs(self.I.conjugate(copy=False), self.I) + + def test_conjugate_copy(self): + """Test that `conjugate()` simply returns a copy of itself""" + It = self.I.conjugate() + self.assertIsNot(self.I, It) + identity_matrix_assert_equal(self.I, It) + + def test_H(self): + """Test that `I.H` (Hermitian transpose) simply returns a copy of itself""" + Ih = self.I.H + self.assertIsNot(self.I, Ih) + identity_matrix_assert_equal(self.I, Ih) From d771f50a30c9df35c669d986da26cc0789ece087 Mon Sep 17 00:00:00 2001 From: "LAPTOP-6UUF43EH\\mathi" Date: Fri, 16 Jul 2021 23:41:09 +0200 Subject: [PATCH 2/4] Keep sparse matrix structure for `weights_per_obs` --- src/pylogit/choice_calcs.py | 5 +++-- src/pylogit/mixed_logit_calcs.py | 5 +++-- src/pylogit/nested_choice_calcs.py | 8 ++++++-- 3 files changed, 12 insertions(+), 6 deletions(-) diff --git a/src/pylogit/choice_calcs.py b/src/pylogit/choice_calcs.py index 5e24be5..fd4adda 100644 --- a/src/pylogit/choice_calcs.py +++ b/src/pylogit/choice_calcs.py @@ -927,8 +927,9 @@ def calc_fisher_info_matrix(beta, # Calculate the weights for the sample if weights is None: weights = np.ones(design.shape[0]) - weights_per_obs =\ - np.max(rows_to_obs.toarray() * weights[:, None], axis=0) + + M = rows_to_obs.multiply(weights.reshape(-1,1)) + weights_per_obs = np.max(M, axis=0).toarray().reshape(-1) ########## # Get the required matrices diff --git a/src/pylogit/mixed_logit_calcs.py b/src/pylogit/mixed_logit_calcs.py index 1973eba..c908453 100755 --- a/src/pylogit/mixed_logit_calcs.py +++ b/src/pylogit/mixed_logit_calcs.py @@ -662,8 +662,9 @@ def calc_bhhh_hessian_approximation_mixed_logit(params, # Calculate the weights for the sample if weights is None: weights = np.ones(design_3d.shape[0]) - weights_per_obs =\ - np.max(rows_to_mixers.toarray() * weights[:, None], axis=0) + + M = rows_to_obs.multiply(weights.reshape(-1,1)) + weights_per_obs = np.max(M, axis=0).toarray().reshape(-1) # Calculate the regular probability array. Note the implicit assumption # that params == index coefficients. prob_array = general_calc_probabilities(params, diff --git a/src/pylogit/nested_choice_calcs.py b/src/pylogit/nested_choice_calcs.py index c949a38..f2c125b 100644 --- a/src/pylogit/nested_choice_calcs.py +++ b/src/pylogit/nested_choice_calcs.py @@ -566,7 +566,9 @@ def calc_nested_gradient(orig_nest_coefs, # Calculate the weights for the sample if weights is None: weights = np.ones(design.shape[0]) - weights_per_obs = np.max(rows_to_obs.toarray() * weights[:, None], axis=0) + + M = rows_to_obs.multiply(weights.reshape(-1,1)) + weights_per_obs = np.max(M, axis=0).toarray().reshape(-1) # Transform the nest coefficients into their "always positive" versions nest_coefs = naturalize_nest_coefs(orig_nest_coefs) @@ -743,7 +745,9 @@ def calc_bhhh_hessian_approximation(orig_nest_coefs, # Calculate the weights for the sample if weights is None: weights = np.ones(design.shape[0]) - weights_per_obs = np.max(rows_to_obs.toarray() * weights[:, None], axis=0) + + M = rows_to_obs.multiply(weights.reshape(-1,1)) + weights_per_obs = np.max(M, axis=0).toarray().reshape(-1) # Transform the nest coefficients into their "always positive" versions nest_coefs = naturalize_nest_coefs(orig_nest_coefs) From 2c395b35e1dc99aed206ef43df144627824f18fb Mon Sep 17 00:00:00 2001 From: "LAPTOP-6UUF43EH\\mathi" Date: Fri, 16 Jul 2021 23:42:38 +0200 Subject: [PATCH 3/4] Use identity_matrix rather than csr_matrix for dh_dv This means that dh_dv.dot(design) will be an instant calculation. --- src/pylogit/conditional_logit.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/pylogit/conditional_logit.py b/src/pylogit/conditional_logit.py index 9688b26..add9d81 100644 --- a/src/pylogit/conditional_logit.py +++ b/src/pylogit/conditional_logit.py @@ -16,6 +16,8 @@ import numpy as np from scipy.sparse import diags +from pylogit.scipy_utils import identity_matrix + from . import choice_calcs as cc from . import base_multinomial_cm_v2 as base_mcm from .estimation import LogitTypeEstimator @@ -177,7 +179,7 @@ class MNLEstimator(LogitTypeEstimator): def set_derivatives(self): # Pre-calculate the derivative of the transformation vector with # respect to the vector of systematic utilities - dh_dv = diags(np.ones(self.design.shape[0]), 0, format='csr') + dh_dv = identity_matrix(self.design.shape[0]) # Create a function to calculate dh_dv which will return the # pre-calculated result when called From 2e7a06907d11b6fe02d3f3f9df91d374ed8a0c6d Mon Sep 17 00:00:00 2001 From: "LAPTOP-6UUF43EH\\mathi" Date: Mon, 19 Jul 2021 22:46:50 +0200 Subject: [PATCH 4/4] Add tests of identity_matrix constructor In the process, add checking for negative `arg1` argument. --- src/pylogit/scipy_utils.py | 5 +- tests/test_scipy_utils.py | 97 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 101 insertions(+), 1 deletion(-) diff --git a/src/pylogit/scipy_utils.py b/src/pylogit/scipy_utils.py index f8a63d6..077e1e7 100644 --- a/src/pylogit/scipy_utils.py +++ b/src/pylogit/scipy_utils.py @@ -37,6 +37,9 @@ class identity_matrix(spmatrix): def __init__(self, arg1, copy=False): super().__init__() if isinstance(arg1, int): + if arg1 < 0: + raise ValueError(f"The number of rows `n` cannot be negative. " + f"Got a value of {arg1}") self.n = arg1 elif isinstance(arg1, type(self)): self.n = arg1.n @@ -50,7 +53,7 @@ def __init__(self, arg1, copy=False): raise ValueError("Only a square shape is valid.") else: raise TypeError(f"Invalid input to constructor. Expected an object of " - f"type `int` or {type(self)}") + f"type `int`, `tuple`, or {type(self)}") def set_shape(self, shape): return super().set_shape(shape) diff --git a/tests/test_scipy_utils.py b/tests/test_scipy_utils.py index 3332eec..4a3228b 100644 --- a/tests/test_scipy_utils.py +++ b/tests/test_scipy_utils.py @@ -22,6 +22,103 @@ def identity_matrix_assert_equal(a1, a2): assert a1.n == a2.n +class IdentityMatrixConstructorTests(unittest.TestCase): + """ + Contains the tests for the `identity_matrix` constructor + """ + + def setUp(self): + """Create the input data needed to test the `identity_matrix`""" + + # Set number rows/columns of an identity matrix + self.n = 3 + + # Set negative number of rows which will throw an error + self.neg_n = -3 + + # Set shape for an identity matrix + self.shape = (3, 3) + + # Set non-square shape which will lead to a ValueError + self.non_square_shape = (3, 4) + + # Set wrong shape not consisting of ints + self.string_shape = ('hello', 'world') + + # Set wrong shape not consisting of ints but whose elements are the same + # If the elements are the same and the type is not checked, the code + # will act as if it's a square shape because the elements are the same. + self.string_square_shape = ('hello', 'hello') + + # Set 1-element tuple which will raise a ValueError because the tuple has to contain two elements + self.tuple1 = (1,) + + # Set 3-element tuple which will raise a ValueError because the tuple has to contain two elements + self.tuple3 = (1, 2, 3) + + # Set identity_matrix which can be used to create another identity matrix with the same shape + self.I = identity_matrix(2) + + # Set arbitrary object type not supported by `identity_matrix` + self.wrong_arg1 = object() + + def test_construct_int(self): + """Test that constructor works for a valid `int` input""" + I = identity_matrix(self.n) + identity_matrix_assert_equal(I, identity_matrix(self.n)) + + def test_construct_zero(self): + """Test that constructor works for `0` input (this is a trivial 0-by-0 matrix)""" + I = identity_matrix(self.n) + identity_matrix_assert_equal(I, identity_matrix(self.n)) + + def test_construct_neg_int_raises(self): + """Test that constructor throws ValueError upon a negative `int` input""" + with self.assertRaises(ValueError): + I = identity_matrix(self.neg_n) + + def test_construct_shape(self): + """Test that constructor works for a valid square `shape` input""" + I = identity_matrix(self.shape) + identity_matrix_assert_equal(I, identity_matrix(self.shape[0])) + + def test_construct_non_square_shape_raises(self): + """Test that constructor throws a ValueError upon a non-square `shape` input""" + with self.assertRaises(ValueError): + I = identity_matrix(self.non_square_shape) + + def test_construct_non_int_shape(self): + """Test that constructor throws a ValueError upon non-int elements in a `shape` input""" + with self.assertRaises(ValueError): + I = identity_matrix(self.string_shape) + + def test_construct_non_int_square_shape(self): + """Test that constructor throws a ValueError upon non-int elements in a `shape` input, + also when the `shape` is 'square' because the two elements are equal""" + with self.assertRaises(ValueError): + I = identity_matrix(self.string_square_shape) + + def test_construct_1_tuple(self): + """Test that constructor throws a ValueError when a 1-tuple is entered rather than a 2-tuple""" + with self.assertRaises(ValueError): + I = identity_matrix(self.tuple1) + + def test_construct_3_tuple(self): + """Test that constructor throws a ValueError when a 3-tuple is entered rather than a 2-tuple""" + with self.assertRaises(ValueError): + I = identity_matrix(self.tuple3) + + def test_construct_identity_matrix(self): + """Test that constructor works when another `identity_matrix` is entered""" + I = identity_matrix(self.I) + identity_matrix_assert_equal(I, self.I) + + def test_construct_wrong_type(self): + """Test that constructor throws a TypeError when a wrong type is entered""" + with self.assertRaises(TypeError): + I = identity_matrix(self.wrong_arg1) + + class IdentityMatrixTests(unittest.TestCase): """ Contains the tests of the `identity_matrix` class