Merge branch 'develop' into feature/irregularMC

pyerrors/correlators.py

@@ -5,7 +5,7 @@ import matplotlib.pyplot as plt
 import scipy.linalg
 from .pyerrors import Obs, dump_object
 from .fits import standard_fit
-from .linalg import eigh, mat_mat_op
+from .linalg import eigh, inv, cholesky
 from .roots import find_root
 
 
@@ -187,10 +187,10 @@ class Corr:
     def Eigenvalue(self, t0, state=1):
         G = self.smearing_symmetric()
         G0 = G.content[t0]
-        L = mat_mat_op(anp.linalg.cholesky, G0)
-        Li = mat_mat_op(anp.linalg.inv, L)
+        L = cholesky(G0)
+        Li = inv(L)
         LT = L.T
-        LTi = mat_mat_op(anp.linalg.inv, LT)
+        LTi = inv(LT)
         newcontent = []
         for t in range(self.T):
             Gt = G.content[t]

pyerrors/input/bdio.py

@@ -175,7 +175,7 @@ def write_ADerrors(obs_list, file_path, bdio_path='./libbdio.so', **kwargs):
     """
 
     for obs in obs_list:
-        if not obs.e_names:
+        if not hasattr(obs, 'e_names'):
             raise Exception('Run the gamma method first for all obs.')
 
     bdio = ctypes.cdll.LoadLibrary(bdio_path)

pyerrors/input/hadrons.py

@@ -109,3 +109,60 @@ def read_ExternalLeg_hd5(path, filestem, ens_id, order='F'):
         matrix[si, sj, ci, cj].gamma_method()
 
     return Npr_matrix(matrix.swapaxes(1, 2).reshape((12, 12), order=order), mom_in=mom)
+
+
+def read_Bilinear_hd5(path, filestem, ens_id, order='F'):
+    """Read hadrons Bilinear hdf5 file and output an array of CObs
+
+    Parameters
+    -----------------
+    path -- path to the files to read
+    filestem -- namestem of the files to read
+    ens_id -- name of the ensemble, required for internal bookkeeping
+    order -- order in which the array is to be reshaped,
+             'F' for the first index changing fastest (9 4x4 matrices) default.
+             'C' for the last index changing fastest (16 3x3 matrices),
+    """
+
+    files = _get_files(path, filestem)
+
+    mom_in = None
+    mom_out = None
+
+    corr_data = {}
+    for hd5_file in files:
+        file = h5py.File(path + '/' + hd5_file, "r")
+        for i in range(16):
+            name = file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['gamma'][0].decode('UTF-8')
+            if name not in corr_data:
+                corr_data[name] = []
+            raw_data = file['Bilinear/Bilinear_' + str(i) + '/corr'][0][0].view('complex')
+            corr_data[name].append(raw_data)
+            if mom_in is not None:
+                assert np.allclose(mom_in, np.array(str(file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['pIn'])[3:-2].strip().split(' '), dtype=int))
+            else:
+                mom_in = np.array(str(file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['pIn'])[3:-2].strip().split(' '), dtype=int)
+            if mom_out is not None:
+                assert np.allclose(mom_out, np.array(str(file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['pOut'])[3:-2].strip().split(' '), dtype=int))
+            else:
+                mom_out = np.array(str(file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['pOut'])[3:-2].strip().split(' '), dtype=int)
+
+        file.close()
+
+    result_dict = {}
+
+    for key, data in corr_data.items():
+        local_data = np.array(data)
+
+        rolled_array = np.rollaxis(local_data, 0, 5)
+
+        matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
+        for si, sj, ci, cj in np.ndindex(rolled_array.shape[:-1]):
+            real = Obs([rolled_array[si, sj, ci, cj].real], [ens_id])
+            imag = Obs([rolled_array[si, sj, ci, cj].imag], [ens_id])
+            matrix[si, sj, ci, cj] = CObs(real, imag)
+            matrix[si, sj, ci, cj].gamma_method()
+
+        result_dict[key] = Npr_matrix(matrix.swapaxes(1, 2).reshape((12, 12), order=order), mom_in=mom_in, mom_out=mom_out)
+
+    return result_dict

pyerrors/linalg.py

@@ -2,59 +2,174 @@
 # coding: utf-8
 
 import numpy as np
+from autograd import jacobian
 import autograd.numpy as anp  # Thinly-wrapped numpy
-from .pyerrors import derived_observable, CObs
+from .pyerrors import derived_observable, CObs, Obs
 
-
-# This code block is directly taken from the current master branch of autograd and remains
-# only until the new version is released on PyPi
 from functools import partial
 from autograd.extend import defvjp
 
-_dot = partial(anp.einsum, '...ij,...jk->...ik')
+
+def derived_array(func, data, **kwargs):
+    """Construct a derived Obs according to func(data, **kwargs) of matrix value data
+    using automatic differentiation.
+
+    Parameters
+    ----------
+    func -- arbitrary function of the form func(data, **kwargs). For the
+            automatic differentiation to work, all numpy functions have to have
+            the autograd wrapper (use 'import autograd.numpy as anp').
+    data -- list of Obs, e.g. [obs1, obs2, obs3].
+
+    Keyword arguments
+    -----------------
+    man_grad -- manually supply a list or an array which contains the jacobian
+                of func. Use cautiously, supplying the wrong derivative will
+                not be intercepted.
+    """
+
+    data = np.asarray(data)
+    raveled_data = data.ravel()
+
+    # Workaround for matrix operations containing non Obs data
+    for i_data in raveled_data:
+        if isinstance(i_data, Obs):
+            first_name = i_data.names[0]
+            first_shape = i_data.shape[first_name]
+            break
+
+    for i in range(len(raveled_data)):
+        if isinstance(raveled_data[i], (int, float)):
+            raveled_data[i] = Obs([raveled_data[i] + np.zeros(first_shape)], [first_name])
+
+    n_obs = len(raveled_data)
+    new_names = sorted(set([y for x in [o.names for o in raveled_data] for y in x]))
+
+    new_shape = {}
+    for i_data in raveled_data:
+        for name in new_names:
+            tmp = i_data.shape.get(name)
+            if tmp is not None:
+                if new_shape.get(name) is None:
+                    new_shape[name] = tmp
+                else:
+                    if new_shape[name] != tmp:
+                        raise Exception('Shapes of ensemble', name, 'do not match.')
+    if data.ndim == 1:
+        values = np.array([o.value for o in data])
+    else:
+        values = np.vectorize(lambda x: x.value)(data)
+
+    new_values = func(values, **kwargs)
+
+    new_r_values = {}
+    for name in new_names:
+        tmp_values = np.zeros(n_obs)
+        for i, item in enumerate(raveled_data):
+            tmp = item.r_values.get(name)
+            if tmp is None:
+                tmp = item.value
+            tmp_values[i] = tmp
+        tmp_values = np.array(tmp_values).reshape(data.shape)
+        new_r_values[name] = func(tmp_values, **kwargs)
+
+    if 'man_grad' in kwargs:
+        deriv = np.asarray(kwargs.get('man_grad'))
+        if new_values.shape + data.shape != deriv.shape:
+            raise Exception('Manual derivative does not have correct shape.')
+    elif kwargs.get('num_grad') is True:
+        raise Exception('Multi mode currently not supported for numerical derivative')
+    else:
+        deriv = jacobian(func)(values, **kwargs)
+
+    final_result = np.zeros(new_values.shape, dtype=object)
+
+    d_extracted = {}
+    for name in new_names:
+        d_extracted[name] = []
+        for i_dat, dat in enumerate(data):
+            ens_length = dat.ravel()[0].shape[name]
+            d_extracted[name].append(np.array([o.deltas[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):
+        new_deltas = {}
+        for name in new_names:
+            ens_length = d_extracted[name][0].shape[-1]
+            new_deltas[name] = np.zeros(ens_length)
+            for i_dat, dat in enumerate(d_extracted[name]):
+                new_deltas[name] += np.tensordot(deriv[i_val + (i_dat, )], dat)
+
+        new_samples = []
+        new_means = []
+        for name in new_names:
+            new_samples.append(new_deltas[name])
+            new_means.append(new_r_values[name][i_val])
+
+        final_result[i_val] = Obs(new_samples, new_names, means=new_means)
+        final_result[i_val]._value = new_val
+
+    return final_result
 
 
-# batched diag
-def _diag(a):
-    return anp.eye(a.shape[-1]) * a
+def matmul(*operands):
+    """Matrix multiply all operands.
+
+       Supports real and complex valued matrices and is faster compared to
+       standard multiplication via the @ operator.
+    """
+    if any(isinstance(o[0, 0], CObs) for o in operands):
+        extended_operands = []
+        for op in operands:
+            tmp = np.vectorize(lambda x: (np.real(x), np.imag(x)))(op)
+            extended_operands.append(tmp[0])
+            extended_operands.append(tmp[1])
+
+        def multi_dot(operands, part):
+            stack_r = operands[0]
+            stack_i = operands[1]
+            for op_r, op_i in zip(operands[2::2], operands[3::2]):
+                tmp_r = stack_r @ op_r - stack_i @ op_i
+                tmp_i = stack_r @ op_i + stack_i @ op_r
+
+                stack_r = tmp_r
+                stack_i = tmp_i
+
+            if part == 'Real':
+                return stack_r
+            else:
+                return stack_i
+
+        def multi_dot_r(operands):
+            return multi_dot(operands, 'Real')
+
+        def multi_dot_i(operands):
+            return multi_dot(operands, 'Imag')
+
+        Nr = derived_array(multi_dot_r, extended_operands)
+        Ni = derived_array(multi_dot_i, extended_operands)
+
+        res = np.empty_like(Nr)
+        for (n, m), entry in np.ndenumerate(Nr):
+            res[n, m] = CObs(Nr[n, m], Ni[n, m])
+
+        return res
+    else:
+        def multi_dot(operands):
+            stack = operands[0]
+            for op in operands[1:]:
+                stack = stack @ op
+            return stack
+        return derived_array(multi_dot, operands)
 
 
-# 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 inv(x):
+    """Inverse of Obs or CObs valued matrices."""
+    return _mat_mat_op(anp.linalg.inv, x)
 
 
-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
+def cholesky(x):
+    """Cholesky decompostion of Obs or CObs valued matrices."""
+    return _mat_mat_op(anp.linalg.cholesky, x)
 
 
 def scalar_mat_op(op, obs, **kwargs):
@@ -82,7 +197,7 @@ def scalar_mat_op(op, obs, **kwargs):
     return derived_observable(_mat, raveled_obs, **kwargs)
 
 
-def mat_mat_op(op, obs, **kwargs):
+def _mat_mat_op(op, obs, **kwargs):
     """Computes the matrix to matrix operation op to a given matrix of Obs."""
     # Use real representation to calculate matrix operations for complex matrices
     if isinstance(obs.ravel()[0], CObs):
@@ -99,15 +214,18 @@ def mat_mat_op(op, obs, **kwargs):
         if kwargs.get('num_grad') is True:
             op_big_matrix = _num_diff_mat_mat_op(op, big_matrix, **kwargs)
         else:
-            op_big_matrix = derived_observable(lambda x, **kwargs: op(x), big_matrix)
+            op_big_matrix = derived_array(lambda x, **kwargs: op(x), [big_matrix])[0]
         dim = op_big_matrix.shape[0]
         op_A = op_big_matrix[0: dim // 2, 0: dim // 2]
         op_B = op_big_matrix[dim // 2:, 0: dim // 2]
-        return (1 + 0j) * op_A + 1j * op_B
+        res = np.empty_like(op_A)
+        for (n, m), entry in np.ndenumerate(op_A):
+            res[n, m] = CObs(op_A[n, m], op_B[n, m])
+        return res
     else:
         if kwargs.get('num_grad') is True:
             return _num_diff_mat_mat_op(op, obs, **kwargs)
-        return derived_observable(lambda x, **kwargs: op(x), obs)
+        return derived_array(lambda x, **kwargs: op(x), [obs])[0]
 
 
 def eigh(obs, **kwargs):
@@ -146,7 +264,7 @@ def svd(obs, **kwargs):
     return (u, s, vh)
 
 
-def slog_det(obs, **kwargs):
+def slogdet(obs, **kwargs):
     """Computes the determinant of a matrix of Obs via np.linalg.slogdet."""
     def _mat(x):
         dim = int(np.sqrt(len(x)))
@@ -375,3 +493,51 @@ 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

pyerrors/npr.py

@@ -1,7 +1,6 @@
 import warnings
 import numpy as np
-import autograd.numpy as anp
-from .linalg import mat_mat_op
+from .linalg import inv, matmul
 from .dirac import gamma, gamma5
 
 
@@ -27,10 +26,9 @@ class Npr_matrix(np.ndarray):
         if self.shape != (12, 12):
             raise Exception('g5H only works for 12x12 matrices.')
         extended_g5 = np.kron(np.eye(3, dtype=int), gamma5)
-        new_matrix = extended_g5 @ self.conj().T @ extended_g5
-        new_matrix.mom_in = self.mom_out
-        new_matrix.mom_out = self.mom_in
-        return new_matrix
+        return Npr_matrix(matmul(extended_g5, self.conj().T, extended_g5),
+                          mom_in=self.mom_out,
+                          mom_out=self.mom_in)
 
     def _propagate_mom(self, other, name):
         s_mom = getattr(self, name, None)
@@ -69,7 +67,7 @@ def inv_propagator(prop):
     """ Inverts a 12x12 quark propagator"""
     if prop.shape != (12, 12):
         raise Exception("Only 12x12 propagators can be inverted.")
-    return Npr_matrix(mat_mat_op(anp.linalg.inv, prop), prop.mom_in)
+    return Npr_matrix(inv(prop), prop.mom_in)
 
 
 def Zq(inv_prop, fermion='Wilson'):
@@ -87,10 +85,18 @@ def Zq(inv_prop, fermion='Wilson'):
 
     if fermion == 'Wilson':
         p_slash = -1j * (sin_mom[0] * gamma[0] + sin_mom[1] * gamma[1] + sin_mom[2] * gamma[2] + sin_mom[3] * gamma[3]) / np.sum(sin_mom ** 2)
+    elif fermion == 'Continuum':
+        p_mom = 2 * np.pi / L * mom
+        p_slash = -1j * (p_mom[0] * gamma[0] + p_mom[1] * gamma[1] + p_mom[2] * gamma[2] + p_mom[3] * gamma[3]) / np.sum(p_mom ** 2)
+    elif fermion == 'DWF':
+        W = np.sum(1 - np.cos(2 * np.pi / L * mom))
+        s2 = np.sum(sin_mom ** 2)
+        p_slash = -1j * (sin_mom[0] * gamma[0] + sin_mom[1] * gamma[1] + sin_mom[2] * gamma[2] + sin_mom[3] * gamma[3])
+        p_slash /= 2 * (W - 1 + np.sqrt((1 - W) ** 2 + s2))
     else:
         raise Exception("Fermion type '" + fermion + "' not implemented")
 
-    res = 1 / 12. * np.trace(inv_prop @ np.kron(np.eye(3, dtype=int), p_slash))
+    res = 1 / 12. * np.trace(matmul(inv_prop, np.kron(np.eye(3, dtype=int), p_slash)))
     res.gamma_method()
 
     if not res.imag.is_zero_within_error(5):

pyerrors/pyerrors.py

@@ -34,11 +34,11 @@ class Obs:
                     ensemble.
     N_sigma_global -- Standard value for N_sigma (default 1.0)
     """
