From c6f0d7e3dd9b9cb53c8e7b7aeef91838d9f9ef5d Mon Sep 17 00:00:00 2001
From: DWesl <22566757+DWesl@users.noreply.github.com>
Date: Sat, 14 Sep 2019 10:06:34 -0400
Subject: [PATCH 01/12] Add and test determinant methods for the easy classes.

---
 src/atmos_flux_inversion/linalg.py | 48 ++++++++++++++++++++++-
 src/atmos_flux_inversion/tests.py  | 61 ++++++++++++++++++++++++++++++
 2 files changed, 107 insertions(+), 2 deletions(-)

diff --git a/src/atmos_flux_inversion/linalg.py b/src/atmos_flux_inversion/linalg.py
index ba9e2c8..b2a68d5 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,18 @@ def quadratic_form(self, mat):
                 operator2.dot(chunk))
         return result
 
+    def det(self):
+        """Find the determinant of the operator.
+
+        Returns
+        -------
+        float
+        """
+        return (
+            matrix_determinant(self._operator1) *
+            matrix_determinant(self._operator2)
+        )
+
 
 class SchmidtKroneckerProduct(DaskLinearOperator):
     """Kronecker product of two operators using Schmidt decomposition.
@@ -705,5 +735,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 9e61bf5..84c2ee5 100644
--- a/src/atmos_flux_inversion/tests.py
+++ b/src/atmos_flux_inversion/tests.py
@@ -3142,5 +3142,66 @@ 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)),
+        ]
+
+        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(self):
+        """Test determinants of Kronecker products."""
+
+        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)
+
+
 if __name__ == "__main__":
     unittest2.main()

From da49d949fbfc9c22a438abf89a82ae0fb10baa1d Mon Sep 17 00:00:00 2001
From: DWesl <22566757+DWesl@users.noreply.github.com>
Date: Sat, 14 Sep 2019 10:26:37 -0400
Subject: [PATCH 02/12] Test determinants of non-simple operators and fix
 functions.

---
 src/atmos_flux_inversion/linalg.py | 30 ++++++++++++++++++++++++++--
 src/atmos_flux_inversion/tests.py  | 32 ++++++++++++++++++++++++++++--
 2 files changed, 58 insertions(+), 4 deletions(-)

diff --git a/src/atmos_flux_inversion/linalg.py b/src/atmos_flux_inversion/linalg.py
index b2a68d5..d7430bf 100644
--- a/src/atmos_flux_inversion/linalg.py
+++ b/src/atmos_flux_inversion/linalg.py
@@ -516,9 +516,16 @@ def det(self):
         -------
         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(self._operator1) *
-            matrix_determinant(self._operator2)
+            matrix_determinant(op1) ** op2.shape[0] *
+            matrix_determinant(op2) ** op1.shape[0]
         )
 
 
@@ -610,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.
diff --git a/src/atmos_flux_inversion/tests.py b/src/atmos_flux_inversion/tests.py
index 84c2ee5..9b0c834 100644
--- a/src/atmos_flux_inversion/tests.py
+++ b/src/atmos_flux_inversion/tests.py
@@ -3177,8 +3177,8 @@ def test_diagonal_determinant(self):
                 )
                 self.assertEqual(op.det(), np.prod(test_case))
 
-    def test_kronecker_determinant(self):
-        """Test determinants of Kronecker products."""
+    def test_kronecker_determinant_simple(self):
+        """Test determinants of Kronecker products of simple inputs."""
 
         kron_classes = (
             atmos_flux_inversion.linalg.DaskKroneckerProductOperator,
@@ -3202,6 +3202,34 @@ def test_kronecker_determinant(self):
                     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)
+
 
 if __name__ == "__main__":
     unittest2.main()

From 268caf7771d3443a8c8bccb3433c3454161396a9 Mon Sep 17 00:00:00 2001
From: DWesl <22566757+DWesl@users.noreply.github.com>
Date: Sat, 14 Sep 2019 14:49:23 -0400
Subject: [PATCH 03/12] Fix a bug in
 HomogeneousIsotropicCorrelation.from_array(is_cyclic=False).

It was dropping the first and last elements along the first axis every
time, instead of along the axis being concatenated.
---
 src/atmos_flux_inversion/correlations.py |  6 +++++-
 src/atmos_flux_inversion/tests.py        | 10 ++++++++++
 2 files changed, 15 insertions(+), 1 deletion(-)

diff --git a/src/atmos_flux_inversion/correlations.py b/src/atmos_flux_inversion/correlations.py
index d5ff55a..59ff61b 100644
--- a/src/atmos_flux_inversion/correlations.py
+++ b/src/atmos_flux_inversion/correlations.py
@@ -284,8 +284,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
diff --git a/src/atmos_flux_inversion/tests.py b/src/atmos_flux_inversion/tests.py
index 9b0c834..b624ea8 100644
--- a/src/atmos_flux_inversion/tests.py
+++ b/src/atmos_flux_inversion/tests.py
@@ -1201,6 +1201,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.

From 5644bba5ada2deb88f7943e007fb259fd4ffd754 Mon Sep 17 00:00:00 2001
From: DWesl <22566757+DWesl@users.noreply.github.com>
Date: Sat, 14 Sep 2019 14:51:28 -0400
Subject: [PATCH 04/12] WIP: Try to find determinant of
 HomogeneousIsotropicCorrelations.

The primary data stored by this class is a list of eigenvalues, so
this should be easy.  For some reason it's broken for all but trivial
cases.  I can get the determinant to be one for diagonal correlations,
but this is not the primary class for diagonal operators, especially
not ones proportional to the identity matrix.  I have no idea why this
doesn't work.  This may be enough for the purposes of this package,
since it's usually only off by 50% or so, but I really don't like
calling the result a `determinant`.
---
 src/atmos_flux_inversion/correlations.py | 17 +++++++++-
 src/atmos_flux_inversion/tests.py        | 42 ++++++++++++++++++++++++
 2 files changed, 58 insertions(+), 1 deletion(-)

diff --git a/src/atmos_flux_inversion/correlations.py b/src/atmos_flux_inversion/correlations.py
index 59ff61b..ad1f017 100644
--- a/src/atmos_flux_inversion/correlations.py
+++ b/src/atmos_flux_inversion/correlations.py
@@ -17,7 +17,7 @@
 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
@@ -443,6 +443,21 @@ def kron(self, other):
             other._fourier_near_zero[other_index])
         return newinst
 
+    def det(self):
+        """Find the determinant of the operator.
+
+        Returns
+        -------
+        float
+        """
+        if self._is_cyclic:
+            return prod(self._corr_fourier.real)
+        index = tuple(
+            slice(1, None, 2)
+            for _ in self._underlying_shape
+        )
+        return prod(self._corr_fourier[index].real)
+
 
 def make_matrix(corr_func, shape):
     """Make a correlation matrix for a domain with shape `shape`.
