From 0644ecf9aaeba0581605d6ecce038c54b135ef79 Mon Sep 17 00:00:00 2001
From: Fabian Joswig <fabian.joswig@ed.ac.uk>
Date: Fri, 12 Nov 2021 09:39:18 +0000
Subject: [PATCH 01/27] bugs in derived_array and gamma_method fixed which
 caused odd behaviour in connection with irregular monte carlo chains

---
 pyerrors/linalg.py | 2 +-
 pyerrors/obs.py    | 4 ++--
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/pyerrors/linalg.py b/pyerrors/linalg.py
index 18a6a98c..c1a042fb 100644
--- a/pyerrors/linalg.py
+++ b/pyerrors/linalg.py
@@ -86,7 +86,7 @@ def derived_array(func, data, **kwargs):
     for name in new_names:
         d_extracted[name] = []
         for i_dat, dat in enumerate(data):
-            ens_length = new_idl_d[name][-1] - new_idl_d[name][0] + 1
+            ens_length = len(new_idl_d[name])
             d_extracted[name].append(np.array([_expand_deltas_for_merge(o.deltas[name], o.idl[name], o.shape[name], new_idl_d[name]) for o in dat.reshape(np.prod(dat.shape))]).reshape(dat.shape + (ens_length, )))
 
     for i_val, new_val in np.ndenumerate(new_values):
diff --git a/pyerrors/obs.py b/pyerrors/obs.py
index 6a2b414a..085c9a95 100644
--- a/pyerrors/obs.py
+++ b/pyerrors/obs.py
@@ -248,7 +248,7 @@ class Obs:
 
             r_length = []
             for r_name in e_content[e_name]:
-                if self.idl[r_name] is range:
+                if isinstance(self.idl[r_name], range):
                     r_length.append(len(self.idl[r_name]))
                 else:
                     r_length.append((self.idl[r_name][-1] - self.idl[r_name][0] + 1))
@@ -339,7 +339,7 @@ class Obs:
         idx     -- List or range of configs on which the deltas are defined.
         shape   -- Number of configs in idx.
         """
-        if type(idx) is range:
+        if isinstance(idx, range):
             return deltas
         else:
             ret = np.zeros(idx[-1] - idx[0] + 1)

From a8a52beadd9b53bbb7fad49d7f9e9817fd1339a5 Mon Sep 17 00:00:00 2001
From: Fabian Joswig <fabian.joswig@ed.ac.uk>
Date: Fri, 12 Nov 2021 09:41:17 +0000
Subject: [PATCH 02/27] same bug fixed in covariance2

---
 pyerrors/obs.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyerrors/obs.py b/pyerrors/obs.py
index 085c9a95..7e98ae8c 100644
--- a/pyerrors/obs.py
+++ b/pyerrors/obs.py
@@ -1307,7 +1307,7 @@ def covariance2(obs1, obs2, correlation=False, **kwargs):
             if r_name not in obs2.e_content[e_name]:
                 continue
             idl_d[r_name] = _merge_idx([obs1.idl[r_name], obs2.idl[r_name]])
-            if idl_d[r_name] is range:
+            if isinstance(idl_d[r_name], range):
                 r_length.append(len(idl_d[r_name]))
             else:
                 r_length.append((idl_d[r_name][-1] - idl_d[r_name][0] + 1))

From b80b746fae7d87b1a67ea6fa79230644e460147d Mon Sep 17 00:00:00 2001
From: Fabian Joswig <fabian.joswig@ed.ac.uk>
Date: Fri, 12 Nov 2021 09:56:07 +0000
Subject: [PATCH 03/27] test for irregular matrix inverse added

---
 tests/linalg_test.py | 20 ++++++++++++++++++++
 1 file changed, 20 insertions(+)

diff --git a/tests/linalg_test.py b/tests/linalg_test.py
index f2bd3e60..f4b32231 100644
--- a/tests/linalg_test.py
+++ b/tests/linalg_test.py
@@ -74,6 +74,26 @@ def test_matmul_irregular_histories():
         assert np.all([o.is_merged for o in t2.ravel()])
 
 
+def test_irregular_matrix_inverse():
+    dim = 3
+    length = 500
+
+    for idl in [range(8, 508, 10), range(250, 273), [2, 8, 19, 20, 78, 99, 828, 10548979]]:
+        irregular_array = []
+        for i in range(dim ** 2):
+            irregular_array.append(pe.Obs([np.random.normal(1.1, 0.2, len(idl)), np.random.normal(0.25, 0.1, 10)], ['ens1', 'ens2'], idl=[idl, range(1, 11)]))
+        irregular_matrix = np.array(irregular_array).reshape((dim, dim))
+
+        invertible_irregular_matrix = np.identity(dim) + irregular_matrix @ irregular_matrix.T
+
+        inverse = pe.linalg.inv(invertible_irregular_matrix)
+
+        assert np.allclose(np.linalg.inv(np.vectorize(lambda x: x.value)(invertible_irregular_matrix)) - np.vectorize(lambda x: x.value)(inverse), 0.0)
+
+        check = pe.linalg.matmul(invertible_irregular_matrix, inverse)
+        assert np.all([o.is_zero() for o in (check - np.identity(dim)).ravel()])
+
+
 def test_matrix_inverse():
     content = []
     for t in range(9):

From d26bbdc3bd72b48985956beeded61ad65c756379 Mon Sep 17 00:00:00 2001
From: Fabian Joswig <fabian.joswig@ed.ac.uk>
Date: Fri, 12 Nov 2021 11:41:23 +0000
Subject: [PATCH 04/27] expand_deltas renamed to _expand_deltas and made a
 function instead of class method

---
 pyerrors/obs.py      | 44 ++++++++++++++++++++++++--------------------
 tests/linalg_test.py |  6 ++++--
 2 files changed, 28 insertions(+), 22 deletions(-)

diff --git a/pyerrors/obs.py b/pyerrors/obs.py
index 7e98ae8c..bcf542f6 100644
--- a/pyerrors/obs.py
+++ b/pyerrors/obs.py
@@ -329,24 +329,6 @@ class Obs:
             self.ddvalue = np.sqrt(self.ddvalue) / self.dvalue
         return
 
-    def expand_deltas(self, deltas, idx, shape):
-        """Expand deltas defined on idx to a regular, contiguous range, where holes are filled by 0.
-           If idx is of type range, the deltas are not changed
-
-        Parameters
-        ----------
-        deltas  -- List of fluctuations
-        idx     -- List or range of configs on which the deltas are defined.
-        shape   -- Number of configs in idx.
-        """
-        if isinstance(idx, range):
-            return deltas
-        else:
-            ret = np.zeros(idx[-1] - idx[0] + 1)
-            for i in range(shape):
-                ret[idx[i] - idx[0]] = deltas[i]
-            return ret
-
     def calc_gamma(self, deltas, idx, shape, w_max, fft):
         """Calculate Gamma_{AA} from the deltas, which are defined on idx.
            idx is assumed to be a contiguous range (possibly with a stepsize != 1)
@@ -361,7 +343,7 @@ class Obs:
                    the computation of the autocorrelation function
         """
         gamma = np.zeros(w_max)
-        deltas = self.expand_deltas(deltas, idx, shape)
+        deltas = _expand_deltas(deltas, idx, shape)
         new_shape = len(deltas)
         if fft:
             max_gamma = min(new_shape, w_max)
@@ -505,7 +487,7 @@ class Obs:
             tmp = []
             for r, r_name in enumerate(self.e_content[e_name]):
                 if expand:
-                    tmp.append(self.expand_deltas(self.deltas[r_name], self.idl[r_name], self.shape[r_name]) + self.r_values[r_name])
+                    tmp.append(_expand_deltas(self.deltas[r_name], self.idl[r_name], self.shape[r_name]) + self.r_values[r_name])
                 else:
                     tmp.append(self.deltas[r_name] + self.r_values[r_name])
                 r_length.append(len(tmp[-1]))
@@ -829,6 +811,28 @@ class CObs:
         return 'CObs[' + str(self) + ']'
 
 
+def _expand_deltas(deltas, idx, shape):
+    """Expand deltas defined on idx to a regular, contiguous range, where holes are filled by 0.
+       If idx is of type range, the deltas are not changed
+
+    Parameters
+    ----------
+    deltas : list
+        List of fluctuations
+    idx : list
+        List or range of configs on which the deltas are defined.
+    shape : int
+        Number of configs in idx.
+    """
+    if isinstance(idx, range):
+        return deltas
+    else:
+        ret = np.zeros(idx[-1] - idx[0] + 1)
+        for i in range(shape):
+            ret[idx[i] - idx[0]] = deltas[i]
+        return ret
+
+
 def _merge_idx(idl):
     """Returns the union of all lists in idl
 
diff --git a/tests/linalg_test.py b/tests/linalg_test.py
index f4b32231..2b6f200c 100644
--- a/tests/linalg_test.py
+++ b/tests/linalg_test.py
@@ -90,8 +90,10 @@ def test_irregular_matrix_inverse():
 
         assert np.allclose(np.linalg.inv(np.vectorize(lambda x: x.value)(invertible_irregular_matrix)) - np.vectorize(lambda x: x.value)(inverse), 0.0)
 
-        check = pe.linalg.matmul(invertible_irregular_matrix, inverse)
-        assert np.all([o.is_zero() for o in (check - np.identity(dim)).ravel()])
+        check1 = pe.linalg.matmul(invertible_irregular_matrix, inverse)
+        assert np.all([o.is_zero() for o in (check1 - np.identity(dim)).ravel()])
+        check2 = invertible_irregular_matrix @ inverse
+        assert np.all([o.is_zero() for o in (check2 - np.identity(dim)).ravel()])
 
 
 def test_matrix_inverse():