-    # __slots__ = ['names', 'shape', 'r_values', 'deltas', 'N', 'value', 'dvalue',
-    #              'ddvalue', 'reweighted', 'S', 'tau_exp', 'N_sigma', 'e_names',
-    #              'e_content', 'e_dvalue', 'e_ddvalue', 'e_tauint', 'e_dtauint',
-    #              'e_windowsize', 'e_rho', 'e_drho', 'e_n_tauint', 'e_n_dtauint',
-    #              'tag']
+    __slots__ = ['names', 'shape', 'r_values', 'deltas', 'N', '_value', '_dvalue',
+                 'ddvalue', 'reweighted', 'S', 'tau_exp', 'N_sigma', 'e_names',
+                 'e_content', 'e_dvalue', 'e_ddvalue', 'e_tauint', 'e_dtauint',
+                 'e_windowsize', 'e_rho', 'e_drho', 'e_n_tauint', 'e_n_dtauint',
+                 'tag', '__dict__']
 
     e_tag_global = 0
     S_global = 2.0
@@ -111,23 +111,6 @@ class Obs:
         self.ddvalue = 0.0
         self.reweighted = 0
 
-        self.S = {}
-        self.tau_exp = {}
-        self.N_sigma = 0
-
-        self.e_names = {}
-        self.e_content = {}
-
-        self.e_dvalue = {}
-        self.e_ddvalue = {}
-        self.e_tauint = {}
-        self.e_dtauint = {}
-        self.e_windowsize = {}
-        self.e_rho = {}
-        self.e_drho = {}
-        self.e_n_tauint = {}
-        self.e_n_dtauint = {}
-
         self.tag = None
 
     @property