diff --git a/src/atmos_flux_inversion/tests.py b/src/atmos_flux_inversion/tests.py
index b624ea8..68d5751 100644
--- a/src/atmos_flux_inversion/tests.py
+++ b/src/atmos_flux_inversion/tests.py
@@ -3162,6 +3162,25 @@ def test_matrix_determinant(self):
             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:
@@ -3240,6 +3259,29 @@ def test_kronecker_determinant_hard(self):
                     expected = la.det(np.kron(left, right))
                     self.assertAlmostEqual(result, expected)
 
+    def test_fourier_determinants(self):
+        """Test determinants of 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, is_cyclic in itertools.product(
+                test_dists,
+                atmos_flux_inversion.correlations.
+                DistanceCorrelationFunction.__subclasses__(),
+                test_shapes,
+                (True, False)
+        ):
+            with self.subTest(test_dist=test_dist,
+                              corr_class=corr_class.__name__,
+                              shape=shape,
+                              is_cyclic=is_cyclic):
+                op = from_function(corr_class(test_dist), shape, is_cyclic)
+                mat = op.dot(np.eye(*op.shape))
+
+                self.assertAlmostEqual(op.det(), la.det(mat), 4)
+
 
 if __name__ == "__main__":
     unittest2.main()

From da13bd81e186e1bddf2636e5e49a308518f5326d Mon Sep 17 00:00:00 2001
From: DWesl <22566757+DWesl@users.noreply.github.com>
Date: Thu, 17 Oct 2019 18:34:48 -0400
Subject: [PATCH 05/12] Fix determinant for cyclic Homogeneous correlations.

Non-cyclic homogeneous correlations are still broken.
I still can't figure out how to get those working.
---
 src/atmos_flux_inversion/correlations.py | 14 ++++++++------
 1 file changed, 8 insertions(+), 6 deletions(-)

diff --git a/src/atmos_flux_inversion/correlations.py b/src/atmos_flux_inversion/correlations.py
index ad1f017..cef89c5 100644
--- a/src/atmos_flux_inversion/correlations.py
+++ b/src/atmos_flux_inversion/correlations.py
@@ -21,6 +21,7 @@
 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
@@ -450,13 +451,14 @@ def det(self):
         -------
         float
         """
