Merge branch 'develop' into documentation

.github/workflows/pytest.yml

@@ -29,6 +29,7 @@ jobs:
       - name: Install
         run: |
           python -m pip install --upgrade pip
+          pip install wheel
           pip install .
           pip install pytest
           pip install pytest-cov

pyerrors/linalg.py

@@ -2,9 +2,6 @@ import numpy as np
 import autograd.numpy as anp  # Thinly-wrapped numpy
 from .obs import derived_observable, CObs, Obs, import_jackknife
 
-from functools import partial
-from autograd.extend import defvjp
-
 
 def matmul(*operands):
     """Matrix multiply all operands.
@@ -527,51 +524,3 @@ def _num_diff_svd(obs, **kwargs):
         res_mat2.append(row)
 
     return (np.array(res_mat0) @ np.identity(mid_index), np.array(res_mat1) @ np.identity(mid_index), np.array(res_mat2) @ np.identity(shape[1]))
-
-
-# This code block is directly taken from the current master branch of autograd and remains
-# only until the new version is released on PyPi
-_dot = partial(anp.einsum, '...ij,...jk->...ik')
-
-
-# batched diag
-def _diag(a):
-    return anp.eye(a.shape[-1]) * a
-
-
-# batched diagonal, similar to matrix_diag in tensorflow
-def _matrix_diag(a):
-    reps = anp.array(a.shape)
-    reps[:-1] = 1
-    reps[-1] = a.shape[-1]
-    newshape = list(a.shape) + [a.shape[-1]]
-    return _diag(anp.tile(a, reps).reshape(newshape))
-
-# https://arxiv.org/pdf/1701.00392.pdf Eq(4.77)
-# Note the formula from Sec3.1 in https://people.maths.ox.ac.uk/gilesm/files/NA-08-01.pdf is incomplete
-
-
-def grad_eig(ans, x):
-    """Gradient of a general square (complex valued) matrix"""
-    e, u = ans  # eigenvalues as 1d array, eigenvectors in columns
-    n = e.shape[-1]
-
-    def vjp(g):
-        ge, gu = g
-        ge = _matrix_diag(ge)
-        f = 1 / (e[..., anp.newaxis, :] - e[..., :, anp.newaxis] + 1.e-20)
-        f -= _diag(f)
-        ut = anp.swapaxes(u, -1, -2)
-        r1 = f * _dot(ut, gu)
-        r2 = -f * (_dot(_dot(ut, anp.conj(u)), anp.real(_dot(ut, gu)) * anp.eye(n)))
-        r = _dot(_dot(anp.linalg.inv(ut), ge + r1 + r2), ut)
-        if not anp.iscomplexobj(x):
-            r = anp.real(r)
-            # the derivative is still complex for real input (imaginary delta is allowed), real output
-            # but the derivative should be real in real input case when imaginary delta is forbidden
-        return r
-    return vjp
-
-
-defvjp(anp.linalg.eig, grad_eig)
-# End of the code block from autograd.master

setup.py

@@ -9,5 +9,5 @@ setup(name='pyerrors',
       author_email='fabian.joswig@ed.ac.uk',
       packages=find_packages(),
       python_requires='>=3.6.0',
-      install_requires=['numpy>=1.16', 'autograd>=1.2', 'numdifftools', 'matplotlib>=3.3', 'scipy', 'iminuit>=2', 'h5py']
+      install_requires=['numpy>=1.16', 'autograd @ git+https://github.com/HIPS/autograd.git', 'numdifftools', 'matplotlib>=3.3', 'scipy', 'iminuit>=2', 'h5py']
      )

tests/io_test.py

@@ -1,7 +1,8 @@
+import os
+import gzip
+import numpy as np
 import pyerrors as pe
 import pyerrors.input.json as jsonio
-import numpy as np
-import os
 
 
 def test_jsonio():
@@ -70,6 +71,15 @@ def test_jsonio():
 
 def test_json_string_reconstruction():
     my_obs = pe.Obs([np.random.rand(100)], ['name'])
+
     json_string = pe.input.json.create_json_string(my_obs)
-    reconstructed_obs = pe.input.json.import_json_string(json_string)
+    reconstructed_obs1 = pe.input.json.import_json_string(json_string)
-    assert my_obs == reconstructed_obs
+    assert my_obs == reconstructed_obs1
+
+    compressed_string = gzip.compress(json_string.encode('utf-8'))
+
+    reconstructed_string = gzip.decompress(compressed_string).decode('utf-8')
+    reconstructed_obs2 = pe.input.json.import_json_string(reconstructed_string)
+
+    assert reconstructed_string == json_string
+    assert my_obs == reconstructed_obs2

tests/linalg_test.py

@@ -300,6 +300,10 @@ def test_matrix_functions():
         for j in range(dim):
             assert tmp[j].is_zero()
 
+    # Check eig function
+    e2 = pe.linalg.eig(sym)
+    assert np.all(np.sort(e) == np.sort(e2))
+
     # Check svd
     u, v, vh = pe.linalg.svd(sym)
     diff = sym - u @ np.diag(v) @ vh