@@ -392,33 +375,35 @@ class Obs:
             else:
                 percentage = np.abs(self.dvalue / self.value) * 100
             print('Result\t %3.8e +/- %3.8e +/- %3.8e (%3.3f%%)' % (self.value, self.dvalue, self.ddvalue, percentage))
-            if len(self.e_names) > 1:
-                print(' Ensemble errors:')
-            for e_name in self.e_names:
+            if hasattr(self, 'e_names'):
                 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))
-                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 level > 1:
-                print(self.N, 'samples in', len(self.e_names), 'ensembles:')
+                    print(' Ensemble errors:')
                 for e_name in self.e_names:
-                    print(e_name, ':', self.e_content[e_name])
+                    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))
+                    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 level > 1:
+                    print(self.N, 'samples in', len(self.e_names), 'ensembles:')
+                    for e_name in self.e_names:
+                        print(e_name, ':', self.e_content[e_name])
 
     def is_zero_within_error(self, sigma=1):
-        """ Checks whether the observable is zero within 'sigma' standard errors.
+        """Checks whether the observable is zero within 'sigma' standard errors.
 
         Works only properly when the gamma method was run.
         """
-        return np.abs(self.value) <= sigma * self.dvalue
+        return self.is_zero() or np.abs(self.value) <= sigma * self.dvalue
 
     def is_zero(self):
+        """Checks whether the observable is zero within machine precision."""
         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."""