From 119ddba5a86a1cc29290fbdcf3d525f860f2bf50 Mon Sep 17 00:00:00 2001
From: Fabian Joswig <fabian.joswig@ed.ac.uk>
Date: Fri, 12 Nov 2021 11:50:14 +0000
Subject: [PATCH 05/27] removed calc_gamma from global namespace

---
 pyerrors/obs.py | 23 ++++++++++++++---------
 1 file changed, 14 insertions(+), 9 deletions(-)

diff --git a/pyerrors/obs.py b/pyerrors/obs.py
index bcf542f6..47e80504 100644
--- a/pyerrors/obs.py
+++ b/pyerrors/obs.py
@@ -260,11 +260,11 @@ class Obs:
             self.e_drho[e_name] = np.zeros(w_max)
 
             for r_name in e_content[e_name]:
-                e_gamma[e_name] += self.calc_gamma(self.deltas[r_name], self.idl[r_name], self.shape[r_name], w_max, fft)
+                e_gamma[e_name] += self._calc_gamma(self.deltas[r_name], self.idl[r_name], self.shape[r_name], w_max, fft)
 
             gamma_div = np.zeros(w_max)
             for r_name in e_content[e_name]:
-                gamma_div += self.calc_gamma(np.ones((self.shape[r_name])), self.idl[r_name], self.shape[r_name], w_max, fft)
+                gamma_div += self._calc_gamma(np.ones((self.shape[r_name])), self.idl[r_name], self.shape[r_name], w_max, fft)
             e_gamma[e_name] /= gamma_div[:w_max]
 
             if np.abs(e_gamma[e_name][0]) < 10 * np.finfo(float).tiny:  # Prevent division by zero
@@ -329,18 +329,23 @@ class Obs:
             self.ddvalue = np.sqrt(self.ddvalue) / self.dvalue
         return
 
-    def calc_gamma(self, deltas, idx, shape, w_max, fft):
+    def _calc_gamma(self, deltas, idx, shape, w_max, fft):
         """Calculate Gamma_{AA} from the deltas, which are defined on idx.
            idx is assumed to be a contiguous range (possibly with a stepsize != 1)
 
         Parameters
         ----------
-        deltas  -- List of fluctuations
-        idx     -- List or range of configs on which the deltas are defined.
-        shape   -- Number of configs in idx.
-        w_max   -- Upper bound for the summation window
-        fft     -- boolean, which determines whether the fft algorithm is used for
-                   the computation of the autocorrelation function
+        deltas : list
+            List of fluctuations
+        idx : list
+            List or range of configs on which the deltas are defined.
+        shape : int
+            Number of configs in idx.
+        w_max : int
+            Upper bound for the summation window.
+        fft : bool
+            determines whether the fft algorithm is used for the computation
+            of the autocorrelation function.
         """
         gamma = np.zeros(w_max)
         deltas = _expand_deltas(deltas, idx, shape)

From a628df7e5710ef028c1e9d806d4aa80945b3b0d1 Mon Sep 17 00:00:00 2001
From: Fabian Joswig <fabian.joswig@ed.ac.uk>
Date: Fri, 12 Nov 2021 11:59:40 +0000
Subject: [PATCH 06/27] docstring extended

---
 pyerrors/obs.py | 54 ++++++++++++++++++++++++++++++++++++++++++++-----
 1 file changed, 49 insertions(+), 5 deletions(-)

diff --git a/pyerrors/obs.py b/pyerrors/obs.py
index 47e80504..59ffab69 100644
--- a/pyerrors/obs.py
+++ b/pyerrors/obs.py
@@ -367,7 +367,13 @@ class Obs:
         self.details(level > 1)
 
     def details(self, ens_content=True):
-        """Output detailed properties of the Obs."""
+        """Output detailed properties of the Obs.
+
+        Parameters
+        ----------
+        ens_content : bool
+            print details about the ensembles and replica if true.
+        """
         if self.value == 0.0:
             percentage = np.nan
         else:
@@ -396,6 +402,11 @@ class Obs:
     def is_zero_within_error(self, sigma=1):
         """Checks whether the observable is zero within 'sigma' standard errors.
 
+        Parameters
+        ----------
+        sigma : int
+            Number of standard errors used for the check.
+
         Works only properly when the gamma method was run.
         """
         return self.is_zero() or np.abs(self.value) <= sigma * self.dvalue
@@ -405,7 +416,13 @@ class Obs:
         return np.isclose(0.0, self.value) and all(np.allclose(0.0, delta) for delta in self.deltas.values())
 
     def plot_tauint(self, save=None):
-        """Plot integrated autocorrelation time for each ensemble."""
+        """Plot integrated autocorrelation time for each ensemble.
+
+        Parameters
+        ----------
+        save : str
+            saves the figure to a file named 'save' if.
+        """
         if not hasattr(self, 'e_names'):
             raise Exception('Run the gamma method first.')
 
@@ -482,7 +499,13 @@ class Obs:
             plt.draw()
 
     def plot_history(self, expand=True):
-        """Plot derived Monte Carlo history for each ensemble."""
+        """Plot derived Monte Carlo history for each ensemble
+
+        Parameters
+        ----------
+        expand : bool
+            show expanded history for irregular Monte Carlo chains (default: True).
+        """
         if not hasattr(self, 'e_names'):
             raise Exception('Run the gamma method first.')
 
@@ -525,6 +548,8 @@ class Obs:
 
         Parameters
         ----------
+        name : str
+            name of the file to be saved.
         path : str
             specifies a custom path for the file (default '.')
         """
@@ -1201,6 +1226,10 @@ def covariance(obs1, obs2, correlation=False, **kwargs):
 
     Parameters
     ----------
+    obs1 : Obs
+        First Obs
+    obs2 : Obs
+        Second Obs
     correlation : bool
         if true the correlation instead of the covariance is
         returned (default False)
@@ -1434,7 +1463,16 @@ def covariance3(obs1, obs2, correlation=False, **kwargs):
 def pseudo_Obs(value, dvalue, name, samples=1000):
     """Generate a pseudo Obs with given value, dvalue and name
 