+        correlations = self._ifft(self._corr_fourier)
         if self._is_cyclic:
-            return prod(self._corr_fourier.real)
-        index = tuple(
-            slice(1, None, 2)
-            for _ in self._underlying_shape
-        )
-        return prod(self._corr_fourier[index].real)
+            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):

From 15900c55e9c6270aff4db4b463f6f6e0c491d409 Mon Sep 17 00:00:00 2001
From: DWesl <22566757+DWesl@users.noreply.github.com>
Date: Sat, 14 Sep 2019 10:06:34 -0400
Subject: [PATCH 06/12] Add and test determinant methods for the easy classes.

---
 src/atmos_flux_inversion/linalg.py | 48 ++++++++++++++++++++++-
 src/atmos_flux_inversion/tests.py  | 61 ++++++++++++++++++++++++++++++
 2 files changed, 107 insertions(+), 2 deletions(-)

diff --git a/src/atmos_flux_inversion/linalg.py b/src/atmos_flux_inversion/linalg.py
index ba9e2c8..b2a68d5 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,18 @@ def quadratic_form(self, mat):
                 operator2.dot(chunk))
         return result
 
+    def det(self):
+        """Find the determinant of the operator.
+
+        Returns
+        -------
+        float
+        """
+        return (
+            matrix_determinant(self._operator1) *
+            matrix_determinant(self._operator2)
+        )
+
 
 class SchmidtKroneckerProduct(DaskLinearOperator):
     """Kronecker product of two operators using Schmidt decomposition.
@@ -705,5 +735,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 9e61bf5..84c2ee5 100644
--- a/src/atmos_flux_inversion/tests.py
+++ b/src/atmos_flux_inversion/tests.py
@@ -3142,5 +3142,66 @@ 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)),
+        ]
+
+        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(self):
+        """Test determinants of Kronecker products."""
+
+        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)
+
+
 if __name__ == "__main__":
     unittest2.main()

From 5e10a2bd1f42b7ae8de7b189540f2e6ec04a6242 Mon Sep 17 00:00:00 2001
From: DWesl <22566757+DWesl@users.noreply.github.com>
Date: Sat, 14 Sep 2019 10:26:37 -0400
Subject: [PATCH 07/12] Test determinants of non-simple operators and fix
 functions.

---
 src/atmos_flux_inversion/linalg.py | 30 ++++++++++++++++++++++++++--
 src/atmos_flux_inversion/tests.py  | 32 ++++++++++++++++++++++++++++--
 2 files changed, 58 insertions(+), 4 deletions(-)

diff --git a/src/atmos_flux_inversion/linalg.py b/src/atmos_flux_inversion/linalg.py
index b2a68d5..d7430bf 100644
--- a/src/atmos_flux_inversion/linalg.py
+++ b/src/atmos_flux_inversion/linalg.py
@@ -516,9 +516,16 @@ def det(self):
         -------
         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(self._operator1) *
-            matrix_determinant(self._operator2)
+            matrix_determinant(op1) ** op2.shape[0] *
+            matrix_determinant(op2) ** op1.shape[0]
         )
 
 
@@ -610,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.
diff --git a/src/atmos_flux_inversion/tests.py b/src/atmos_flux_inversion/tests.py
index 84c2ee5..9b0c834 100644
--- a/src/atmos_flux_inversion/tests.py
+++ b/src/atmos_flux_inversion/tests.py
@@ -3177,8 +3177,8 @@ def test_diagonal_determinant(self):
                 )
                 self.assertEqual(op.det(), np.prod(test_case))
 
-    def test_kronecker_determinant(self):
-        """Test determinants of Kronecker products."""
+    def test_kronecker_determinant_simple(self):
+        """Test determinants of Kronecker products of simple inputs."""
 
         kron_classes = (
             atmos_flux_inversion.linalg.DaskKroneckerProductOperator,
@@ -3202,6 +3202,34 @@ def test_kronecker_determinant(self):
                     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)
+
 
 if __name__ == "__main__":
     unittest2.main()

From 981a831f7bfcf6ce474c0f23117d3bee7d18a3bd Mon Sep 17 00:00:00 2001
From: DWesl <22566757+DWesl@users.noreply.github.com>
Date: Sat, 14 Sep 2019 14:49:23 -0400
Subject: [PATCH 08/12] Fix a bug in
 HomogeneousIsotropicCorrelation.from_array(is_cyclic=False).

It was dropping the first and last elements along the first axis every
time, instead of along the axis being concatenated.
---
 src/atmos_flux_inversion/correlations.py |  6 +++++-
 src/atmos_flux_inversion/tests.py        | 10 ++++++++++
 2 files changed, 15 insertions(+), 1 deletion(-)

diff --git a/src/atmos_flux_inversion/correlations.py b/src/atmos_flux_inversion/correlations.py
index d5ff55a..59ff61b 100644
--- a/src/atmos_flux_inversion/correlations.py
+++ b/src/atmos_flux_inversion/correlations.py
@@ -284,8 +284,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
diff --git a/src/atmos_flux_inversion/tests.py b/src/atmos_flux_inversion/tests.py
index 9b0c834..b624ea8 100644
--- a/src/atmos_flux_inversion/tests.py
+++ b/src/atmos_flux_inversion/tests.py
@@ -1201,6 +1201,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.

From 81f5b84b69e4f4ef72ba26cf562828b7535d5e0f Mon Sep 17 00:00:00 2001
From: DWesl <22566757+DWesl@users.noreply.github.com>
Date: Sat, 14 Sep 2019 14:51:28 -0400
Subject: [PATCH 09/12] WIP: Try to find determinant of
 HomogeneousIsotropicCorrelations.

The primary data stored by this class is a list of eigenvalues, so
this should be easy.  For some reason it's broken for all but trivial
cases.  I can get the determinant to be one for diagonal correlations,
but this is not the primary class for diagonal operators, especially
not ones proportional to the identity matrix.  I have no idea why this
doesn't work.  This may be enough for the purposes of this package,
since it's usually only off by 50% or so, but I really don't like
calling the result a `determinant`.
---
 src/atmos_flux_inversion/correlations.py | 17 +++++++++-
 src/atmos_flux_inversion/tests.py        | 42 ++++++++++++++++++++++++
 2 files changed, 58 insertions(+), 1 deletion(-)

diff --git a/src/atmos_flux_inversion/correlations.py b/src/atmos_flux_inversion/correlations.py
index 59ff61b..ad1f017 100644
--- a/src/atmos_flux_inversion/correlations.py
+++ b/src/atmos_flux_inversion/correlations.py
@@ -17,7 +17,7 @@
 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
@@ -443,6 +443,21 @@ def kron(self, other):
             other._fourier_near_zero[other_index])
         return newinst
 
+    def det(self):
+        """Find the determinant of the operator.
+
+        Returns
+        -------
+        float
+        """
+        if self._is_cyclic:
+            return prod(self._corr_fourier.real)
+        index = tuple(
+            slice(1, None, 2)
+            for _ in self._underlying_shape
+        )
+        return prod(self._corr_fourier[index].real)
+
 
 def make_matrix(corr_func, shape):
     """Make a correlation matrix for a domain with shape `shape`.