-        if not self.e_names:
+        if not hasattr(self, 'e_names'):
             raise Exception('Run the gamma method first.')
 
         fig = plt.figure()
@@ -451,7 +436,7 @@ class Obs:
 
     def plot_rho(self):
         """Plot normalized autocorrelation function time for each ensemble."""
-        if not self.e_names:
+        if not hasattr(self, 'e_names'):
             raise Exception('Run the gamma method first.')
         for e, e_name in enumerate(self.e_names):
             plt.xlabel('W')
@@ -473,7 +458,7 @@ class Obs:
 
     def plot_rep_dist(self):
         """Plot replica distribution for each ensemble with more than one replicum."""
-        if not self.e_names:
+        if not hasattr(self, 'e_names'):
             raise Exception('Run the gamma method first.')
         for e, e_name in enumerate(self.e_names):
             if len(self.e_content[e_name]) == 1:
@@ -495,7 +480,7 @@ class Obs:
 
     def plot_history(self, expand=True):
         """Plot derived Monte Carlo history for each ensemble."""
-        if not self.e_names:
+        if not hasattr(self, 'e_names'):
             raise Exception('Run the gamma method first.')
 
         for e, e_name in enumerate(self.e_names):
@@ -519,7 +504,7 @@ class Obs:
     def plot_piechart(self):
         """Plot piechart which shows the fractional contribution of each
         ensemble to the error and returns a dictionary containing the fractions."""