-    The standard number of samples is a 1000. This can be adjusted.
+    Parameters
+    ----------
+    value : float
+        central value of the Obs to be generated.
+    dvalue : float
+        error of the Obs to be generated.
+    name : str
+        name of the ensemble for which the Obs is to be generated.
+    samples: int
+        number of samples for the Obs (default 1000).
     """
     if dvalue <= 0.0:
         return Obs([np.zeros(samples) + value], [name])
@@ -1475,7 +1513,13 @@ def dump_object(obj, name, **kwargs):
 
 
 def load_object(path):
-    """Load object from pickle file. """
+    """Load object from pickle file.
+
+    Parameters
+    ----------
+    path : str
+        path to the file
+    """
     with open(path, 'rb') as file:
         return pickle.load(file)
 

From cf79bcc0898db72b5a3383e25bc390f4d4c26c72 Mon Sep 17 00:00:00 2001
From: Fabian Joswig <fabian.joswig@ed.ac.uk>
Date: Fri, 12 Nov 2021 13:55:03 +0000
Subject: [PATCH 07/27] parsing routine for kwargs in gamma_method added

---
 pyerrors/obs.py | 17 +++++++++++++++++
 1 file changed, 17 insertions(+)

diff --git a/pyerrors/obs.py b/pyerrors/obs.py
index 59ffab69..b0bcb7a2 100644
--- a/pyerrors/obs.py
+++ b/pyerrors/obs.py
@@ -189,6 +189,23 @@ class Obs:
         else:
             fft = True
 
+        def _parse_kwarg(kwarg_name, **kwargs):
+            if kwarg_name in kwargs:
+                tmp = kwargs.get(kwarg_name)
+                if isinstance(tmp, (int, float)):
+                    if tmp <= 0:
+                        raise Exception(kwarg_name + ' has to be larger than 0.')
+                    for e, e_name in enumerate(self.e_names):
+                        getattr(self, kwarg_name)[e_name] = tmp
+                else:
+                    raise TypeError(kwarg_name + ' is not in proper format.')
+            else:
+                for e, e_name in enumerate(self.e_names):
+                    if e_name in getattr(Obs, kwarg_name + '_dict'):
+                        getattr(self, kwarg_name)[e_name] = getattr(Obs, kwarg_name + '_dict')[e_name]
+                    else:
+                        getattr(self, kwarg_name)[e_name] = getattr(Obs, kwarg_name + '_global')
+
         if 'S' in kwargs:
             tmp = kwargs.get('S')
             if isinstance(tmp, list):

From 3d4aee703e4709591343c7e51a21f5d0458cd2da Mon Sep 17 00:00:00 2001
From: Fabian Joswig <fabian.joswig@ed.ac.uk>
Date: Fri, 12 Nov 2021 14:00:57 +0000
Subject: [PATCH 08/27] S and tau_exp kwargs in gamma method now parsed in the
 same way

---
 pyerrors/obs.py | 51 +++----------------------------------------------
 1 file changed, 3 insertions(+), 48 deletions(-)

diff --git a/pyerrors/obs.py b/pyerrors/obs.py
index b0bcb7a2..265c0649 100644
--- a/pyerrors/obs.py
+++ b/pyerrors/obs.py
@@ -189,7 +189,7 @@ class Obs:
         else:
             fft = True
 
-        def _parse_kwarg(kwarg_name, **kwargs):
+        def _parse_kwarg(kwarg_name):
             if kwarg_name in kwargs:
                 tmp = kwargs.get(kwarg_name)
                 if isinstance(tmp, (int, float)):
@@ -206,53 +206,8 @@ class Obs:
                     else:
                         getattr(self, kwarg_name)[e_name] = getattr(Obs, kwarg_name + '_global')
 
-        if 'S' in kwargs:
-            tmp = kwargs.get('S')
-            if isinstance(tmp, list):
-                if len(tmp) != len(self.e_names):
-                    raise Exception('Length of S array does not match ensembles.')
-                for e, e_name in enumerate(self.e_names):
-                    if tmp[e] <= 0:
-                        raise Exception('S has to be larger than 0.')
-                    self.S[e_name] = tmp[e]
-            else:
-                if isinstance(tmp, (int, float)):
-                    if tmp <= 0:
-                        raise Exception('S has to be larger than 0.')
-                    for e, e_name in enumerate(self.e_names):
-                        self.S[e_name] = tmp
-                else:
-                    raise TypeError('S is not in proper format.')
-        else:
-            for e, e_name in enumerate(self.e_names):
-                if e_name in Obs.S_dict:
-                    self.S[e_name] = Obs.S_dict[e_name]
-                else:
-                    self.S[e_name] = Obs.S_global
-
-        if 'tau_exp' in kwargs:
-            tmp = kwargs.get('tau_exp')
-            if isinstance(tmp, list):
-                if len(tmp) != len(self.e_names):
-                    raise Exception('Length of tau_exp array does not match ensembles.')
-                for e, e_name in enumerate(self.e_names):
-                    if tmp[e] < 0:
-                        raise Exception('tau_exp smaller than 0.')
-                    self.tau_exp[e_name] = tmp[e]
-            else:
-                if isinstance(tmp, (int, float)):
-                    if tmp < 0:
-                        raise Exception('tau_exp smaller than 0.')
-                    for e, e_name in enumerate(self.e_names):
-                        self.tau_exp[e_name] = tmp
-                else:
-                    raise TypeError('tau_exp is not in proper format.')
-        else:
-            for e, e_name in enumerate(self.e_names):
-                if e_name in Obs.tau_exp_dict:
-                    self.tau_exp[e_name] = Obs.tau_exp_dict[e_name]
-                else:
-                    self.tau_exp[e_name] = Obs.tau_exp_global
+        _parse_kwarg('S')
+        _parse_kwarg('tau_exp')
 
         if 'N_sigma' in kwargs:
             self.N_sigma = kwargs.get('N_sigma')

From 56e1425835949a1c4935e03413d0ca4ace0e563e Mon Sep 17 00:00:00 2001
From: Fabian Joswig <fabian.joswig@ed.ac.uk>
Date: Fri, 12 Nov 2021 14:06:06 +0000
Subject: [PATCH 09/27] Bug in kwarg parsing fixed, N_sigma_dict added

---
 pyerrors/obs.py | 23 +++++++++++------------
 1 file changed, 11 insertions(+), 12 deletions(-)

diff --git a/pyerrors/obs.py b/pyerrors/obs.py
index 265c0649..1c98d2a3 100644
--- a/pyerrors/obs.py
+++ b/pyerrors/obs.py
@@ -28,11 +28,14 @@ class Obs:
         exists this overwrites the standard value for that ensemble.
     tau_exp_global : float
         Standard value for tau_exp (default 0.0)
-    tau_exp_dict :dict
+    tau_exp_dict : dict
         Dictionary for tau_exp values. If an entry for a given ensemble exists
         this overwrites the standard value for that ensemble.
     N_sigma_global : float
         Standard value for N_sigma (default 1.0)
+    N_sigma_dict : dict
+        Dictionary for N_sigma values. If an entry for a given ensemble exists
+        this overwrites the standard value for that ensemble.
     """
     __slots__ = ['names', 'shape', 'r_values', 'deltas', 'N', '_value', '_dvalue',
                  'ddvalue', 'reweighted', 'S', 'tau_exp', 'N_sigma',
@@ -45,6 +48,7 @@ class Obs:
     tau_exp_global = 0.0
     tau_exp_dict = {}
     N_sigma_global = 1.0
+    N_sigma_dict = {}
     filter_eps = 1e-10
 
     def __init__(self, samples, names, idl=None, means=None, **kwargs):
@@ -183,6 +187,7 @@ class Obs:
 
         self.S = {}
         self.tau_exp = {}
+        self.N_sigma = {}
 
         if kwargs.get('fft') is False:
             fft = False
@@ -193,8 +198,8 @@ class Obs:
             if kwarg_name in kwargs:
                 tmp = kwargs.get(kwarg_name)
                 if isinstance(tmp, (int, float)):
-                    if tmp <= 0:
-                        raise Exception(kwarg_name + ' has to be larger than 0.')
+                    if tmp < 0:
+                        raise Exception(kwarg_name + ' has to be larger or equal to 0.')
                     for e, e_name in enumerate(self.e_names):
                         getattr(self, kwarg_name)[e_name] = tmp
                 else:
@@ -208,13 +213,7 @@ class Obs:
 
         _parse_kwarg('S')
         _parse_kwarg('tau_exp')
-
-        if 'N_sigma' in kwargs:
-            self.N_sigma = kwargs.get('N_sigma')
-            if not isinstance(self.N_sigma, (int, float)):
-                raise TypeError('N_sigma is not a number.')
-        else:
-            self.N_sigma = Obs.N_sigma_global
+        _parse_kwarg('N_sigma')
 
         for e, e_name in enumerate(self.e_names):
 
@@ -267,7 +266,7 @@ class Obs:
                 # Critical slowing down analysis
                 for n in range(1, w_max // 2):
                     _compute_drho(n + 1)
-                    if (self.e_rho[e_name][n] - self.N_sigma * self.e_drho[e_name][n]) < 0 or n >= w_max // 2 - 2:
+                    if (self.e_rho[e_name][n] - self.N_sigma[e_name] * self.e_drho[e_name][n]) < 0 or n >= w_max // 2 - 2:
                         # Bias correction hep-lat/0306017 eq. (49) included
                         self.e_tauint[e_name] = self.e_n_tauint[e_name][n] * (1 + (2 * n + 1) / e_N) / (1 + 1 / e_N) + texp * np.abs(self.e_rho[e_name][n + 1])  # The absolute makes sure, that the tail contribution is always positive
                         self.e_dtauint[e_name] = np.sqrt(self.e_n_dtauint[e_name][n] ** 2 + texp ** 2 * self.e_drho[e_name][n + 1] ** 2)
@@ -358,7 +357,7 @@ class Obs:
                 if len(self.e_names) > 1:
                     print('', e_name, '\t %3.8e +/- %3.8e' % (self.e_dvalue[e_name], self.e_ddvalue[e_name]))
                 if self.tau_exp[e_name] > 0:
-                    print(' t_int\t %3.8e +/- %3.8e tau_exp = %3.2f,  N_sigma = %1.0i' % (self.e_tauint[e_name], self.e_dtauint[e_name], self.tau_exp[e_name], self.N_sigma))
+                    print(' t_int\t %3.8e +/- %3.8e tau_exp = %3.2f,  N_sigma = %1.0i' % (self.e_tauint[e_name], self.e_dtauint[e_name], self.tau_exp[e_name], self.N_sigma[e_name]))
                 else:
                     print(' t_int\t %3.8e +/- %3.8e S = %3.2f' % (self.e_tauint[e_name], self.e_dtauint[e_name], self.S[e_name]))
         if self.tag is not None:

From b90caa4cdc0b77fa32e1c146ea4291fed96b0d39 Mon Sep 17 00:00:00 2001
From: Fabian Joswig <fabian.joswig@ed.ac.uk>
Date: Fri, 12 Nov 2021 14:25:51 +0000
Subject: [PATCH 10/27] further docstrings updated

---
 pyerrors/correlators.py | 21 +++++++++++++++------
 pyerrors/obs.py         |  8 ++++----
 2 files changed, 19 insertions(+), 10 deletions(-)

diff --git a/pyerrors/correlators.py b/pyerrors/correlators.py
index c2ae3eda..9c1ef29c 100644
--- a/pyerrors/correlators.py
+++ b/pyerrors/correlators.py
@@ -79,16 +79,16 @@ class Corr:
         else:
             raise Exception("Reweighting status of correlator corrupted.")
 
-    def gamma_method(self):
+    def gamma_method(self, **kwargs):
         """Apply the gamma method to the content of the Corr."""
         for item in self.content:
             if not(item is None):
                 if self.N == 1:
-                    item[0].gamma_method()
+                    item[0].gamma_method(**kwargs)
                 else:
                     for i in range(self.N):
                         for j in range(self.N):
-                            item[i, j].gamma_method()
+                            item[i, j].gamma_method(**kwargs)
 
     # We need to project the Correlator with a Vector to get a single value at each timeslice.
     # The method can use one or two vectors.
@@ -238,7 +238,15 @@ class Corr:
         return Corr(self.content[::-1])
 
     def correlate(self, partner):
-        """Correlate the correlator with another correlator or Obs"""
+        """Correlate the correlator with another correlator or Obs
+
+        Parameters
+        ----------
+        partner : Obs or Corr
+            partner to correlate the correlator with.
+            Can either be an Obs which is correlated with all entries of the
+            correlator or a Corr of same length.
+        """
         new_content = []
         for x0, t_slice in enumerate(self.content):
             if t_slice is None:
@@ -350,10 +358,11 @@ class Corr:
         Parameters
         ----------
         variant : str
-            log: uses the standard effective mass log(C(t) / C(t+1))
-            cosh : Use periodicitiy of the correlator by solving C(t) / C(t+1) = cosh(m * (t - T/2)) / cosh(m * (t + 1 - T/2)) for m.
+            log : uses the standard effective mass log(C(t) / C(t+1))
+            cosh, periodic : Use periodicitiy of the correlator by solving C(t) / C(t+1) = cosh(m * (t - T/2)) / cosh(m * (t + 1 - T/2)) for m.
             sinh : Use anti-periodicitiy of the correlator by solving C(t) / C(t+1) = sinh(m * (t - T/2)) / sinh(m * (t + 1 - T/2)) for m.
             See, e.g., arXiv:1205.5380
+            arccosh : Uses the explicit form of the symmetrized correlator (not recommended)
         guess : float
             guess for the root finder, only relevant for the root variant
         """
diff --git a/pyerrors/obs.py b/pyerrors/obs.py
index 1c98d2a3..2b4398cf 100644
--- a/pyerrors/obs.py
+++ b/pyerrors/obs.py
@@ -333,10 +333,6 @@ class Obs:
 
         return gamma
 
-    def print(self, level=1):
-        warnings.warn("Method 'print' renamed to 'details'", DeprecationWarning)
-        self.details(level > 1)
-
     def details(self, ens_content=True):
         """Output detailed properties of the Obs.
 
@@ -370,6 +366,10 @@ class Obs:
             m = max(map(len, list(self.e_content.keys()))) + 1
             print('\n'.join(['  ' + key.rjust(m) + ': ' + str(value) for key, value in sorted(self.e_content.items())]))
 
+    def print(self, level=1):
+        warnings.warn("Method 'print' renamed to 'details'", DeprecationWarning)
+        self.details(level > 1)
+
     def is_zero_within_error(self, sigma=1):
         """Checks whether the observable is zero within 'sigma' standard errors.
 

From 2ba4e747601037cd8d378a79fc80150e30c6455c Mon Sep 17 00:00:00 2001
From: Fabian Joswig <fabian.joswig@ed.ac.uk>
Date: Fri, 12 Nov 2021 14:46:45 +0000
Subject: [PATCH 11/27] tests for gamma method kwargs added

---
 tests/{pyerrors_test.py => obs_test.py} | 59 +++++++++++++++++++++++++
 1 file changed, 59 insertions(+)
 rename tests/{pyerrors_test.py => obs_test.py} (89%)

diff --git a/tests/pyerrors_test.py b/tests/obs_test.py
similarity index 89%
rename from tests/pyerrors_test.py
rename to tests/obs_test.py
index 8f2ca2bd..f0c63ac7 100644
--- a/tests/pyerrors_test.py
+++ b/tests/obs_test.py
@@ -117,6 +117,65 @@ def test_gamma_method_persistance():
     assert ddvalue == my_obs.ddvalue
 
 
+def test_gamma_method_kwargs():
+
+    my_obs = pe.Obs([np.random.normal(1, 0.8, 5)], ['ens'], idl=[[1, 2, 3, 6, 17]])
+
+    pe.Obs.S_dict['ens13.7'] = 3
+
+    my_obs.gamma_method()
+    assert my_obs.S['ens'] == pe.Obs.S_global
+    assert my_obs.tau_exp['ens'] == pe.Obs.tau_exp_global
+    assert my_obs.N_sigma['ens'] == pe.Obs.N_sigma_global
+
+    my_obs.gamma_method(S=3.71)
+    assert my_obs.S['ens'] == 3.71
+    assert my_obs.tau_exp['ens'] == pe.Obs.tau_exp_global
+    assert my_obs.N_sigma['ens'] == pe.Obs.N_sigma_global
+
+    my_obs.gamma_method(tau_exp=17)
+    assert my_obs.S['ens'] == pe.Obs.S_global
+    assert my_obs.tau_exp['ens'] == 17
+    assert my_obs.N_sigma['ens'] == pe.Obs.N_sigma_global
+
+    my_obs.gamma_method(tau_exp=1.7, N_sigma=2.123)
+    assert my_obs.S['ens'] == pe.Obs.S_global
+    assert my_obs.tau_exp['ens'] == 1.7
+    assert my_obs.N_sigma['ens'] == 2.123
+
+    pe.Obs.S_dict['ens'] = 3
+    pe.Obs.S_dict['ens|23'] = 7
+
+    my_obs.gamma_method()
+    assert my_obs.S['ens'] == pe.Obs.S_dict['ens'] == 3
+    assert my_obs.tau_exp['ens'] == pe.Obs.tau_exp_global
+    assert my_obs.N_sigma['ens'] == pe.Obs.N_sigma_global
+
+    pe.Obs.tau_exp_dict['ens'] = 4
+    pe.Obs.N_sigma_dict['ens'] = 4
+
+    my_obs.gamma_method()
+    assert my_obs.S['ens'] == pe.Obs.S_dict['ens'] == 3
+    assert my_obs.tau_exp['ens'] == pe.Obs.tau_exp_dict['ens'] == 4
+    assert my_obs.N_sigma['ens'] == pe.Obs.N_sigma_dict['ens'] == 4
+
+    my_obs.gamma_method(S=1.1, tau_exp=1.2, N_sigma=1.3)
+    assert my_obs.S['ens'] == 1.1
+    assert my_obs.tau_exp['ens'] == 1.2
+    assert my_obs.N_sigma['ens'] == 1.3
+
+    pe.Obs.S_dict = {}
+    pe.Obs.tau_exp_dict = {}
+    pe.Obs.N_sigma_dict = {}
+
+    my_obs = pe.Obs([np.random.normal(1, 0.8, 5)], ['ens'])
+
+    my_obs.gamma_method()
+    assert my_obs.S['ens'] == pe.Obs.S_global
+    assert my_obs.tau_exp['ens'] == pe.Obs.tau_exp_global
+    assert my_obs.N_sigma['ens'] == pe.Obs.N_sigma_global
+
+
 def test_covariance_is_variance():
     value = np.random.normal(5, 10)
     dvalue = np.abs(np.random.normal(0, 1))

From db94a49d76539033e23f147d13c73153a5eb71ad Mon Sep 17 00:00:00 2001
From: Fabian Joswig <fabian.joswig@ed.ac.uk>
Date: Fri, 12 Nov 2021 15:13:50 +0000
Subject: [PATCH 12/27] added special case for tau_exp error analysis and less
 than 8 samples

---
 pyerrors/obs.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/pyerrors/obs.py b/pyerrors/obs.py
index 2b4398cf..ed3c23a2 100644
--- a/pyerrors/obs.py
+++ b/pyerrors/obs.py
@@ -264,6 +264,8 @@ class Obs:
                 # if type(self.idl[e_name]) is range: # scale tau_exp according to step size
                 #    texp /= self.idl[e_name].step
                 # Critical slowing down analysis
+                if w_max // 2 <= 1:
+                    raise Exception("Need at least 8 samples for tau_exp error analysis")
                 for n in range(1, w_max // 2):
                     _compute_drho(n + 1)
                     if (self.e_rho[e_name][n] - self.N_sigma[e_name] * self.e_drho[e_name][n]) < 0 or n >= w_max // 2 - 2:

From 1d345476fa46c606f29dde6da72eb62b9f81c7f6 Mon Sep 17 00:00:00 2001
From: Fabian Joswig <fabian.joswig@ed.ac.uk>
Date: Fri, 12 Nov 2021 15:26:53 +0000
Subject: [PATCH 13/27] When initializing Obs with less than 5 samples the
 correct exception message is not given, additional check on length of idl
 added.

---
 pyerrors/obs.py | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/pyerrors/obs.py b/pyerrors/obs.py
index ed3c23a2..87cbf092 100644
--- a/pyerrors/obs.py
+++ b/pyerrors/obs.py
@@ -70,12 +70,15 @@ class Obs:
         if means is None:
             if len(samples) != len(names):
                 raise Exception('Length of samples and names incompatible.')
+            if idl is not None:
+                if len(idl) != len(names):
+                    raise Exception('Length of idl incompatible with samples and names.')
             if len(names) != len(set(names)):
-                raise Exception('Names are not unique.')
+                raise Exception('names are not unique.')
             if not all(isinstance(x, str) for x in names):
                 raise TypeError('All names have to be strings.')
             if min(len(x) for x in samples) <= 4:
-                raise Exception('Samples have to have at least 4 entries.')
+                raise Exception('Samples have to have at least 5 entries.')
 
         self.names = sorted(names)
         self.shape = {}

From 3a265c03ba56cc9bcb13412f5640a58b97bd2c8b Mon Sep 17 00:00:00 2001
From: Fabian Joswig <fabian.joswig@ed.ac.uk>
Date: Fri, 12 Nov 2021 16:21:17 +0000
Subject: [PATCH 14/27] Output of Obs.details now informs the user about
 irregular monte carlo histories

---
 pyerrors/obs.py | 25 +++++++++++++++++++++++--
 1 file changed, 23 insertions(+), 2 deletions(-)

diff --git a/pyerrors/obs.py b/pyerrors/obs.py
index 87cbf092..f98d2a6e 100644
--- a/pyerrors/obs.py
+++ b/pyerrors/obs.py
@@ -368,8 +368,29 @@ class Obs:
                 print(self.N, 'samples in', len(self.e_names), 'ensemble:')
             else:
                 print(self.N, 'samples in', len(self.e_names), 'ensembles:')
-            m = max(map(len, list(self.e_content.keys()))) + 1
-            print('\n'.join(['  ' + key.rjust(m) + ': ' + str(value) for key, value in sorted(self.e_content.items())]))
+            my_string_list = []
+            for key, value in sorted(self.e_content.items()):
+                my_string = '  ' + "\u00B7 Ensemble '" + key + "' "
+                if len(value) == 1:
+                    my_string += f': {self.shape[value[0]]} configurations'
+                    if isinstance(self.idl[value[0]], range):
+                        my_string += f' (from {self.idl[value[0]].start} to {self.idl[value[0]][-1]}' + int(self.idl[value[0]].step != 1) * f' in steps of {self.idl[value[0]].step}' + ')'
+                    else:
+                        my_string += ' (irregular range)'
+                else:
+                    sublist = []
+                    for v in value:
+                        my_substring = '    ' + "\u00B7 Replicum '" + v[len(key) + 1:] + "' "
+                        my_substring += f': {self.shape[v]} configurations'
+                        if isinstance(self.idl[v], range):
+                            my_substring += f' (from {self.idl[v].start} to {self.idl[v][-1]}' + int(self.idl[v].step != 1) * f' in steps of {self.idl[v].step}' + ')'
+                        else:
+                            my_substring += ' (irregular range)'
+                        sublist.append(my_substring)
+
+                    my_string += '\n' + '\n'.join(sublist)
+                my_string_list.append(my_string)
+            print('\n'.join(my_string_list))
 
     def print(self, level=1):
         warnings.warn("Method 'print' renamed to 'details'", DeprecationWarning)

From ee9d43ad26b53d349d0874fe3598fc48e203bbad Mon Sep 17 00:00:00 2001
From: Fabian Joswig <fabian.joswig@ed.ac.uk>
Date: Fri, 12 Nov 2021 16:24:13 +0000
Subject: [PATCH 15/27] description moved to top in Obs.details

---
 pyerrors/obs.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/pyerrors/obs.py b/pyerrors/obs.py
index f98d2a6e..8ef64353 100644
--- a/pyerrors/obs.py
+++ b/pyerrors/obs.py
@@ -346,6 +346,8 @@ class Obs:
         ens_content : bool
             print details about the ensembles and replica if true.
         """
+        if self.tag is not None:
+            print("Description:", self.tag)
         if self.value == 0.0:
             percentage = np.nan
         else:
@@ -361,8 +363,6 @@ class Obs:
                     print(' t_int\t %3.8e +/- %3.8e tau_exp = %3.2f,  N_sigma = %1.0i' % (self.e_tauint[e_name], self.e_dtauint[e_name], self.tau_exp[e_name], self.N_sigma[e_name]))
                 else:
                     print(' t_int\t %3.8e +/- %3.8e S = %3.2f' % (self.e_tauint[e_name], self.e_dtauint[e_name], self.S[e_name]))
-        if self.tag is not None:
-            print("Description:", self.tag)
         if ens_content is True:
             if len(self.e_names) == 1:
                 print(self.N, 'samples in', len(self.e_names), 'ensemble:')

From 1ca3055460cd680839c516b608edd4fa7eb86d06 Mon Sep 17 00:00:00 2001
From: Fabian Joswig <fabian.joswig@ed.ac.uk>
Date: Mon, 15 Nov 2021 09:54:42 +0000
Subject: [PATCH 16/27] docs: output of Obs.details updated

---
 pyerrors/__init__.py | 17 +++++++++--------
 1 file changed, 9 insertions(+), 8 deletions(-)

diff --git a/pyerrors/__init__.py b/pyerrors/__init__.py
index c639c887..74f87b14 100644
--- a/pyerrors/__init__.py
+++ b/pyerrors/__init__.py
@@ -76,11 +76,10 @@ obs2 = pe.Obs([samples2], ['ensemble2'])
 
 my_sum = obs1 + obs2
 my_sum.details()
-> Result	 2.00596631e+00 +/- 0.00000000e+00 +/- 0.00000000e+00 (0.000%)
+> Result   2.00697958e+00 +/- 0.00000000e+00 +/- 0.00000000e+00 (0.000%)
 > 1500 samples in 2 ensembles:
->    ensemble1: ['ensemble1']
->    ensemble2: ['ensemble2']
-
+>   · Ensemble 'ensemble1' : 1000 configurations (from 1 to 1000)
+>   · Ensemble 'ensemble2' : 500 configurations (from 1 to 500)
 ```
 
 `pyerrors` identifies multiple replica (independent Markov chains with identical simulation parameters) by the vertical bar `|` in the name of the dataset.
@@ -90,11 +89,13 @@ Example:
 obs1 = pe.Obs([samples1], ['ensemble1|r01'])
 obs2 = pe.Obs([samples2], ['ensemble1|r02'])
 
-my_sum = obs1 + obs2
-my_sum.details()
-> Result	 2.00596631e+00 +/- 0.00000000e+00 +/- 0.00000000e+00 (0.000%)
+> my_sum = obs1 + obs2
+> my_sum.details()
+> Result   2.00697958e+00 +/- 0.00000000e+00 +/- 0.00000000e+00 (0.000%)
 > 1500 samples in 1 ensemble:
->    ensemble1: ['ensemble1|r01', 'ensemble1|r02']
+>   · Ensemble 'ensemble1'
+>     · Replicum 'r01' : 1000 configurations (from 1 to 1000)
+>     · Replicum 'r02' : 500 configurations (from 1 to 500)
 ```
 ## Irregular Monte Carlo chains
 

From dfd5eafe125ffc9c7c1c95b29a09fd03e75896c4 Mon Sep 17 00:00:00 2001
From: Fabian Joswig <fabian.joswig@ed.ac.uk>
Date: Mon, 15 Nov 2021 10:11:15 +0000
Subject: [PATCH 17/27] docs: gamma_method doc extended

---
 pyerrors/__init__.py | 46 +++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 45 insertions(+), 1 deletion(-)

diff --git a/pyerrors/__init__.py b/pyerrors/__init__.py
index 74f87b14..d121f1f5 100644
--- a/pyerrors/__init__.py
+++ b/pyerrors/__init__.py
@@ -61,9 +61,41 @@ my_m_eff = np.log(my_obs1 / my_obs2)
 
 The error propagation is based on the gamma method introduced in [arXiv:hep-lat/0306017](https://arxiv.org/abs/hep-lat/0306017).
 
+After having arrived at
+
+Example:
+```python
+my_sum.gamma_method()
+my_sum.details()
+```
+
+The standard value for the automatic windowing procedure is $S=2$. Other values for $S$ can be passed to the `gamma_method` as parameter.
+
+Example:
+```python
+my_sum.gamma_method(S=3.0)
+my_sum.details()
+```
+
+The integrated autocorrelation time $\tau_\mathrm{int}$ and the autocorrelation function $\rho(W)$ can be monitored via the methods ´pyerrors.obs.Obs.plot_tauint` and ´pyerrors.obs.Obs.plot_tauint`.
+
+Example:
+```python
+my_sum.plot_tauint()
+my_sum.plot_rho()
+```
+
+### Exponential tails
+
+Slow modes in the Monte Carlo history can be accounted for by attaching and exponntial tail to the autocorrelation function $\rho$ as suggested in [arXiv:1009.5228](https://arxiv.org/abs/1009.5228). The longest autocorrelation time in the history, $\tau_\mathrm{exp}$, can be passed to the `gamma_method` as parameter. In this case the automatic windowing procedure is vacated and the parameter $S$ does not affect the error estimate.
+
+Example:
+```python
+my_sum.gamma_method(tau_exp=4.2)
+my_sum.details()
+```
 
 For the full API see `pyerrors.obs.Obs.gamma_method`
-### Exponential tails
 
 ## Multiple ensembles/replica
 
@@ -97,6 +129,18 @@ obs2 = pe.Obs([samples2], ['ensemble1|r02'])
 >     · Replicum 'r01' : 1000 configurations (from 1 to 1000)
 >     · Replicum 'r02' : 500 configurations (from 1 to 500)
 ```
+
+### Error estimation for multiple ensembles
+
+In order to keep track of different error analyis parameters for different ensembles one can make use of global dictionaries as detailed in the following example.
+
+Example:
+```python
+pe.Obs.S_dict['ensemble1'] = 2.5
+pe.Obs.tau_exp_dict['ensemble2'] = 8.0
+pe.Obs.tau_exp_dict['ensemble3'] = 2.0
+```
+
 ## Irregular Monte Carlo chains
 
 Irregular Monte Carlo chains can be initilized with the parameter `idl`.

From ab19abf84fae35cd99a6f6961883fbad357d9f67 Mon Sep 17 00:00:00 2001
From: Fabian Joswig <fabian.joswig@ed.ac.uk>
Date: Mon, 15 Nov 2021 10:42:34 +0000
Subject: [PATCH 18/27] docs: gamma_method explanation extended

---
 pyerrors/__init__.py | 19 ++++++++++++++++++-
 1 file changed, 18 insertions(+), 1 deletion(-)

diff --git a/pyerrors/__init__.py b/pyerrors/__init__.py
index d121f1f5..51d778c4 100644
--- a/pyerrors/__init__.py
+++ b/pyerrors/__init__.py
@@ -67,6 +67,10 @@ Example:
 ```python
 my_sum.gamma_method()
 my_sum.details()
+> Result	 1.70000000e+00 +/- 3.89934513e+00 +/- 5.84901770e-01 (229.373%)
+>  t_int	 3.72133617e+00 +/- 9.81032454e-01 S = 2.00
+> 1000 samples in 1 ensemble:
+>   · Ensemble 'ensemble_name' : 1000 configurations (from 1 to 1000)
 ```
 
 The standard value for the automatic windowing procedure is $S=2$. Other values for $S$ can be passed to the `gamma_method` as parameter.
@@ -75,6 +79,11 @@ Example:
 ```python
 my_sum.gamma_method(S=3.0)
 my_sum.details()
+> Result	 1.70000000e+00 +/- 3.77151850e+00 +/- 6.47779576e-01 (221.854%)
+>  t_int	 3.48135280e+00 +/- 1.06547679e+00 S = 3.00
+> 1000 samples in 1 ensemble:
+>   · Ensemble 'ensemble_name' : 1000 configurations (from 1 to 1000)
+
 ```
 
 The integrated autocorrelation time $\tau_\mathrm{int}$ and the autocorrelation function $\rho(W)$ can be monitored via the methods ´pyerrors.obs.Obs.plot_tauint` and ´pyerrors.obs.Obs.plot_tauint`.
@@ -91,8 +100,12 @@ Slow modes in the Monte Carlo history can be accounted for by attaching and expo
 
 Example:
 ```python
-my_sum.gamma_method(tau_exp=4.2)
+my_sum.gamma_method(tau_exp=7.2)
 my_sum.details()
+> Result	 1.70000000e+00 +/- 3.77806247e+00 +/- 3.48320149e-01 (222.239%)
+>  t_int	 3.49344429e+00 +/- 7.62747210e-01 tau_exp = 7.20,  N_sigma = 1
+> 1000 samples in 1 ensemble:
+>   · Ensemble 'ensemble_name' : 1000 configurations (from 1 to 1000)
 ```
 
 For the full API see `pyerrors.obs.Obs.gamma_method`
@@ -141,6 +154,10 @@ pe.Obs.tau_exp_dict['ensemble2'] = 8.0
 pe.Obs.tau_exp_dict['ensemble3'] = 2.0
 ```
 
+In case the `gamma_method` is called without any parameters it will use the values specified in the dictionaries for the respective ensembles.
+Passing arguments to the `gamma_method` still dominates over the dictionaries.
+
+
 ## Irregular Monte Carlo chains
 
 Irregular Monte Carlo chains can be initilized with the parameter `idl`.

From 24ef100197eb6f542e4087c79960496f2e6f5bac Mon Sep 17 00:00:00 2001
From: Fabian Joswig <fabian.joswig@ed.ac.uk>
Date: Mon, 15 Nov 2021 10:56:25 +0000
Subject: [PATCH 19/27] docs: gamma_method explanation extended, typos fixed

---
 pyerrors/__init__.py | 18 +++++++++---------
 1 file changed, 9 insertions(+), 9 deletions(-)

diff --git a/pyerrors/__init__.py b/pyerrors/__init__.py
index 51d778c4..be67e231 100644
--- a/pyerrors/__init__.py
+++ b/pyerrors/__init__.py
@@ -5,7 +5,7 @@ It is based on the **gamma method** [arXiv:hep-lat/0306017](https://arxiv.org/ab
 - **automatic differentiation** as suggested in [arXiv:1809.01289](https://arxiv.org/abs/1809.01289) (partly based on the [autograd](https://github.com/HIPS/autograd) package)
 - **treatment of slow modes** in the simulation as suggested in [arXiv:1009.5228](https://arxiv.org/abs/1009.5228)
 - coherent **error propagation** for data from **different Markov chains**
-- **non-linear fits with x- and y-errors** and exact linear error propagation based on automatic differentiation as introduced in [arXiv:1809.01289]
+- **non-linear fits with x- and y-errors** and exact linear error propagation based on automatic differentiation as introduced in [arXiv:1809.01289](https://arxiv.org/abs/1809.01289)
 - **real and complex matrix operations** and their error propagation based on automatic differentiation (cholesky decomposition, calculation of eigenvalues and eigenvectors, singular value decomposition...)
 
 ## Getting started
@@ -23,7 +23,7 @@ print(my_new_obs)
 # The `Obs` class
 
 `pyerrors` introduces a new datatype, `Obs`, which simplifies error propagation and estimation for auto- and cross-correlated data.
-An `Obs` object can be initialized with two arguments, the first is a list containining the samples for an Observable from a Monte Carlo chain.
+An `Obs` object can be initialized with two arguments, the first is a list containing the samples for an Observable from a Monte Carlo chain.
 The samples can either be provided as python list or as numpy array.
 The second argument is a list containing the names of the respective Monte Carlo chains as strings. These strings uniquely identify a Monte Carlo chain/ensemble.
 
@@ -59,9 +59,9 @@ my_m_eff = np.log(my_obs1 / my_obs2)
 
 ## Error estimation
 
-The error propagation is based on the gamma method introduced in [arXiv:hep-lat/0306017](https://arxiv.org/abs/hep-lat/0306017).
+The error estimation within `pyerrors` is based on the gamma method introduced in [arXiv:hep-lat/0306017](https://arxiv.org/abs/hep-lat/0306017).
 
-After having arrived at
+After having arrived at the derived quantity of interest the `gamma_method` can be called as detailed in the following example.
 
 Example:
 ```python
@@ -96,7 +96,7 @@ my_sum.plot_rho()
 
 ### Exponential tails
 
-Slow modes in the Monte Carlo history can be accounted for by attaching and exponntial tail to the autocorrelation function $\rho$ as suggested in [arXiv:1009.5228](https://arxiv.org/abs/1009.5228). The longest autocorrelation time in the history, $\tau_\mathrm{exp}$, can be passed to the `gamma_method` as parameter. In this case the automatic windowing procedure is vacated and the parameter $S$ does not affect the error estimate.
+Slow modes in the Monte Carlo history can be accounted for by attaching an exponential tail to the autocorrelation function $\rho$ as suggested in [arXiv:1009.5228](https://arxiv.org/abs/1009.5228). The longest autocorrelation time in the history, $\tau_\mathrm{exp}$, can be passed to the `gamma_method` as parameter. In this case the automatic windowing procedure is vacated and the parameter $S$ does not affect the error estimate.
 
 Example:
 ```python
@@ -112,7 +112,7 @@ For the full API see `pyerrors.obs.Obs.gamma_method`
 
 ## Multiple ensembles/replica
 
-Error propagation for multiple ensembles (Markov chains with different simulation parameters) is handeled automatically. Ensembles are uniquely identified by their `name`.
+Error propagation for multiple ensembles (Markov chains with different simulation parameters) is handled automatically. Ensembles are uniquely identified by their `name`.
 
 Example:
 ```python
@@ -145,7 +145,7 @@ obs2 = pe.Obs([samples2], ['ensemble1|r02'])
 
 ### Error estimation for multiple ensembles
 
-In order to keep track of different error analyis parameters for different ensembles one can make use of global dictionaries as detailed in the following example.
+In order to keep track of different error analysis parameters for different ensembles one can make use of global dictionaries as detailed in the following example.
 
 Example:
 ```python
@@ -160,7 +160,7 @@ Passing arguments to the `gamma_method` still dominates over the dictionaries.
 
 ## Irregular Monte Carlo chains
 
-Irregular Monte Carlo chains can be initilized with the parameter `idl`.
+Irregular Monte Carlo chains can be initialized with the parameter `idl`.
 
 Example:
 ```python
@@ -175,7 +175,7 @@ obs3 = pe.Obs([samples3], ['ensemble1'], idl=[[2, 9, 28, 29, 501]])
 ```
 
 **Warning:** Irregular Monte Carlo chains can result in odd patterns in the autocorrelation functions.
-Make sure to check the with e.g. `pyerrors.obs.Obs.plot_rho` or `pyerrors.obs.Obs.plot_tauint`.
+Make sure to check the autocorrelation time with e.g. `pyerrors.obs.Obs.plot_rho` or `pyerrors.obs.Obs.plot_tauint`.
 
 For the full API see `pyerrors.obs.Obs`
 

From b937ef04bbc2b32f3c628bf1dfcdcee88bda261e Mon Sep 17 00:00:00 2001
From: Fabian Joswig <fabian.joswig@ed.ac.uk>
Date: Mon, 15 Nov 2021 11:13:12 +0000
Subject: [PATCH 20/27] docs: typos in docstrings corrected

---
 pyerrors/correlators.py | 12 ++++++------
 pyerrors/fits.py        |  4 ++--
 pyerrors/linalg.py      |  2 +-
 pyerrors/obs.py         |  8 ++++----
 4 files changed, 13 insertions(+), 13 deletions(-)

diff --git a/pyerrors/correlators.py b/pyerrors/correlators.py
index 9c1ef29c..7e333cdb 100644
--- a/pyerrors/correlators.py
+++ b/pyerrors/correlators.py
@@ -14,7 +14,7 @@ class Corr:
 
     Everything, this class does, can be achieved using lists or arrays of Obs.
     But it is simply more convenient to have a dedicated object for correlators.
-    One often wants to add or multiply correlators of the same length at every timeslice and it is inconvinient
+    One often wants to add or multiply correlators of the same length at every timeslice and it is inconvenient
     to iterate over all timeslices for every operation. This is especially true, when dealing with smearing matrices.
 
     The correlator can have two types of content: An Obs at every timeslice OR a GEVP
@@ -317,7 +317,7 @@ class Corr:
         Parameters
         ----------
         symmetric : bool
-            decides whether symmertic of simple finite differences are used. Default: True
+            decides whether symmetric of simple finite differences are used. Default: True
         """
         if not symmetric:
             newcontent = []
@@ -415,7 +415,7 @@ class Corr:
             return np.arccosh(Corr(newcontent, padding_back=1, padding_front=1))
 
         else:
-            raise Exception('Unkown variant.')
+            raise Exception('Unknown variant.')
 
     def fit(self, function, fitrange=None, silent=False, **kwargs):
         """Fits function to the data
@@ -433,7 +433,7 @@ class Corr:
         if self.N != 1:
             raise Exception("Correlator must be projected before fitting")
 
-        # The default behaviour is:
+        # The default behavior is:
         # 1 use explicit fitrange
         # if none is provided, use the range of the corr
         # if this is also not set, use the whole length of the corr (This could come with a warning!)
@@ -451,7 +451,7 @@ class Corr:
         return result
 
     def plateau(self, plateau_range=None, method="fit"):
-        """ Extract a plateu value from a Corr object
+        """ Extract a plateau value from a Corr object
 
         Parameters
         ----------
@@ -582,7 +582,7 @@ class Corr:
         return
 
     def dump(self, filename):
-        """Dumps the Corr into a pickel file
+        """Dumps the Corr into a pickle file
 
         Parameters
         ----------
diff --git a/pyerrors/fits.py b/pyerrors/fits.py
index 686ed81f..291a5568 100644
--- a/pyerrors/fits.py
+++ b/pyerrors/fits.py
@@ -21,7 +21,7 @@ class Fit_result(Sequence):
     ----------
     fit_parameters : list
         results for the individual fit parameters,
-        also accesible via indices.
+        also accessible via indices.
     """
 
     def __init__(self):
@@ -427,7 +427,7 @@ def _standard_fit(x, y, func, silent=False, **kwargs):
         raise Exception('x and y input have to have the same length')
 
     if len(x.shape) > 2:
-        raise Exception('Unkown format for x values')
+        raise Exception('Unknown format for x values')
 
     if not callable(func):
         raise TypeError('func has to be a function.')
diff --git a/pyerrors/linalg.py b/pyerrors/linalg.py
index c1a042fb..d993f291 100644
--- a/pyerrors/linalg.py
+++ b/pyerrors/linalg.py
@@ -175,7 +175,7 @@ def inv(x):
 
 
 def cholesky(x):
-    """Cholesky decompostion of Obs or CObs valued matrices."""
+    """Cholesky decomposition of Obs or CObs valued matrices."""
     return _mat_mat_op(anp.linalg.cholesky, x)
 
 
diff --git a/pyerrors/obs.py b/pyerrors/obs.py
index 8ef64353..6f29d01d 100644
--- a/pyerrors/obs.py
+++ b/pyerrors/obs.py
@@ -59,7 +59,7 @@ class Obs:
         samples : list
             list of numpy arrays containing the Monte Carlo samples
         names : list
-            list of strings labeling the indivdual samples
+            list of strings labeling the individual samples
         idl : list, optional
             list of ranges or lists on which the samples are defined
         means : list, optional
@@ -314,9 +314,9 @@ class Obs:
         deltas : list
             List of fluctuations
         idx : list
-            List or range of configs on which the deltas are defined.
+            List or range of configurations on which the deltas are defined.
         shape : int
-            Number of configs in idx.
+            Number of configurations in idx.
         w_max : int
             Upper bound for the summation window.
         fft : bool
@@ -1109,7 +1109,7 @@ def _reduce_deltas(deltas, idx_old, idx_new):
         Has to be a subset of idx_old.
     """
     if not len(deltas) == len(idx_old):
-        raise Exception('Lenght of deltas and idx_old have to be the same: %d != %d' % (len(deltas), len(idx_old)))
+        raise Exception('Length of deltas and idx_old have to be the same: %d != %d' % (len(deltas), len(idx_old)))
     if type(idx_old) is range and type(idx_new) is range:
         if idx_old == idx_new:
             return deltas

From 02a77f4924d7a8c103e790ad94e6070140f6fe8a Mon Sep 17 00:00:00 2001
From: Fabian Joswig <fabian.joswig@ed.ac.uk>
Date: Mon, 15 Nov 2021 13:12:02 +0000
Subject: [PATCH 21/27] docs: getting started updated

---
 pyerrors/__init__.py | 10 ++++++++--
 1 file changed, 8 insertions(+), 2 deletions(-)

diff --git a/pyerrors/__init__.py b/pyerrors/__init__.py
index be67e231..8eb05e35 100644
--- a/pyerrors/__init__.py
+++ b/pyerrors/__init__.py
@@ -15,11 +15,17 @@ import numpy as np
 import pyerrors as pe
 
 my_obs = pe.Obs([samples], ['ensemble_name'])
-my_new_obs = 2 * np.log(my_obs) / my_obs
+my_new_obs = 2 * np.log(my_obs) / my_obs ** 2
 my_new_obs.gamma_method()
-my_new_obs.details()
 print(my_new_obs)
+> 0.31498(72)
+
+iamzero = my_new_obs - my_new_obs
+iamzero.gamma_method()
+print(iamzero)
+> 0.0
 ```
+
 # The `Obs` class
 
 `pyerrors` introduces a new datatype, `Obs`, which simplifies error propagation and estimation for auto- and cross-correlated data.

From 69d786c1554b66707b8503b06210de5c695cc34a Mon Sep 17 00:00:00 2001
From: Fabian Joswig <fabian.joswig@ed.ac.uk>
Date: Mon, 15 Nov 2021 14:11:45 +0000
Subject: [PATCH 22/27] feat: export_jackknife method added to Obs class

---
 pyerrors/obs.py | 27 +++++++++++++++++++++++++++
 1 file changed, 27 insertions(+)

diff --git a/pyerrors/obs.py b/pyerrors/obs.py
index 6f29d01d..d73f92e3 100644
--- a/pyerrors/obs.py
+++ b/pyerrors/obs.py
@@ -557,6 +557,33 @@ class Obs:
         with open(file_name, 'wb') as fb:
             pickle.dump(self, fb)
 
+    def export_jackknife(self):
+        """Export jackknife samples from the Obs
+
+        Returns
+        -------
+        np.ndarray
+            Returns a numpy array of length N + 1 where N is the number of samples
+            for the given ensemble and replicum. The zeroth entry of the array contains
+            the mean value of the Obs, entries 1 to N contain the N jackknife samples
+            derived from the Obs. The current implementation only works for observables
+            defined on exactly one ensemble and replicum.
+        """
+
+        if len(self.names) != 1:
+            raise Exception("'export_jackknife' is only implemented for Obs defined on one ensemble and replicum.")
+
+        name = self.names[0]
+        full_data = self.deltas[name] + self.r_values[name]
+        n = full_data.size
+        mean = np.mean(full_data)
+        tmp_jacks = np.zeros(n + 1)
+        tmp_jacks[0] = self.value
+        for i in range(n):
+            tmp_jacks[i + 1] = (n * mean - full_data[i]) / (n - 1)
+
+        return tmp_jacks
+
     def __float__(self):
         return float(self.value)
 

From 331d7cb4716ec2dc4601b50d4aaf24634831d6dc Mon Sep 17 00:00:00 2001
From: Fabian Joswig <fabian.joswig@ed.ac.uk>
Date: Mon, 15 Nov 2021 14:18:26 +0000
Subject: [PATCH 23/27] test_jackknife added

---
 pyerrors/obs.py   |  2 +-
 tests/obs_test.py | 15 +++++++++++++++
 2 files changed, 16 insertions(+), 1 deletion(-)

diff --git a/pyerrors/obs.py b/pyerrors/obs.py
index d73f92e3..2204a429 100644
--- a/pyerrors/obs.py
+++ b/pyerrors/obs.py
@@ -562,7 +562,7 @@ class Obs:
 
         Returns
         -------
-        np.ndarray
+        numpy.ndarray
             Returns a numpy array of length N + 1 where N is the number of samples
             for the given ensemble and replicum. The zeroth entry of the array contains
             the mean value of the Obs, entries 1 to N contain the N jackknife samples
diff --git a/tests/obs_test.py b/tests/obs_test.py
index f0c63ac7..7ce035ab 100644
--- a/tests/obs_test.py
+++ b/tests/obs_test.py
@@ -504,3 +504,18 @@ def test_covariance2_symmetry():
     cov_ba = pe.covariance2(a, test_obs1)
     assert np.abs(cov_ab - cov_ba) <= 10 * np.finfo(np.float64).eps
     assert np.abs(cov_ab) < test_obs1.dvalue * test_obs2.dvalue * (1 + 10 * np.finfo(np.float64).eps)
+
+
+def test_jackknife():
+    full_data = np.random.normal(1.1, 0.87, 5487)
+
+    my_obs = pe.Obs([full_data], ['test'])
+
+    n = full_data.size
+    mean = np.mean(full_data)
+    tmp_jacks = np.zeros(n + 1)
+    tmp_jacks[0] = mean
+    for i in range(n):
+        tmp_jacks[i + 1] = (n * mean - full_data[i]) / (n - 1)
+
+    assert np.allclose(tmp_jacks, my_obs.export_jackknife())

From 1fabdb07519f0c44e11c35ae4874583d87789dad Mon Sep 17 00:00:00 2001
From: Fabian Joswig <fabian.joswig@ed.ac.uk>
Date: Mon, 15 Nov 2021 14:29:39 +0000
Subject: [PATCH 24/27] docs: export_jackknife added to docs mainpage

---
 pyerrors/__init__.py | 9 ++++++++-
 1 file changed, 8 insertions(+), 1 deletion(-)

diff --git a/pyerrors/__init__.py b/pyerrors/__init__.py
index 8eb05e35..a298f982 100644
--- a/pyerrors/__init__.py
+++ b/pyerrors/__init__.py
@@ -133,7 +133,7 @@ my_sum.details()
 >   · Ensemble 'ensemble2' : 500 configurations (from 1 to 500)
 ```
 
-`pyerrors` identifies multiple replica (independent Markov chains with identical simulation parameters) by the vertical bar `|` in the name of the dataset.
+`pyerrors` identifies multiple replica (independent Markov chains with identical simulation parameters) by the vertical bar `|` in the name of the data set.
 
 Example:
 ```python
@@ -198,6 +198,13 @@ For the full API see `pyerrors.correlators.Corr`
 # Matrix operations
 `pyerrors.linalg`
 
+# Export data
+The preferred exported file format within `pyerrors` is
+
+## Jackknife samples
+For comparison with other analysis workflows `pyerrors` can generate jackknife samples from an `Obs` object.
+See `pyerrors.obs.Obs.export_jackknife` for details.
+
 # Input
 `pyerrors.input`
 '''

From 85f37a7abfbc3bf304ca4ef3593dbaeda1e9bd11 Mon Sep 17 00:00:00 2001
From: Fabian Joswig <fabian.joswig@ed.ac.uk>
Date: Mon, 15 Nov 2021 14:38:21 +0000
Subject: [PATCH 25/27] docs: details about jackknife conversion added

---
 pyerrors/obs.py   | 3 ++-
 tests/obs_test.py | 3 +++
 2 files changed, 5 insertions(+), 1 deletion(-)

diff --git a/pyerrors/obs.py b/pyerrors/obs.py
index 2204a429..4fbda90f 100644
--- a/pyerrors/obs.py
+++ b/pyerrors/obs.py
@@ -567,7 +567,8 @@ class Obs:
             for the given ensemble and replicum. The zeroth entry of the array contains
             the mean value of the Obs, entries 1 to N contain the N jackknife samples
             derived from the Obs. The current implementation only works for observables
-            defined on exactly one ensemble and replicum.
+            defined on exactly one ensemble and replicum. The derived jackknife samples
+            should agree with samples from a full jackknife analysis up to O(1/N).
         """
 
         if len(self.names) != 1:
diff --git a/tests/obs_test.py b/tests/obs_test.py
index 7ce035ab..f0adcc49 100644
--- a/tests/obs_test.py
+++ b/tests/obs_test.py
@@ -519,3 +519,6 @@ def test_jackknife():
         tmp_jacks[i + 1] = (n * mean - full_data[i]) / (n - 1)
 
     assert np.allclose(tmp_jacks, my_obs.export_jackknife())
+    my_new_obs = my_obs + pe.Obs([full_data], ['test2'])
+    with pytest.raises(Exception):
+        my_new_obs.export_jackknife()

From 8bb1a910c613329c9015cf7fbddb913909497d05 Mon Sep 17 00:00:00 2001
From: Fabian Joswig <fabian.joswig@ed.ac.uk>
Date: Mon, 15 Nov 2021 14:57:58 +0000
Subject: [PATCH 26/27] docs: documentation of CObs extended

---
 pyerrors/__init__.py | 27 ++++++++++++++++++++++++---
 1 file changed, 24 insertions(+), 3 deletions(-)

diff --git a/pyerrors/__init__.py b/pyerrors/__init__.py
index a298f982..b7f8f1ed 100644
--- a/pyerrors/__init__.py
+++ b/pyerrors/__init__.py
@@ -22,8 +22,8 @@ print(my_new_obs)
 
 iamzero = my_new_obs - my_new_obs
 iamzero.gamma_method()
-print(iamzero)
-> 0.0
+print(iamzero == 0.0)
+> True
 ```
 
 # The `Obs` class
@@ -189,7 +189,28 @@ For the full API see `pyerrors.obs.Obs`
 For the full API see `pyerrors.correlators.Corr`
 
 # Complex observables
-`pyerrors.obs.CObs`
+
+`pyerrors` can handle complex valued observables via the class `pyerrors.obs.CObs`.
+`CObs` are initialized with a real and an imaginary part which both can be `Obs` valued.
+
+Example:
+```python
+my_real_part = pe.Obs([samples1], ['ensemble1'])
+my_imag_part = pe.Obs([samples2], ['ensemble1'])
+
+my_cobs = pe.CObs(my_real_part, my_imag_part)
+my_cobs.gamma_method()
+print(my_cobs)
+> (0.9959(91)+0.659(28)j)
+```
+
+Elementary mathematical operations are overloaded and samples are properly propagated as for the `Obs` class.
+```python
+my_derived_cobs = (my_cobs + my_cobs.conjugate()) / np.abs(my_cobs)
+my_derived_cobs.gamma_method()
+print(my_derived_cobs)
+> (1.668(23)+0.0j)
+```
 
 # Optimization / fits / roots
 `pyerrors.fits`

From 78c3405fb76a9b9713b2470252a05fa6a8f0c751 Mon Sep 17 00:00:00 2001
From: Fabian Joswig <fabian.joswig@ed.ac.uk>
Date: Mon, 15 Nov 2021 15:16:26 +0000
Subject: [PATCH 27/27] docs: typos fixed

---
 pyerrors/__init__.py | 18 ++++++++++--------
 1 file changed, 10 insertions(+), 8 deletions(-)

diff --git a/pyerrors/__init__.py b/pyerrors/__init__.py
index b7f8f1ed..2f194514 100644
--- a/pyerrors/__init__.py
+++ b/pyerrors/__init__.py
@@ -6,7 +6,7 @@ It is based on the **gamma method** [arXiv:hep-lat/0306017](https://arxiv.org/ab
 - **treatment of slow modes** in the simulation as suggested in [arXiv:1009.5228](https://arxiv.org/abs/1009.5228)
 - coherent **error propagation** for data from **different Markov chains**
 - **non-linear fits with x- and y-errors** and exact linear error propagation based on automatic differentiation as introduced in [arXiv:1809.01289](https://arxiv.org/abs/1809.01289)
-- **real and complex matrix operations** and their error propagation based on automatic differentiation (cholesky decomposition, calculation of eigenvalues and eigenvectors, singular value decomposition...)
+- **real and complex matrix operations** and their error propagation based on automatic differentiation (Cholesky decomposition, calculation of eigenvalues and eigenvectors, singular value decomposition...)
 
 ## Getting started
 
@@ -66,17 +66,19 @@ my_m_eff = np.log(my_obs1 / my_obs2)
 ## Error estimation
 
 The error estimation within `pyerrors` is based on the gamma method introduced in [arXiv:hep-lat/0306017](https://arxiv.org/abs/hep-lat/0306017).
-
 After having arrived at the derived quantity of interest the `gamma_method` can be called as detailed in the following example.
 
 Example:
 ```python
 my_sum.gamma_method()
+print(my_sum)
+> 1.70(57)
 my_sum.details()
-> Result	 1.70000000e+00 +/- 3.89934513e+00 +/- 5.84901770e-01 (229.373%)
->  t_int	 3.72133617e+00 +/- 9.81032454e-01 S = 2.00
+> Result	 1.70000000e+00 +/- 5.72046658e-01 +/- 7.56746598e-02 (33.650%)
+>  t_int	 2.71422900e+00 +/- 6.40320983e-01 S = 2.00
 > 1000 samples in 1 ensemble:
 >   · Ensemble 'ensemble_name' : 1000 configurations (from 1 to 1000)
+
 ```
 
 The standard value for the automatic windowing procedure is $S=2$. Other values for $S$ can be passed to the `gamma_method` as parameter.
@@ -85,8 +87,8 @@ Example:
 ```python
 my_sum.gamma_method(S=3.0)
 my_sum.details()
-> Result	 1.70000000e+00 +/- 3.77151850e+00 +/- 6.47779576e-01 (221.854%)
->  t_int	 3.48135280e+00 +/- 1.06547679e+00 S = 3.00
+> Result	 1.70000000e+00 +/- 6.30675201e-01 +/- 1.04585650e-01 (37.099%)
+>  t_int	 3.29909703e+00 +/- 9.77310102e-01 S = 3.00
 > 1000 samples in 1 ensemble:
 >   · Ensemble 'ensemble_name' : 1000 configurations (from 1 to 1000)
 
@@ -108,8 +110,8 @@ Example:
 ```python
 my_sum.gamma_method(tau_exp=7.2)
 my_sum.details()
-> Result	 1.70000000e+00 +/- 3.77806247e+00 +/- 3.48320149e-01 (222.239%)
->  t_int	 3.49344429e+00 +/- 7.62747210e-01 tau_exp = 7.20,  N_sigma = 1
+> Result	 1.70000000e+00 +/- 6.28097762e-01 +/- 5.79077524e-02 (36.947%)
+>  t_int	 3.27218667e+00 +/- 7.99583654e-01 tau_exp = 7.20,  N_sigma = 1
 > 1000 samples in 1 ensemble:
 >   · Ensemble 'ensemble_name' : 1000 configurations (from 1 to 1000)
 ```