diff --git a/src/atmos_flux_inversion/tests.py b/src/atmos_flux_inversion/tests.py
index b624ea8..68d5751 100644
--- a/src/atmos_flux_inversion/tests.py
+++ b/src/atmos_flux_inversion/tests.py
@@ -3162,6 +3162,25 @@ def test_matrix_determinant(self):
             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:
@@ -3240,6 +3259,29 @@ def test_kronecker_determinant_hard(self):
                     expected = la.det(np.kron(left, right))
                     self.assertAlmostEqual(result, expected)
 
+    def test_fourier_determinants(self):
+        """Test determinants of 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, is_cyclic in itertools.product(
+                test_dists,
+                atmos_flux_inversion.correlations.
+                DistanceCorrelationFunction.__subclasses__(),
+                test_shapes,
+                (True, False)
+        ):
+            with self.subTest(test_dist=test_dist,
+                              corr_class=corr_class.__name__,
+                              shape=shape,
+                              is_cyclic=is_cyclic):
+                op = from_function(corr_class(test_dist), shape, is_cyclic)
+                mat = op.dot(np.eye(*op.shape))
+
+                self.assertAlmostEqual(op.det(), la.det(mat), 4)
+
 
 if __name__ == "__main__":
     unittest2.main()

From 0b17702c06f4913384acd92ccceea88ff2f6aea9 Mon Sep 17 00:00:00 2001
From: DWesl <22566757+DWesl@users.noreply.github.com>
Date: Thu, 17 Oct 2019 18:34:48 -0400
Subject: [PATCH 10/12] Fix determinant for cyclic Homogeneous correlations.

Non-cyclic homogeneous correlations are still broken.
I still can't figure out how to get those working.
---
 src/atmos_flux_inversion/correlations.py | 14 ++++++++------
 1 file changed, 8 insertions(+), 6 deletions(-)

diff --git a/src/atmos_flux_inversion/correlations.py b/src/atmos_flux_inversion/correlations.py
index ad1f017..cef89c5 100644
--- a/src/atmos_flux_inversion/correlations.py
+++ b/src/atmos_flux_inversion/correlations.py
@@ -21,6 +21,7 @@
 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
@@ -450,13 +451,14 @@ def det(self):
         -------
         float
         """