-        if not self.e_names:
+        if not hasattr(self, 'e_names'):
             raise Exception('Run the gamma method first.')
         if self.dvalue == 0.0:
             raise Exception('Error is 0.0')
@@ -737,19 +722,23 @@ class CObs:
         return self._imag
 
     def gamma_method(self, **kwargs):
+        """Executes the gamma_method for the real and the imaginary part."""
         if isinstance(self.real, Obs):
             self.real.gamma_method(**kwargs)
         if isinstance(self.imag, Obs):
             self.imag.gamma_method(**kwargs)
 
     def is_zero(self):
+        """Checks whether both real and imaginary part are zero within machine precision."""
         return self.real == 0.0 and self.imag == 0.0
 
     def conjugate(self):
         return CObs(self.real, -self.imag)
 
     def __add__(self, other):
-        if hasattr(other, 'real') and hasattr(other, 'imag'):
+        if isinstance(other, np.ndarray):
+            return other + self
+        elif hasattr(other, 'real') and hasattr(other, 'imag'):
             return CObs(self.real + other.real,
                         self.imag + other.imag)
         else:
@@ -759,7 +748,9 @@ class CObs:
         return self + y
 
     def __sub__(self, other):
-        if hasattr(other, 'real') and hasattr(other, 'imag'):
+        if isinstance(other, np.ndarray):
+            return -1 * (other - self)
+        elif hasattr(other, 'real') and hasattr(other, 'imag'):
             return CObs(self.real - other.real, self.imag - other.imag)
         else:
             return CObs(self.real - other, self.imag)
@@ -768,29 +759,43 @@ class CObs:
         return -1 * (self - other)
 
     def __mul__(self, other):
