diff --git a/src/atmos_flux_inversion/correlations.py b/src/atmos_flux_inversion/correlations.py index e0ad9fd..d474971 100644 --- a/src/atmos_flux_inversion/correlations.py +++ b/src/atmos_flux_inversion/correlations.py @@ -17,10 +17,11 @@ from numpy import fromfunction, asarray, hstack, flip from numpy import exp, square, fmin, sqrt from numpy import logical_or, concatenate, isnan -from numpy import where +from numpy import where, prod from numpy import sum as array_sum from scipy.special import gamma, kv as K_nu from scipy.sparse.linalg.interface import LinearOperator +from scipy.fftpack import dctn, fftn import pyfftw.interfaces.cache from pyfftw import next_fast_len @@ -284,8 +285,12 @@ def from_array(cls, corr_array, is_cyclic=True): self = cls(shape, computational_shape) for axis in reversed(range(ndims)): + sub_index = [ + slice(None) for i in range(ndims) + ] + sub_index[axis] = slice(1, -1) corr_array = concatenate( - [corr_array, flip(corr_array[1:-1], axis)], + [corr_array, flip(corr_array[tuple(sub_index)], axis)], axis=axis) # Advantages over dctn: guaranteed same format and gets @@ -439,6 +444,22 @@ def kron(self, other): other._fourier_near_zero[other_index]) return newinst + def det(self): + """Find the determinant of the operator. + + Returns + ------- + float + """ + correlations = self._ifft(self._corr_fourier) + if self._is_cyclic: + spectrum = fftn(correlations, shape=self._underlying_shape).real + # The order is different from la.eigvalsh, but that + # doesn't matter + return prod(spectrum) + spectrum = dctn(correlations, type=1, shape=self._underlying_shape) + return prod(spectrum) + def make_matrix(corr_func, shape): """Make a correlation matrix for a domain with shape `shape`. diff --git a/src/atmos_flux_inversion/linalg.py b/src/atmos_flux_inversion/linalg.py index ba9e2c8..d7430bf 100644 --- a/src/atmos_flux_inversion/linalg.py +++ b/src/atmos_flux_inversion/linalg.py @@ -14,7 +14,7 @@ # implicitly called on the operator2.dot(chunk) that usually follows # this. from numpy import einsum -from numpy import concatenate, zeros, nonzero +from numpy import concatenate, zeros, nonzero, prod from numpy import asarray, atleast_2d, stack, where, sqrt from scipy.sparse.linalg import lgmres @@ -278,6 +278,24 @@ def matrix_sqrt(mat): cls=type(mat))) +def matrix_determinant(operator): + """Find the determinant of the operator. + + Parameters + ---------- + operator: LinearOperator + + Returns + ------- + determinant: float + """ + if hasattr(operator, "det"): + return operator.det() + if isinstance(operator, DaskMatrixLinearOperator): + operator = operator.A + return la.det(operator) + + class DaskKroneckerProductOperator(DaskLinearOperator): """Operator for Kronecker product. @@ -491,6 +509,25 @@ def quadratic_form(self, mat): operator2.dot(chunk)) return result + def det(self): + """Find the determinant of the operator. + + Returns + ------- + float + """ + op1 = self._operator1 + op2 = self._operator2 + if ( + self.shape[0] != self.shape[1] or + op1.shape[0] != op1.shape[1] + ): + raise ValueError("Determinant only defined for square operators.") + return ( + matrix_determinant(op1) ** op2.shape[0] * + matrix_determinant(op2) ** op1.shape[0] + ) + class SchmidtKroneckerProduct(DaskLinearOperator): """Kronecker product of two operators using Schmidt decomposition. @@ -580,6 +617,25 @@ def _matvec(self, vector): return asarray(result) + def det(self): + """Find the determinant of the operator. + + Returns + ------- + float + """ + op1 = self._operator1 + op2 = self._operator2 + if ( + self.shape[0] != self.shape[1] or + op1.shape[0] != op1.shape[1] + ): + raise ValueError("Determinant only defined for square operators.") + return ( + matrix_determinant(op1) ** op2.shape[0] * + matrix_determinant(op2) ** op1.shape[0] + ) + class SelfAdjointLinearOperator(DaskLinearOperator): """Self-adjoint linear operators. @@ -705,5 +761,19 @@ def solve(self, vector): return where(self._diag_near_zero, 0, result) def sqrt(self): - """Find S such that S.T @ S == self.""" + """Find S such that S.T @ S == self. + + Returns + ------- + LinearOperator + """ return DiagonalOperator(sqrt(self._diag)) + + def det(self): + """Find the determinant of the operator. + + Returns + ------- + float + """ + return prod(self._diag) diff --git a/src/atmos_flux_inversion/tests.py b/src/atmos_flux_inversion/tests.py index 7f8edac..e593e0f 100644 --- a/src/atmos_flux_inversion/tests.py +++ b/src/atmos_flux_inversion/tests.py @@ -1222,6 +1222,16 @@ def test_acyclic_from_array(self): np_tst.assert_allclose(op.dot(np.eye(*mat.shape)), mat) + def test_acyclic_from_array2d(self): + """Test from_array with correlations assumed acyclic.""" + array = [[1, .5, .25, .125, .0625, .03125], + [.5, .25, .125, .0625, .03125, .015625]] + # At one point this operation crashed, so this is an + # "assertDoesNotRaise" + (atmos_flux_inversion.correlations. + HomogeneousIsotropicCorrelation. + from_array(array, False)) + class TestSchmidtKroneckerProduct(unittest2.TestCase): """Test the Schmidt Kronecker product implementation for LinearOperators. @@ -3163,5 +3173,155 @@ def test_simple_variances(self): np_tst.assert_allclose(corr_matrix, np.diag(series.values)) +class TestDeterminants(unittest2.TestCase): + """Test the determinant-finding abilities of operators.""" + + def test_matrix_determinant(self): + """Test the matrix_determinant function.""" + test_ops = [ + np.eye(3), + np.eye(10, dtype=DTYPE), + atmos_flux_inversion.linalg.DiagonalOperator(np.ones(10)), + atmos_flux_inversion.correlations. + HomogeneousIsotropicCorrelation.from_array([1, 0, 0, 0, 0]), + atmos_flux_inversion.correlations. + HomogeneousIsotropicCorrelation.from_array([1, 0, 0, 0, 0], False), + atmos_flux_inversion.correlations. + HomogeneousIsotropicCorrelation.from_array( + [[1, 0, 0, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 0]] + ), + atmos_flux_inversion.correlations. + HomogeneousIsotropicCorrelation.from_array( + [[1, 0, 0, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 0]], + False + ), + ] + + for test_op in test_ops: + with self.subTest(test_op=test_op): + result = atmos_flux_inversion.linalg.matrix_determinant( + test_op + ) + self.assertIsInstance(result, float) + self.assertEqual(result, 1.0) + + def test_diagonal_determinant(self): + """Test determinants of DiagonalOperators.""" + test_cases = ( + np.arange(10), + np.arange(1, 10), + np.ones(15), + ) + + for test_case in test_cases: + with self.subTest(test_case=test_case): + op = atmos_flux_inversion.linalg.DiagonalOperator( + test_case + ) + self.assertEqual(op.det(), np.prod(test_case)) + + def test_kronecker_determinant_simple(self): + """Test determinants of Kronecker products of simple inputs.""" + kron_classes = ( + atmos_flux_inversion.linalg.DaskKroneckerProductOperator, + atmos_flux_inversion.correlations.SchmidtKroneckerProduct, + ) + test_mats = (np.eye(3), np.eye(10, dtype=DTYPE)) + test_ops = ( + atmos_flux_inversion.linalg.DiagonalOperator(np.ones(10)), + ) + + for kron_class in kron_classes: + for left, right in itertools.product( + test_mats, + test_mats + test_ops + ): + with self.subTest( + kron_class=kron_class, + left=left, right=right + ): + kron_op = kron_class(left, right) + result = kron_op.det() + self.assertEqual(result, 1) + + def test_kronecker_determinant_hard(self): + """Test determinants of KroneckerProducts of other inputs.""" + kron_classes = ( + atmos_flux_inversion.linalg.DaskKroneckerProductOperator, + atmos_flux_inversion.correlations.SchmidtKroneckerProduct, + ) + test_mats = ( + np.eye(3), + np.eye(10, dtype=DTYPE), + np.arange(4).reshape(2, 2), + np.random.rand(4, 4), + np.random.rand(5, 5), + ) + + for kron_class in kron_classes: + for left, right in itertools.product( + test_mats, + test_mats + ): + with self.subTest( + kron_class=kron_class, + left=left, right=right + ): + kron_op = kron_class(left, right) + result = kron_op.det() + expected = la.det(np.kron(left, right)) + self.assertAlmostEqual(result, expected) + + def test_cyclic_fourier_determinants(self): + """Test determinants of periodic HomogeneousIsotropicCorrelation.""" + from_function = (atmos_flux_inversion.correlations. + HomogeneousIsotropicCorrelation.from_function) + test_dists = (.3, 1, 3,) + test_shapes = (10, 15, (4, 5)) + + for test_dist, corr_class, shape in itertools.product( + test_dists, + atmos_flux_inversion.correlations. + DistanceCorrelationFunction.__subclasses__(), + test_shapes + ): + with self.subTest(test_dist=test_dist, + corr_class=corr_class.__name__, + shape=shape): + op = from_function(corr_class(test_dist), shape, + is_cyclic=True) + mat = op.dot(np.eye(*op.shape)) + + self.assertAlmostEqual(op.det(), la.det(mat), 6) + + def test_acyclic_fourier_determinants(self): + """Test determinants of aperiodic HomogeneousIsotropicCorrelation.""" + from_function = (atmos_flux_inversion.correlations. + HomogeneousIsotropicCorrelation.from_function) + test_dists = (.3, 1, 3,) + test_shapes = (10, 15, (4, 5)) + + for test_dist, corr_class, shape in itertools.product( + test_dists, + atmos_flux_inversion.correlations. + DistanceCorrelationFunction.__subclasses__(), + test_shapes + ): + with self.subTest(test_dist=test_dist, + corr_class=corr_class.__name__, + shape=shape): + op = from_function(corr_class(test_dist), shape, + is_cyclic=False) + mat = op.dot(np.eye(*op.shape)) + + self.assertAlmostEqual(op.det(), la.det(mat), 4) + + if __name__ == "__main__": unittest2.main()