From 751d28711f9ff70bd07494b9dd4648cce2c25103 Mon Sep 17 00:00:00 2001
From: Simon Kuberski <simon.kuberski@uni-muenster.de>
Date: Fri, 25 Feb 2022 17:39:49 +0100
Subject: [PATCH] feat: Check for symmetry and positive-semidefiniteness of
 covariance matrices in the initialization of covobs

---
 pyerrors/covobs.py   | 29 +++++++++++++++++++++++++++++
 tests/covobs_test.py |  9 ++++++---
 2 files changed, 35 insertions(+), 3 deletions(-)

diff --git a/pyerrors/covobs.py b/pyerrors/covobs.py
index 2c8f81fc..49953983 100644
--- a/pyerrors/covobs.py
+++ b/pyerrors/covobs.py
@@ -45,6 +45,16 @@ class Covobs:
         return float(np.dot(np.transpose(self.grad), np.dot(self.cov, self.grad)))
 
     def _set_cov(self, cov):
+        """ Set the covariance matrix of the covobs
+
+        Parameters
+        ----------
+        cov : list or array
+            Has to be either of:
+            0 dimensional number: variance of a single covobs,
+            1 dimensional list or array of lenght N: variances of multiple covobs
+            2 dimensional list or array (N x N): Symmetric, positive-semidefinite covariance matrix
+        """
         self._cov = np.array(cov)
         if self._cov.ndim == 0:
             self.N = 1
@@ -59,7 +69,26 @@ class Covobs:
         else:
             raise Exception('Covariance matrix has to be a 2 dimensional square matrix!')
 
+        for i in range(self.N):
+            for j in range(i):
+                if not self._cov[i][j] == self._cov[j][i]:
+                    raise Exception('Covariance matrix is non-symmetric for (%d, %d' % (i, j))
+
+        evals = np.linalg.eigvalsh(self._cov)
+        for ev in evals:
+            if ev < 0:
+                raise Exception('Covariance matrix is not positive-semidefinite!')
+
     def _set_grad(self, grad):
+        """ Set the gradient of the covobs
+
+        Parameters
+        ----------
+        grad : list or array
+            Has to be either of:
+            0 dimensional number: gradient w.r.t. a single covobs,
+            1 dimensional list or array of lenght N: gradient w.r.t. multiple covobs
+        """
         self._grad = np.array(grad)
         if self._grad.ndim in [0, 1]:
             self._grad = np.reshape(self._grad, (self.N, 1))
diff --git a/tests/covobs_test.py b/tests/covobs_test.py
index 215c75d3..418cb83e 100644
--- a/tests/covobs_test.py
+++ b/tests/covobs_test.py
@@ -78,9 +78,8 @@ def test_covobs_init():
     covobs = pe.cov_Obs(0.5, 0.002, 'test')
     covobs = pe.cov_Obs([1, 2], [0.1, 0.2], 'test')
     covobs = pe.cov_Obs([1, 2], np.array([0.1, 0.2]), 'test')
-    covobs = pe.cov_Obs([1, 2], [[0.1, 0.2], [0.1, 0.2]], 'test')
-    covobs = pe.cov_Obs([1, 2], np.array([[0.1, 0.2], [0.1, 0.2]]), 'test')
-
+    covobs = pe.cov_Obs([1, 2], [[0.21, 0.2], [0.2, 0.21]], 'test')
+    covobs = pe.cov_Obs([1, 2], np.array([[0.21, 0.2], [0.2, 0.21]]), 'test')
 
 
 def test_covobs_exceptions():
@@ -92,3 +91,7 @@ def test_covobs_exceptions():
         covobs = pe.cov_Obs([0.5, 0.1], np.array([[2, 1, 3], [1, 2, 3]]), 'test')
     with pytest.raises(Exception):
         covobs = pe.cov_Obs([0.5, 0.1], np.random.random((2, 2, 2)), 'test')
+    with pytest.raises(Exception):
+        covobs = pe.cov_Obs([1.5, 0.1], [[1., .2,], [.3, .5]] , 'test')
+    with pytest.raises(Exception):
+        covobs = pe.cov_Obs([1.5, 0.1], [[8, 4,], [4, -2]] , 'test')