-        if all(isinstance(i, Obs) for i in [self.real, self.imag, other.real, other.imag]):
-            return CObs(derived_observable(lambda x, **kwargs: x[0] * x[1] - x[2] * x[3],
-                                           [self.real, other.real, self.imag, other.imag],
-                                           man_grad=[other.real.value, self.real.value, -other.imag.value, -self.imag.value]),
-                        derived_observable(lambda x, **kwargs: x[2] * x[1] + x[0] * x[3],
-                                           [self.real, other.real, self.imag, other.imag],
-                                           man_grad=[other.imag.value, self.imag.value, other.real.value, self.real.value]))
-        elif hasattr(other, 'real') and getattr(other, 'imag', 0) != 0:
-            return CObs(self.real * other.real - self.imag * other.imag,
-                        self.imag * other.real + self.real * other.imag)
+        if isinstance(other, np.ndarray):
+            return other * self
+        elif hasattr(other, 'real') and hasattr(other, 'imag'):
+            if all(isinstance(i, Obs) for i in [self.real, self.imag, other.real, other.imag]):
+                return CObs(derived_observable(lambda x, **kwargs: x[0] * x[1] - x[2] * x[3],
+                                               [self.real, other.real, self.imag, other.imag],
+                                               man_grad=[other.real.value, self.real.value, -other.imag.value, -self.imag.value]),
+                            derived_observable(lambda x, **kwargs: x[2] * x[1] + x[0] * x[3],
+                                               [self.real, other.real, self.imag, other.imag],
+                                               man_grad=[other.imag.value, self.imag.value, other.real.value, self.real.value]))
+            elif getattr(other, 'imag', 0) != 0:
+                return CObs(self.real * other.real - self.imag * other.imag,
+                            self.imag * other.real + self.real * other.imag)
+            else:
+                return CObs(self.real * other.real, self.imag * other.real)
         else:
-            return CObs(self.real * np.real(other), self.imag * np.real(other))
+            return CObs(self.real * other, self.imag * other)
 
     def __rmul__(self, other):
         return self * other
 
     def __truediv__(self, other):
-        if hasattr(other, 'real') and hasattr(other, 'imag'):
+        if isinstance(other, np.ndarray):
+            return 1 / (other / self)
+        elif hasattr(other, 'real') and hasattr(other, 'imag'):
             r = other.real ** 2 + other.imag ** 2
             return CObs((self.real * other.real + self.imag * other.imag) / r, (self.imag * other.real - self.real * other.imag) / r)
         else:
             return CObs(self.real / other, self.imag / other)
 
+    def __rtruediv__(self, other):
+        r = self.real ** 2 + self.imag ** 2
+        if hasattr(other, 'real') and hasattr(other, 'imag'):
+            return CObs((self.real * other.real + self.imag * other.imag) / r, (self.real * other.imag - self.imag * other.real) / r)
+        else:
+            return CObs(self.real * other / r, -self.imag * other / r)
+
     def __abs__(self):
         return np.sqrt(self.real**2 + self.imag**2)
 
@@ -1148,7 +1153,7 @@ def covariance(obs1, obs2, correlation=False, **kwargs):
                 (1 != len(set([len(idx) for idx in [obs1.idl[name], obs2.idl[name], merge_idx([obs1.idl[name], obs2.idl[name]])]])))):
             raise Exception('Shapes of ensemble', name, 'do not fit')
 
-    if obs1.e_names == {} or obs2.e_names == {}:
+    if not hasattr(obs1, 'e_names') or not hasattr(obs2, 'e_names'):
         raise Exception('The gamma method has to be applied to both Obs first.')
 
     dvalue = 0
@@ -1232,7 +1237,7 @@ def covariance2(obs1, obs2, correlation=False, **kwargs):
 
         return gamma
 
-    if obs1.e_names == {} or obs2.e_names == {}:
+    if not hasattr(obs1, 'e_names') or not hasattr(obs2, 'e_names'):
         raise Exception('The gamma method has to be applied to both Obs first.')
 
     dvalue = 0
@@ -1322,7 +1327,7 @@ def covariance3(obs1, obs2, correlation=False, **kwargs):
                 (1 != len(set([len(idx) for idx in [obs1.idl[name], obs2.idl[name], merge_idx([obs1.idl[name], obs2.idl[name]])]])))):
             raise Exception('Shapes of ensemble', name, 'do not fit')
 
-    if obs1.e_names == {} or obs2.e_names == {}:
+    if not hasattr(obs1, 'e_names') or not hasattr(obs2, 'e_names'):
         raise Exception('The gamma method has to be applied to both Obs first.')
 
     tau_exp = []