+        correlations = self._ifft(self._corr_fourier)
         if self._is_cyclic:
-            return prod(self._corr_fourier.real)
-        index = tuple(
-            slice(1, None, 2)
-            for _ in self._underlying_shape
-        )
-        return prod(self._corr_fourier[index].real)
+            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):

From 508e6438504b103a5ad94cf2409fa0da16a3bfd9 Mon Sep 17 00:00:00 2001
From: DWesl <22566757+DWesl@users.noreply.github.com>
Date: Fri, 8 Nov 2019 11:52:36 -0500
Subject: [PATCH 11/12] Rename test environment `dist` to `distrib`.

Gets tox to try it again.
---
 tox.ini | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/tox.ini b/tox.ini
index 936f9fb..213c73d 100644
--- a/tox.ini
+++ b/tox.ini
@@ -4,7 +4,7 @@
 # and then run "tox" from this directory.
 
 [tox]
-envlist = cov_erase, py2, py3, docs, dist, codestyle, coverage
+envlist = cov_erase, py2, py3, docs, distrib, codestyle, coverage
 
 [testenv]
 commands =
@@ -75,7 +75,7 @@ commands=
 skip_install=true
 depends=
 
-[testenv:dist]
+[testenv:distrib]
 deps=
     setuptools
     twine

From e5f5fb5be8b2f0cbb20a698884cdf6e4a149ba17 Mon Sep 17 00:00:00 2001
From: DWesl <22566757+DWesl@users.noreply.github.com>
Date: Fri, 8 Nov 2019 11:56:55 -0500
Subject: [PATCH 12/12] Split the now-working periodic determinant tests from
 the still-broken apreriodic tests.

At this point just under half of the aperiodic tests are broken.  Some
of them are close and would pass if I required only three-decimal
precision, others differ by an order of magnitude.
---
 src/atmos_flux_inversion/tests.py | 39 +++++++++++++++++++++++--------
 1 file changed, 29 insertions(+), 10 deletions(-)

diff --git a/src/atmos_flux_inversion/tests.py b/src/atmos_flux_inversion/tests.py
index 68d5751..1cd7d7d 100644
--- a/src/atmos_flux_inversion/tests.py
+++ b/src/atmos_flux_inversion/tests.py
@@ -3157,7 +3157,6 @@ class TestDeterminants(unittest2.TestCase):
 
     def test_matrix_determinant(self):
         """Test the matrix_determinant function."""
-
         test_ops = [
             np.eye(3),
             np.eye(10, dtype=DTYPE),
@@ -3208,7 +3207,6 @@ def test_diagonal_determinant(self):
 
     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,
@@ -3259,25 +3257,46 @@ def test_kronecker_determinant_hard(self):
                     expected = la.det(np.kron(left, right))
                     self.assertAlmostEqual(result, expected)
 
-    def test_fourier_determinants(self):
-        """Test determinants of HomogeneousIsotropicCorrelation."""
+    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, is_cyclic in itertools.product(
+        for test_dist, corr_class, shape in itertools.product(
                 test_dists,
                 atmos_flux_inversion.correlations.
                 DistanceCorrelationFunction.__subclasses__(),
-                test_shapes,
-                (True, False)
+                test_shapes
         ):
             with self.subTest(test_dist=test_dist,
                               corr_class=corr_class.__name__,
-                              shape=shape,
-                              is_cyclic=is_cyclic):
-                op = from_function(corr_class(test_dist), shape, is_cyclic)
+                              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)