pyerrors.correlators

   1import itertools
   2import warnings
   3from itertools import permutations
   4
   5import autograd.numpy as anp
   6import matplotlib.pyplot as plt
   7import numpy as np
   8import scipy.linalg
   9
  10from . import linalg
  11from .fits import least_squares
  12from .misc import _assert_equal_properties, dump_object
  13from .obs import CObs, Obs, correlate, reweight
  14from .roots import find_root
  15
  16
  17class Corr:
  18    r"""The class for a correlator (time dependent sequence of pe.Obs).
  19
  20    Everything, this class does, can be achieved using lists or arrays of Obs.
  21    But it is simply more convenient to have a dedicated object for correlators.
  22    One often wants to add or multiply correlators of the same length at every timeslice and it is inconvenient
  23    to iterate over all timeslices for every operation. This is especially true, when dealing with matrices.
  24
  25    The correlator can have two types of content: An Obs at every timeslice OR a matrix at every timeslice.
  26    Other dependency (eg. spatial) are not supported.
  27
  28    The Corr class can also deal with missing measurements or paddings for fixed boundary conditions.
  29    The missing entries are represented via the `None` object.
  30
  31    Initialization
  32    --------------
  33    A simple correlator can be initialized with a list or a one-dimensional array of `Obs` or `Cobs`
  34    ```python
  35    corr11 = pe.Corr([obs1, obs2])
  36    corr11 = pe.Corr(np.array([obs1, obs2]))
  37    ```
  38    A matrix-valued correlator can either be initialized via a two-dimensional array of `Corr` objects
  39    ```python
  40    matrix_corr = pe.Corr(np.array([[corr11, corr12], [corr21, corr22]]))
  41    ```
  42    or alternatively via a three-dimensional array of `Obs` or `CObs` of shape (T, N, N) where T is
  43    the temporal extent of the correlator and N is the dimension of the matrix.
  44    """
  45
  46    __slots__ = ["N", "T", "content", "prange", "tag"]
  47
  48    def __init__(self, data_input, padding=None, prange=None):
  49        """ Initialize a Corr object.
  50
  51        Parameters
  52        ----------
  53        data_input : list or array
  54            list of Obs or list of arrays of Obs or array of Corrs (see class docstring for details).
  55        padding : list, optional
  56            List with two entries where the first labels the padding
  57            at the front of the correlator and the second the padding
  58            at the back.
  59        prange : list, optional
  60            List containing the first and last timeslice of the plateau
  61            region identified for this correlator.
  62        """
  63
  64        if padding is None:
  65            padding = [0, 0]
  66
  67        if isinstance(data_input, np.ndarray):
  68            if data_input.ndim == 1:
  69                data_input = list(data_input)
  70            elif data_input.ndim == 2:
  71                if not data_input.shape[0] == data_input.shape[1]:
  72                    raise ValueError("Array needs to be square.")
  73                if not all([isinstance(item, Corr) for item in data_input.flatten()]):
  74                    raise ValueError("If the input is an array, its elements must be of type pe.Corr.")
  75                if not all([item.N == 1 for item in data_input.flatten()]):
  76                    raise ValueError("Can only construct matrix correlator from single valued correlators.")
  77                if not len(set([item.T for item in data_input.flatten()])) == 1:
  78                    raise ValueError("All input Correlators must be defined over the same timeslices.")
  79
  80                T = data_input[0, 0].T
  81                N = data_input.shape[0]
  82                input_as_list = []
  83                for t in range(T):
  84                    if any([(item.content[t] is None) for item in data_input.flatten()]):
  85                        if not all([(item.content[t] is None) for item in data_input.flatten()]):
  86                            warnings.warn("Input ill-defined at different timeslices. Conversion leads to data loss.!", RuntimeWarning, stacklevel=2)
  87                        input_as_list.append(None)
  88                    else:
  89                        array_at_timeslace = np.empty([N, N], dtype="object")
  90                        for i in range(N):
  91                            for j in range(N):
  92                                array_at_timeslace[i, j] = data_input[i, j][t]
  93                        input_as_list.append(array_at_timeslace)
  94                data_input = input_as_list
  95            elif data_input.ndim == 3:
  96                if not data_input.shape[1] == data_input.shape[2]:
  97                    raise ValueError("Array needs to be square.")
  98                data_input = list(data_input)
  99            else:
 100                raise ValueError("Arrays with ndim>3 not supported.")
 101
 102        if isinstance(data_input, list):
 103
 104            if all([isinstance(item, (Obs, CObs)) or item is None for item in data_input]):
 105                _assert_equal_properties([o for o in data_input if o is not None])
 106                self.content = [np.asarray([item]) if item is not None else None for item in data_input]
 107                self.N = 1
 108            elif all([isinstance(item, np.ndarray) or item is None for item in data_input]) and any([isinstance(item, np.ndarray) for item in data_input]):
 109                self.content = data_input
 110                noNull = [a for a in self.content if a is not None]  # To check if the matrices are correct for all undefined elements
 111                self.N = noNull[0].shape[0]
 112                if self.N > 1 and noNull[0].shape[0] != noNull[0].shape[1]:
 113                    raise ValueError("Smearing matrices are not NxN.")
 114                if (not all([item.shape == noNull[0].shape for item in noNull])):
 115                    raise ValueError("Items in data_input are not of identical shape." + str(noNull))
 116            else:
 117                raise TypeError("'data_input' contains item of wrong type.")
 118        else:
 119            raise TypeError("Data input was not given as list or correct array.")
 120
 121        self.tag = None
 122
 123        # An undefined timeslice is represented by the None object
 124        self.content = [None] * padding[0] + self.content + [None] * padding[1]
 125        self.T = len(self.content)
 126        self.prange = prange
 127
 128    def __getitem__(self, idx):
 129        """Return the content of timeslice idx"""
 130        if self.content[idx] is None:
 131            return None
 132        elif len(self.content[idx]) == 1:
 133            return self.content[idx][0]
 134        else:
 135            return self.content[idx]
 136
 137    @property
 138    def reweighted(self):
 139        bool_array = np.array([list(map(lambda x: x.reweighted, o)) for o in [x for x in self.content if x is not None]])
 140        if np.all(bool_array == 1):
 141            return True
 142        elif np.all(bool_array == 0):
 143            return False
 144        else:
 145            raise Exception("Reweighting status of correlator corrupted.")
 146
 147    def gamma_method(self, **kwargs):
 148        """Apply the gamma method to the content of the Corr."""
 149        for item in self.content:
 150            if item is not None:
 151                if self.N == 1:
 152                    item[0].gamma_method(**kwargs)
 153                else:
 154                    for i in range(self.N):
 155                        for j in range(self.N):
 156                            item[i, j].gamma_method(**kwargs)
 157
 158    gm = gamma_method
 159
 160    def projected(self, vector_l=None, vector_r=None, normalize=False):
 161        """We need to project the Correlator with a Vector to get a single value at each timeslice.
 162
 163        The method can use one or two vectors.
 164        If two are specified it returns v1@G@v2 (the order might be very important.)
 165        By default it will return the lowest source, which usually means unsmeared-unsmeared (0,0), but it does not have to
 166        """
 167        if self.N == 1:
 168            raise ValueError("Trying to project a Corr, that already has N=1.")
 169
 170        if vector_l is None:
 171            vector_l, vector_r = np.asarray([1.] + (self.N - 1) * [0.]), np.asarray([1.] + (self.N - 1) * [0.])
 172        elif (vector_r is None):
 173            vector_r = vector_l
 174        if isinstance(vector_l, list) and not isinstance(vector_r, list):
 175            if len(vector_l) != self.T:
 176                raise ValueError("Length of vector list must be equal to T")
 177            vector_r = [vector_r] * self.T
 178        if isinstance(vector_r, list) and not isinstance(vector_l, list):
 179            if len(vector_r) != self.T:
 180                raise ValueError("Length of vector list must be equal to T")
 181            vector_l = [vector_l] * self.T
 182
 183        if not isinstance(vector_l, list):
 184            if not vector_l.shape == vector_r.shape == (self.N,):
 185                raise ValueError("Vectors are of wrong shape!")
 186            if normalize:
 187                vector_l, vector_r = vector_l / np.sqrt(vector_l @ vector_l), vector_r / np.sqrt(vector_r @ vector_r)
 188            newcontent = [None if _check_for_none(self, item) else np.asarray([vector_l.T @ item @ vector_r]) for item in self.content]
 189
 190        else:
 191            # There are no checks here yet. There are so many possible scenarios, where this can go wrong.
 192            if normalize:
 193                for t in range(self.T):
 194                    vector_l[t], vector_r[t] = vector_l[t] / np.sqrt(vector_l[t] @ vector_l[t]), vector_r[t] / np.sqrt(vector_r[t] @ vector_r[t])
 195
 196            newcontent = [None if (_check_for_none(self, self.content[t]) or vector_l[t] is None or vector_r[t] is None) else np.asarray([vector_l[t].T @ self.content[t] @ vector_r[t]]) for t in range(self.T)]
 197        return Corr(newcontent)
 198
 199    def item(self, i, j):
 200        """Picks the element [i,j] from every matrix and returns a correlator containing one Obs per timeslice.
 201
 202        Parameters
 203        ----------
 204        i : int
 205            First index to be picked.
 206        j : int
 207            Second index to be picked.
 208        """
 209        if self.N == 1:
 210            raise ValueError("Trying to pick item from projected Corr")
 211        newcontent = [None if (item is None) else item[i, j] for item in self.content]
 212        return Corr(newcontent)
 213
 214    def plottable(self):
 215        """Outputs the correlator in a plotable format.
 216
 217        Outputs three lists containing the timeslice index, the value on each
 218        timeslice and the error on each timeslice.
 219        """
 220        if self.N != 1:
 221            raise ValueError("Can only make Corr[N=1] plottable")
 222        x_list = [x for x in range(self.T) if self.content[x] is not None]
 223        y_list = [y[0].value for y in self.content if y is not None]
 224        y_err_list = [y[0].dvalue for y in self.content if y is not None]
 225
 226        return x_list, y_list, y_err_list
 227
 228    def symmetric(self):
 229        """ Symmetrize the correlator around x0=0."""
 230        if self.N != 1:
 231            raise ValueError('symmetric cannot be safely applied to multi-dimensional correlators.')
 232        if self.T % 2 != 0:
 233            raise ValueError("Can not symmetrize odd T")
 234
 235        if self.content[0] is not None:
 236            if np.argmax(np.abs([o[0].value if o is not None else 0 for o in self.content])) != 0:
 237                warnings.warn("Correlator does not seem to be symmetric around x0=0.", RuntimeWarning, stacklevel=2)
 238
 239        newcontent = [self.content[0]]
 240        for t in range(1, self.T):
 241            if (self.content[t] is None) or (self.content[self.T - t] is None):
 242                newcontent.append(None)
 243            else:
 244                newcontent.append(0.5 * (self.content[t] + self.content[self.T - t]))
 245        if (all([x is None for x in newcontent])):
 246            raise ValueError("Corr could not be symmetrized: No redundant values")
 247        return Corr(newcontent, prange=self.prange)
 248
 249    def anti_symmetric(self):
 250        """Anti-symmetrize the correlator around x0=0."""
 251        if self.N != 1:
 252            raise TypeError('anti_symmetric cannot be safely applied to multi-dimensional correlators.')
 253        if self.T % 2 != 0:
 254            raise ValueError("Can not symmetrize odd T")
 255
 256        test = 1 * self
 257        test.gamma_method()
 258        if not all([o.is_zero_within_error(3) for o in test.content[0]]):
 259            warnings.warn("Correlator does not seem to be anti-symmetric around x0=0.", RuntimeWarning, stacklevel=2)
 260
 261        newcontent = [self.content[0]]
 262        for t in range(1, self.T):
 263            if (self.content[t] is None) or (self.content[self.T - t] is None):
 264                newcontent.append(None)
 265            else:
 266                newcontent.append(0.5 * (self.content[t] - self.content[self.T - t]))
 267        if (all([x is None for x in newcontent])):
 268            raise ValueError("Corr could not be symmetrized: No redundant values")
 269        return Corr(newcontent, prange=self.prange)
 270
 271    def is_matrix_symmetric(self):
 272        """Checks whether a correlator matrices is symmetric on every timeslice."""
 273        if self.N == 1:
 274            raise TypeError("Only works for correlator matrices.")
 275        for t in range(self.T):
 276            if self[t] is None:
 277                continue
 278            for i in range(self.N):
 279                for j in range(i + 1, self.N):
 280                    if self[t][i, j] is self[t][j, i]:
 281                        continue
 282                    if hash(self[t][i, j]) != hash(self[t][j, i]):
 283                        return False
 284        return True
 285
 286    def trace(self):
 287        """Calculates the per-timeslice trace of a correlator matrix."""
 288        if self.N == 1:
 289            raise ValueError("Only works for correlator matrices.")
 290        newcontent = []
 291        for t in range(self.T):
 292            if _check_for_none(self, self.content[t]):
 293                newcontent.append(None)
 294            else:
 295                newcontent.append(np.trace(self.content[t]))
 296        return Corr(newcontent)
 297
 298    def matrix_symmetric(self):
 299        """Symmetrizes the correlator matrices on every timeslice."""
 300        if self.N == 1:
 301            raise ValueError("Trying to symmetrize a correlator matrix, that already has N=1.")
 302        if self.is_matrix_symmetric():
 303            return 1.0 * self
 304        else:
 305            transposed = [None if _check_for_none(self, G) else G.T for G in self.content]
 306            return 0.5 * (Corr(transposed) + self)
 307
 308    def GEVP(self, t0, ts=None, sort="Eigenvalue", vector_obs=False, **kwargs):
 309        r'''Solve the generalized eigenvalue problem on the correlator matrix and returns the corresponding eigenvectors.
 310
 311        The eigenvectors are sorted according to the descending eigenvalues, the zeroth eigenvector(s) correspond to the
 312        largest eigenvalue(s). The eigenvector(s) for the individual states can be accessed via slicing
 313        ```python
 314        C.GEVP(t0=2)[0]  # Ground state vector(s)
 315        C.GEVP(t0=2)[:3]  # Vectors for the lowest three states
 316        ```
 317
 318        Parameters
 319        ----------
 320        t0 : int
 321            The time t0 for the right hand side of the GEVP according to $G(t)v_i=\lambda_i G(t_0)v_i$
 322        ts : int
 323            fixed time $G(t_s)v_i=\lambda_i G(t_0)v_i$ if sort=None.
 324            If sort="Eigenvector" it gives a reference point for the sorting method.
 325        sort : string
 326            If this argument is set, a list of self.T vectors per state is returned. If it is set to None, only one vector is returned.
 327            - "Eigenvalue": The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice. (default)
 328            - "Eigenvector": Use the method described in arXiv:2004.10472 to find the set of v(t) belonging to the state.
 329              The reference state is identified by its eigenvalue at $t=t_s$.
 330            - None: The GEVP is solved only at ts, no sorting is necessary
 331        vector_obs : bool
 332            If True, uncertainties are propagated in the eigenvector computation (default False).
 333
 334        Other Parameters
 335        ----------------
 336        state : int
 337           Returns only the vector(s) for a specified state. The lowest state is zero.
 338        method : str
 339           Method used to solve the GEVP.
 340           - "eigh": Use scipy.linalg.eigh to solve the GEVP. (default for vector_obs=False)
 341           - "cholesky": Use manually implemented solution via the Cholesky decomposition. Automatically chosen if vector_obs==True.
 342        '''
 343
 344        if self.N == 1:
 345            raise ValueError("GEVP methods only works on correlator matrices and not single correlators.")
 346        if ts is not None:
 347            if (ts <= t0):
 348                raise ValueError("ts has to be larger than t0.")
 349
 350        if "sorted_list" in kwargs:
 351            warnings.warn("Argument 'sorted_list' is deprecated, use 'sort' instead.", DeprecationWarning, stacklevel=2)
 352            sort = kwargs.get("sorted_list")
 353
 354        if self.is_matrix_symmetric():
 355            symmetric_corr = self
 356        else:
 357            symmetric_corr = self.matrix_symmetric()
 358
 359        def _get_mat_at_t(t, vector_obs=vector_obs):
 360            if vector_obs:
 361                return symmetric_corr[t]
 362            else:
 363                return np.vectorize(lambda x: x.value)(symmetric_corr[t])
 364        G0 = _get_mat_at_t(t0)
 365
 366        method = kwargs.get('method', 'eigh')
 367        if vector_obs:
 368            chol = linalg.cholesky(G0)
 369            chol_inv = linalg.inv(chol)
 370            method = 'cholesky'
 371        else:
 372            chol = np.linalg.cholesky(_get_mat_at_t(t0, vector_obs=False))  # Check if matrix G0 is positive-semidefinite.
 373            if method == 'cholesky':
 374                chol_inv = np.linalg.inv(chol)
 375            else:
 376                chol_inv = None
 377
 378        if sort is None:
 379            if (ts is None):
 380                raise ValueError("ts is required if sort=None.")
 381            if (self.content[t0] is None) or (self.content[ts] is None):
 382                raise ValueError("Corr not defined at t0/ts.")
 383            Gt = _get_mat_at_t(ts)
 384            reordered_vecs = _GEVP_solver(Gt, G0, method=method, chol_inv=chol_inv)
 385            if kwargs.get('auto_gamma', False) and vector_obs:
 386                [[o.gm() for o in ev if isinstance(o, Obs)] for ev in reordered_vecs]
 387
 388        elif sort in ["Eigenvalue", "Eigenvector"]:
 389            if sort == "Eigenvalue" and ts is not None:
 390                warnings.warn("ts has no effect when sorting by eigenvalue is chosen.", RuntimeWarning, stacklevel=2)
 391            all_vecs = [None] * (t0 + 1)
 392            for t in range(t0 + 1, self.T):
 393                try:
 394                    Gt = _get_mat_at_t(t)
 395                    all_vecs.append(_GEVP_solver(Gt, G0, method=method, chol_inv=chol_inv))
 396                except Exception:
 397                    all_vecs.append(None)
 398            if sort == "Eigenvector":
 399                if ts is None:
 400                    raise ValueError("ts is required for the Eigenvector sorting method.")
 401                all_vecs = _sort_vectors(all_vecs, ts)
 402
 403            reordered_vecs = [[v[s] if v is not None else None for v in all_vecs] for s in range(self.N)]
 404            if kwargs.get('auto_gamma', False) and vector_obs:
 405                [[[o.gm() for o in evn] for evn in ev if evn is not None] for ev in reordered_vecs]
 406        else:
 407            raise ValueError("Unknown value for 'sort'. Choose 'Eigenvalue', 'Eigenvector' or None.")
 408
 409        if "state" in kwargs:
 410            return reordered_vecs[kwargs.get("state")]
 411        else:
 412            return reordered_vecs
 413
 414    def Eigenvalue(self, t0, ts=None, state=0, sort="Eigenvalue", **kwargs):
 415        """Determines the eigenvalue of the GEVP by solving and projecting the correlator
 416
 417        Parameters
 418        ----------
 419        state : int
 420            The state one is interested in ordered by energy. The lowest state is zero.
 421
 422        All other parameters are identical to the ones of Corr.GEVP.
 423        """
 424        vec = self.GEVP(t0, ts=ts, sort=sort, **kwargs)[state]
 425        return self.projected(vec)
 426
 427    def Hankel(self, N, periodic=False):
 428        """Constructs an NxN Hankel matrix
 429
 430        C(t) c(t+1) ... c(t+n-1)
 431        C(t+1) c(t+2) ... c(t+n)
 432        .................
 433        C(t+(n-1)) c(t+n) ... c(t+2(n-1))
 434
 435        Parameters
 436        ----------
 437        N : int
 438            Dimension of the Hankel matrix
 439        periodic : bool, optional
 440            determines whether the matrix is extended periodically
 441        """
 442
 443        if self.N != 1:
 444            raise NotImplementedError("Multi-operator Prony not implemented!")
 445
 446        array = np.empty([N, N], dtype="object")
 447        new_content = []
 448        for _t in range(self.T):
 449            new_content.append(array.copy())
 450
 451        def wrap(i):
 452            while i >= self.T:
 453                i -= self.T
 454            return i
 455
 456        for t in range(self.T):
 457            for i in range(N):
 458                for j in range(N):
 459                    if periodic:
 460                        new_content[t][i, j] = self.content[wrap(t + i + j)][0]
 461                    elif (t + i + j) >= self.T:
 462                        new_content[t] = None
 463                    else:
 464                        new_content[t][i, j] = self.content[t + i + j][0]
 465
 466        return Corr(new_content)
 467
 468    def roll(self, dt):
 469        """Periodically shift the correlator by dt timeslices
 470
 471        Parameters
 472        ----------
 473        dt : int
 474            number of timeslices
 475        """
 476        return Corr(list(np.roll(np.array(self.content, dtype=object), dt, axis=0)))
 477
 478    def reverse(self):
 479        """Reverse the time ordering of the Corr"""
 480        return Corr(self.content[:: -1])
 481
 482    def thin(self, spacing=2, offset=0):
 483        """Thin out a correlator to suppress correlations
 484
 485        Parameters
 486        ----------
 487        spacing : int
 488            Keep only every 'spacing'th entry of the correlator
 489        offset : int
 490            Offset the equal spacing
 491        """
 492        new_content = []
 493        for t in range(self.T):
 494            if (offset + t) % spacing != 0:
 495                new_content.append(None)
 496            else:
 497                new_content.append(self.content[t])
 498        return Corr(new_content)
 499
 500    def correlate(self, partner):
 501        """Correlate the correlator with another correlator or Obs
 502
 503        Parameters
 504        ----------
 505        partner : Obs or Corr
 506            partner to correlate the correlator with.
 507            Can either be an Obs which is correlated with all entries of the
 508            correlator or a Corr of same length.
 509        """
 510        if self.N != 1:
 511            raise ValueError("Only one-dimensional correlators can be safely correlated.")
 512        new_content = []
 513        for x0, t_slice in enumerate(self.content):
 514            if _check_for_none(self, t_slice):
 515                new_content.append(None)
 516            else:
 517                if isinstance(partner, Corr):
 518                    if _check_for_none(partner, partner.content[x0]):
 519                        new_content.append(None)
 520                    else:
 521                        new_content.append(np.array([correlate(o, partner.content[x0][0]) for o in t_slice]))
 522                elif isinstance(partner, Obs):  # Should this include CObs?
 523                    new_content.append(np.array([correlate(o, partner) for o in t_slice]))
 524                else:
 525                    raise TypeError("Can only correlate with an Obs or a Corr.")
 526
 527        return Corr(new_content)
 528
 529    def reweight(self, weight, **kwargs):
 530        """Reweight the correlator.
 531
 532        Parameters
 533        ----------
 534        weight : Obs
 535            Reweighting factor. An Observable that has to be defined on a superset of the
 536            configurations in obs[i].idl for all i.
 537        all_configs : bool
 538            if True, the reweighted observables are normalized by the average of
 539            the reweighting factor on all configurations in weight.idl and not
 540            on the configurations in obs[i].idl.
 541        """
 542        if self.N != 1:
 543            raise ValueError("Reweighting only implemented for one-dimensional correlators.")
 544        new_content = []
 545        for t_slice in self.content:
 546            if _check_for_none(self, t_slice):
 547                new_content.append(None)
 548            else:
 549                new_content.append(np.array(reweight(weight, t_slice, **kwargs)))
 550        return Corr(new_content)
 551
 552    def T_symmetry(self, partner, parity=+1):
 553        """Return the time symmetry average of the correlator and its partner
 554
 555        Parameters
 556        ----------
 557        partner : Corr
 558            Time symmetry partner of the Corr
 559        parity : int
 560            Parity quantum number of the correlator, can be +1 or -1
 561        """
 562        if self.N != 1:
 563            raise ValueError("T_symmetry only implemented for one-dimensional correlators.")
 564        if not isinstance(partner, Corr):
 565            raise TypeError("T partner has to be a Corr object.")
 566        if parity not in [+1, -1]:
 567            raise ValueError("Parity has to be +1 or -1.")
 568        T_partner = parity * partner.reverse()
 569
 570        t_slices = []
 571        test = (self - T_partner)
 572        test.gamma_method()
 573        for x0, t_slice in enumerate(test.content):
 574            if t_slice is not None:
 575                if not t_slice[0].is_zero_within_error(5):
 576                    t_slices.append(x0)
 577        if t_slices:
 578            warnings.warn("T symmetry partners do not agree within 5 sigma on time slices " + str(t_slices) + ".", RuntimeWarning, stacklevel=2)
 579
 580        return (self + T_partner) / 2
 581
 582    def deriv(self, variant="symmetric"):
 583        """Return the first derivative of the correlator with respect to x0.
 584
 585        Parameters
 586        ----------
 587        variant : str
 588            decides which definition of the finite differences derivative is used.
 589            Available choice: symmetric, forward, backward, improved, log, default: symmetric
 590        """
 591        if self.N != 1:
 592            raise ValueError("deriv only implemented for one-dimensional correlators.")
 593        if variant == "symmetric":
 594            newcontent = []
 595            for t in range(1, self.T - 1):
 596                if (self.content[t - 1] is None) or (self.content[t + 1] is None):
 597                    newcontent.append(None)
 598                else:
 599                    newcontent.append(0.5 * (self.content[t + 1] - self.content[t - 1]))
 600            if (all([x is None for x in newcontent])):
 601                raise ValueError('Derivative is undefined at all timeslices')
 602            return Corr(newcontent, padding=[1, 1])
 603        elif variant == "forward":
 604            newcontent = []
 605            for t in range(self.T - 1):
 606                if (self.content[t] is None) or (self.content[t + 1] is None):
 607                    newcontent.append(None)
 608                else:
 609                    newcontent.append(self.content[t + 1] - self.content[t])
 610            if (all([x is None for x in newcontent])):
 611                raise ValueError("Derivative is undefined at all timeslices")
 612            return Corr(newcontent, padding=[0, 1])
 613        elif variant == "backward":
 614            newcontent = []
 615            for t in range(1, self.T):
 616                if (self.content[t - 1] is None) or (self.content[t] is None):
 617                    newcontent.append(None)
 618                else:
 619                    newcontent.append(self.content[t] - self.content[t - 1])
 620            if (all([x is None for x in newcontent])):
 621                raise ValueError("Derivative is undefined at all timeslices")
 622            return Corr(newcontent, padding=[1, 0])
 623        elif variant == "improved":
 624            newcontent = []
 625            for t in range(2, self.T - 2):
 626                if (self.content[t - 2] is None) or (self.content[t - 1] is None) or (self.content[t + 1] is None) or (self.content[t + 2] is None):
 627                    newcontent.append(None)
 628                else:
 629                    newcontent.append((1 / 12) * (self.content[t - 2] - 8 * self.content[t - 1] + 8 * self.content[t + 1] - self.content[t + 2]))
 630            if (all([x is None for x in newcontent])):
 631                raise ValueError('Derivative is undefined at all timeslices')
 632            return Corr(newcontent, padding=[2, 2])
 633        elif variant == 'log':
 634            newcontent = []
 635            for t in range(self.T):
 636                if (self.content[t] is None) or (self.content[t] <= 0):
 637                    newcontent.append(None)
 638                else:
 639                    newcontent.append(np.log(self.content[t]))
 640            if (all([x is None for x in newcontent])):
 641                raise ValueError("Log is undefined at all timeslices")
 642            logcorr = Corr(newcontent)
 643            return self * logcorr.deriv('symmetric')
 644        else:
 645            raise ValueError("Unknown variant.")
 646
 647    def second_deriv(self, variant="symmetric"):
 648        r"""Return the second derivative of the correlator with respect to x0.
 649
 650        Parameters
 651        ----------
 652        variant : str
 653            decides which definition of the finite differences derivative is used.
 654            Available choice:
 655                - symmetric (default)
 656                    $$\tilde{\partial}^2_0 f(x_0) = f(x_0+1)-2f(x_0)+f(x_0-1)$$
 657                - big_symmetric
 658                    $$\partial^2_0 f(x_0) = \frac{f(x_0+2)-2f(x_0)+f(x_0-2)}{4}$$
 659                - improved
 660                    $$\partial^2_0 f(x_0) = \frac{-f(x_0+2) + 16 * f(x_0+1) - 30 * f(x_0) + 16 * f(x_0-1) - f(x_0-2)}{12}$$
 661                - log
 662                    $$f(x) = \tilde{\partial}^2_0 log(f(x_0))+(\tilde{\partial}_0 log(f(x_0)))^2$$
 663        """
 664        if self.N != 1:
 665            raise ValueError("second_deriv only implemented for one-dimensional correlators.")
 666        if variant == "symmetric":
 667            newcontent = []
 668            for t in range(1, self.T - 1):
 669                if (self.content[t - 1] is None) or (self.content[t + 1] is None):
 670                    newcontent.append(None)
 671                else:
 672                    newcontent.append(self.content[t + 1] - 2 * self.content[t] + self.content[t - 1])
 673            if (all([x is None for x in newcontent])):
 674                raise ValueError("Derivative is undefined at all timeslices")
 675            return Corr(newcontent, padding=[1, 1])
 676        elif variant == "big_symmetric":
 677            newcontent = []
 678            for t in range(2, self.T - 2):
 679                if (self.content[t - 2] is None) or (self.content[t + 2] is None):
 680                    newcontent.append(None)
 681                else:
 682                    newcontent.append((self.content[t + 2] - 2 * self.content[t] + self.content[t - 2]) / 4)
 683            if (all([x is None for x in newcontent])):
 684                raise ValueError("Derivative is undefined at all timeslices")
 685            return Corr(newcontent, padding=[2, 2])
 686        elif variant == "improved":
 687            newcontent = []
 688            for t in range(2, self.T - 2):
 689                if (self.content[t - 2] is None) or (self.content[t - 1] is None) or (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t + 2] is None):
 690                    newcontent.append(None)
 691                else:
 692                    newcontent.append((1 / 12) * (-self.content[t + 2] + 16 * self.content[t + 1] - 30 * self.content[t] + 16 * self.content[t - 1] - self.content[t - 2]))
 693            if (all([x is None for x in newcontent])):
 694                raise ValueError("Derivative is undefined at all timeslices")
 695            return Corr(newcontent, padding=[2, 2])
 696        elif variant == 'log':
 697            newcontent = []
 698            for t in range(self.T):
 699                if (self.content[t] is None) or (self.content[t] <= 0):
 700                    newcontent.append(None)
 701                else:
 702                    newcontent.append(np.log(self.content[t]))
 703            if (all([x is None for x in newcontent])):
 704                raise ValueError("Log is undefined at all timeslices")
 705            logcorr = Corr(newcontent)
 706            return self * (logcorr.second_deriv('symmetric') + (logcorr.deriv('symmetric'))**2)
 707        else:
 708            raise ValueError("Unknown variant.")
 709
 710    def m_eff(self, variant='log', guess=1.0):
 711        """Returns the effective mass of the correlator as correlator object
 712
 713        Parameters
 714        ----------
 715        variant : str
 716            log : uses the standard effective mass log(C(t) / C(t+1))
 717            cosh, periodic : Use periodicity of the correlator by solving C(t) / C(t+1) = cosh(m * (t - T/2)) / cosh(m * (t + 1 - T/2)) for m.
 718            sinh : Use anti-periodicity of the correlator by solving C(t) / C(t+1) = sinh(m * (t - T/2)) / sinh(m * (t + 1 - T/2)) for m.
 719            See, e.g., arXiv:1205.5380
 720            arccosh : Uses the explicit form of the symmetrized correlator (not recommended)
 721            logsym: uses the symmetric effective mass log(C(t-1) / C(t+1))/2
 722        guess : float
 723            guess for the root finder, only relevant for the root variant
 724        """
 725        if self.N != 1:
 726            raise ValueError('Correlator must be projected before getting m_eff')
 727        if variant == 'log':
 728            newcontent = []
 729            for t in range(self.T - 1):
 730                if ((self.content[t] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0):
 731                    newcontent.append(None)
 732                elif self.content[t][0].value / self.content[t + 1][0].value < 0:
 733                    newcontent.append(None)
 734                else:
 735                    newcontent.append(self.content[t] / self.content[t + 1])
 736            if (all([x is None for x in newcontent])):
 737                raise ValueError('m_eff is undefined at all timeslices')
 738
 739            return np.log(Corr(newcontent, padding=[0, 1]))
 740
 741        elif variant == 'logsym':
 742            newcontent = []
 743            for t in range(1, self.T - 1):
 744                if ((self.content[t - 1] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0):
 745                    newcontent.append(None)
 746                elif self.content[t - 1][0].value / self.content[t + 1][0].value < 0:
 747                    newcontent.append(None)
 748                else:
 749                    newcontent.append(self.content[t - 1] / self.content[t + 1])
 750            if (all([x is None for x in newcontent])):
 751                raise ValueError('m_eff is undefined at all timeslices')
 752
 753            return np.log(Corr(newcontent, padding=[1, 1])) / 2
 754
 755        elif variant in ['periodic', 'cosh', 'sinh']:
 756            if variant in ['periodic', 'cosh']:
 757                func = anp.cosh
 758            else:
 759                func = anp.sinh
 760
 761            def root_function(x, d):
 762                return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d
 763
 764            newcontent = []
 765            for t in range(self.T - 1):
 766                if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t + 1][0].value == 0):
 767                    newcontent.append(None)
 768                # Fill the two timeslices in the middle of the lattice with their predecessors
 769                elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]:
 770                    newcontent.append(newcontent[-1])
 771                elif self.content[t][0].value / self.content[t + 1][0].value < 0:
 772                    newcontent.append(None)
 773                else:
 774                    newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess)))
 775            if (all([x is None for x in newcontent])):
 776                raise ValueError('m_eff is undefined at all timeslices')
 777
 778            return Corr(newcontent, padding=[0, 1])
 779
 780        elif variant == 'arccosh':
 781            newcontent = []
 782            for t in range(1, self.T - 1):
 783                if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t - 1] is None) or (self.content[t][0].value == 0):
 784                    newcontent.append(None)
 785                else:
 786                    newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t]))
 787            if (all([x is None for x in newcontent])):
 788                raise ValueError("m_eff is undefined at all timeslices")
 789            return np.arccosh(Corr(newcontent, padding=[1, 1]))
 790
 791        else:
 792            raise ValueError('Unknown variant.')
 793
 794    def fit(self, function, fitrange=None, silent=False, **kwargs):
 795        r'''Fits function to the data
 796
 797        Parameters
 798        ----------
 799        function : obj
 800            function to fit to the data. See fits.least_squares for details.
 801        fitrange : list
 802            Two element list containing the timeslices on which the fit is supposed to start and stop.
 803            Caution: This range is inclusive as opposed to standard python indexing.
 804            `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6.
 805            If not specified, self.prange or all timeslices are used.
 806        silent : bool
 807            Decides whether output is printed to the standard output.
 808        '''
 809        if self.N != 1:
 810            raise ValueError("Correlator must be projected before fitting")
 811
 812        if fitrange is None:
 813            if self.prange:
 814                fitrange = self.prange
 815            else:
 816                fitrange = [0, self.T - 1]
 817        else:
 818            if not isinstance(fitrange, list):
 819                raise TypeError("fitrange has to be a list with two elements")
 820            if len(fitrange) != 2:
 821                raise ValueError("fitrange has to have exactly two elements [fit_start, fit_stop]")
 822
 823        xs = np.array([x for x in range(fitrange[0], fitrange[1] + 1) if self.content[x] is not None])
 824        ys = np.array([self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if self.content[x] is not None])
 825        result = least_squares(xs, ys, function, silent=silent, **kwargs)
 826        return result
 827
 828    def plateau(self, plateau_range=None, method="fit", auto_gamma=False):
 829        """ Extract a plateau value from a Corr object
 830
 831        Parameters
 832        ----------
 833        plateau_range : list
 834            list with two entries, indicating the first and the last timeslice
 835            of the plateau region.
 836        method : str
 837            method to extract the plateau.
 838                'fit' fits a constant to the plateau region
 839                'avg', 'average' or 'mean' just average over the given timeslices.
 840        auto_gamma : bool
 841            apply gamma_method with default parameters to the Corr. Defaults to None
 842        """
 843        if not plateau_range:
 844            if self.prange:
 845                plateau_range = self.prange
 846            else:
 847                raise ValueError("no plateau range provided")
 848        if self.N != 1:
 849            raise ValueError("Correlator must be projected before getting a plateau.")
 850        if (all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])):
 851            raise ValueError("plateau is undefined at all timeslices in plateaurange.")
 852        if auto_gamma:
 853            self.gamma_method()
 854        if method == "fit":
 855            def const_func(a, t):
 856                return a[0]
 857            return self.fit(const_func, plateau_range)[0]
 858        elif method in ["avg", "average", "mean"]:
 859            returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None])
 860            return returnvalue
 861
 862        else:
 863            raise ValueError("Unsupported plateau method: " + method)
 864
 865    def set_prange(self, prange):
 866        """Sets the attribute prange of the Corr object."""
 867        if not len(prange) == 2:
 868            raise ValueError("prange must be a list or array with two values")
 869        if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))):
 870            raise TypeError("Start and end point must be integers")
 871        if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] <= prange[1]):
 872            raise ValueError("Start and end point must define a range in the interval 0,T")
 873
 874        self.prange = prange
 875        return
 876
 877    def show(self, x_range=None, comp=None, y_range=None, logscale=False, plateau=None, fit_res=None, fit_key=None, ylabel=None, save=None, auto_gamma=False, hide_sigma=None, references=None, title=None):
 878        """Plots the correlator using the tag of the correlator as label if available.
 879
 880        Parameters
 881        ----------
 882        x_range : list
 883            list of two values, determining the range of the x-axis e.g. [4, 8].
 884        comp : Corr or list of Corr
 885            Correlator or list of correlators which are plotted for comparison.
 886            The tags of these correlators are used as labels if available.
 887        logscale : bool
 888            Sets y-axis to logscale.
 889        plateau : Obs
 890            Plateau value to be visualized in the figure.
 891        fit_res : Fit_result
 892            Fit_result object to be visualized.
 893        fit_key : str
 894            Key for the fit function in Fit_result.fit_function (for combined fits).
 895        ylabel : str
 896            Label for the y-axis.
 897        save : str
 898            path to file in which the figure should be saved.
 899        auto_gamma : bool
 900            Apply the gamma method with standard parameters to all correlators and plateau values before plotting.
 901        hide_sigma : float
 902            Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors.
 903        references : list
 904            List of floating point values that are displayed as horizontal lines for reference.
 905        title : string
 906            Optional title of the figure.
 907        """
 908        if self.N != 1:
 909            raise ValueError("Correlator must be projected before plotting")
 910
 911        if auto_gamma:
 912            self.gamma_method()
 913
 914        if x_range is None:
 915            x_range = [0, self.T - 1]
 916
 917        fig = plt.figure()
 918        ax1 = fig.add_subplot(111)
 919
 920        x, y, y_err = self.plottable()
 921        if hide_sigma:
 922            hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
 923        else:
 924            hide_from = None
 925        ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag)
 926        if logscale:
 927            ax1.set_yscale('log')
 928        else:
 929            if y_range is None:
 930                try:
 931                    y_min = min([(x[0].value - x[0].dvalue) for x in self.content[x_range[0]: x_range[1] + 1] if (x is not None) and x[0].dvalue < 2 * np.abs(x[0].value)])
 932                    y_max = max([(x[0].value + x[0].dvalue) for x in self.content[x_range[0]: x_range[1] + 1] if (x is not None) and x[0].dvalue < 2 * np.abs(x[0].value)])
 933                    ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)])
 934                except Exception:
 935                    pass
 936            else:
 937                ax1.set_ylim(y_range)
 938        if comp:
 939            if isinstance(comp, (Corr, list)):
 940                for corr in comp if isinstance(comp, list) else [comp]:
 941                    if auto_gamma:
 942                        corr.gamma_method()
 943                    x, y, y_err = corr.plottable()
 944                    if hide_sigma:
 945                        hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
 946                    else:
 947                        hide_from = None
 948                    ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor'])
 949            else:
 950                raise TypeError("'comp' must be a correlator or a list of correlators.")
 951
 952        if plateau:
 953            if isinstance(plateau, Obs):
 954                if auto_gamma:
 955                    plateau.gamma_method()
 956                ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau))
 957                ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-')
 958            else:
 959                raise TypeError("'plateau' must be an Obs")
 960
 961        if references:
 962            if isinstance(references, list):
 963                for ref in references:
 964                    ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--')
 965            else:
 966                raise TypeError("'references' must be a list of floating pint values.")
 967
 968        if self.prange:
 969            ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',', color="black", zorder=0)
 970            ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',', color="black", zorder=0)
 971
 972        if fit_res:
 973            x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05)
 974            if isinstance(fit_res.fit_function, dict):
 975                if fit_key:
 976                    ax1.plot(x_samples, fit_res.fit_function[fit_key]([o.value for o in fit_res.fit_parameters], x_samples), ls='-', marker=',', lw=2)
 977                else:
 978                    raise ValueError("Please provide a 'fit_key' for visualizing combined fits.")
 979            else:
 980                ax1.plot(x_samples, fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples), ls='-', marker=',', lw=2)
 981
 982        ax1.set_xlabel(r'$x_0 / a$')
 983        if ylabel:
 984            ax1.set_ylabel(ylabel)
 985        ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5])
 986
 987        _handles, labels = ax1.get_legend_handles_labels()
 988        if labels:
 989            ax1.legend()
 990
 991        if title:
 992            plt.title(title)
 993
 994        plt.draw()
 995
 996        if save:
 997            if isinstance(save, str):
 998                fig.savefig(save, bbox_inches='tight')
 999            else:
1000                raise TypeError("'save' has to be a string.")
1001
1002    def spaghetti_plot(self, logscale=True):
1003        """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations.
1004
1005        Parameters
1006        ----------
1007        logscale : bool
1008            Determines whether the scale of the y-axis is logarithmic or standard.
1009        """
1010        if self.N != 1:
1011            raise ValueError("Correlator needs to be projected first.")
1012
1013        mc_names = list(set([item for sublist in [list(itertools.chain.from_iterable(map(o[0].e_content.get, o[0].mc_names))) for o in self.content if o is not None] for item in sublist]))
1014        x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content, strict=True) if o is not None]
1015
1016        for name in mc_names:
1017            data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T
1018
1019            fig = plt.figure()
1020            ax = fig.add_subplot(111)
1021            for dat in data:
1022                ax.plot(x0_vals, dat, ls='-', marker='')
1023
1024            if logscale is True:
1025                ax.set_yscale('log')
1026
1027            ax.set_xlabel(r'$x_0 / a$')
1028            plt.title(name)
1029            plt.draw()
1030
1031    def dump(self, filename, datatype="json.gz", **kwargs):
1032        """Dumps the Corr into a file of chosen type
1033        Parameters
1034        ----------
1035        filename : str
1036            Name of the file to be saved.
1037        datatype : str
1038            Format of the exported file. Supported formats include
1039            "json.gz" and "pickle"
1040        path : str
1041            specifies a custom path for the file (default '.')
1042        """
1043        if datatype == "json.gz":
1044            from .input.json import dump_to_json
1045            if 'path' in kwargs:
1046                file_name = kwargs.get('path') + '/' + filename
1047            else:
1048                file_name = filename
1049            dump_to_json(self, file_name)
1050        elif datatype == "pickle":
1051            dump_object(self, filename, **kwargs)
1052        else:
1053            raise ValueError("Unknown datatype " + str(datatype))
1054
1055    def print(self, print_range=None):
1056        print(self.__repr__(print_range))
1057
1058    def __repr__(self, print_range=None):
1059        if print_range is None:
1060            print_range = [0, None]
1061
1062        content_string = ""
1063        content_string += "Corr T=" + str(self.T) + " N=" + str(self.N) + "\n"  # +" filled with"+ str(type(self.content[0][0])) there should be a good solution here
1064
1065        if self.tag is not None:
1066            content_string += "Description: " + self.tag + "\n"
1067        if self.N != 1:
1068            return content_string
1069
1070        if print_range[1]:
1071            print_range[1] += 1
1072        content_string += 'x0/a\tCorr(x0/a)\n------------------\n'
1073        for i, sub_corr in enumerate(self.content[print_range[0]:print_range[1]]):
1074            if sub_corr is None:
1075                content_string += str(i + print_range[0]) + '\n'
1076            else:
1077                content_string += str(i + print_range[0])
1078                for element in sub_corr:
1079                    content_string += f"\t{element:+2}"
1080                content_string += '\n'
1081        return content_string
1082
1083    def __str__(self):
1084        return self.__repr__()
1085
1086    # We define the basic operations, that can be performed with correlators.
1087    # While */+- get defined here, they only work for Corr*Obs and not Obs*Corr.
1088    # This is because Obs*Corr checks Obs.__mul__ first and does not catch an exception.
1089    # One could try and tell Obs to check if the y in __mul__ is a Corr and
1090
1091    __array_priority__ = 10000
1092
1093    def __eq__(self, y):
1094        if isinstance(y, Corr):
1095            comp = np.asarray(y.content, dtype=object)
1096        else:
1097            comp = np.asarray(y)
1098        return np.asarray(self.content, dtype=object) == comp
1099
1100    __hash__ = None
1101
1102    def __add__(self, y):
1103        if isinstance(y, Corr):
1104            if ((self.N != y.N) or (self.T != y.T)):
1105                raise ValueError("Addition of Corrs with different shape")
1106            newcontent = []
1107            for t in range(self.T):
1108                if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]):
1109                    newcontent.append(None)
1110                else:
1111                    newcontent.append(self.content[t] + y.content[t])
1112            return Corr(newcontent)
1113
1114        elif isinstance(y, (Obs, int, float, CObs, complex)):
1115            newcontent = []
1116            for t in range(self.T):
1117                if _check_for_none(self, self.content[t]):
1118                    newcontent.append(None)
1119                else:
1120                    newcontent.append(self.content[t] + y)
1121            return Corr(newcontent, prange=self.prange)
1122        elif isinstance(y, np.ndarray):
1123            if y.shape == (self.T,):
1124                return Corr(list((np.array(self.content).T + y).T))
1125            else:
1126                raise ValueError("operands could not be broadcast together")
1127        else:
1128            raise TypeError("Corr + wrong type")
1129
1130    def __mul__(self, y):
1131        if isinstance(y, Corr):
1132            if not ((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T):
1133                raise ValueError("Multiplication of Corr object requires N=N or N=1 and T=T")
1134            newcontent = []
1135            for t in range(self.T):
1136                if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]):
1137                    newcontent.append(None)
1138                else:
1139                    newcontent.append(self.content[t] * y.content[t])
1140            return Corr(newcontent)
1141
1142        elif isinstance(y, (Obs, int, float, CObs, complex)):
1143            newcontent = []
1144            for t in range(self.T):
1145                if _check_for_none(self, self.content[t]):
1146                    newcontent.append(None)
1147                else:
1148                    newcontent.append(self.content[t] * y)
1149            return Corr(newcontent, prange=self.prange)
1150        elif isinstance(y, np.ndarray):
1151            if y.shape == (self.T,):
1152                return Corr(list((np.array(self.content).T * y).T))
1153            else:
1154                raise ValueError("operands could not be broadcast together")
1155        else:
1156            raise TypeError("Corr * wrong type")
1157
1158    def __matmul__(self, y):
1159        if isinstance(y, np.ndarray):
1160            if y.ndim != 2 or y.shape[0] != y.shape[1]:
1161                raise ValueError("Can only multiply correlators by square matrices.")
1162            if not self.N == y.shape[0]:
1163                raise ValueError("matmul: mismatch of matrix dimensions")
1164            newcontent = []
1165            for t in range(self.T):
1166                if _check_for_none(self, self.content[t]):
1167                    newcontent.append(None)
1168                else:
1169                    newcontent.append(self.content[t] @ y)
1170            return Corr(newcontent)
1171        elif isinstance(y, Corr):
1172            if not self.N == y.N:
1173                raise ValueError("matmul: mismatch of matrix dimensions")
1174            newcontent = []
1175            for t in range(self.T):
1176                if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]):
1177                    newcontent.append(None)
1178                else:
1179                    newcontent.append(self.content[t] @ y.content[t])
1180            return Corr(newcontent)
1181
1182        else:
1183            return NotImplemented
1184
1185    def __rmatmul__(self, y):
1186        if isinstance(y, np.ndarray):
1187            if y.ndim != 2 or y.shape[0] != y.shape[1]:
1188                raise ValueError("Can only multiply correlators by square matrices.")
1189            if not self.N == y.shape[0]:
1190                raise ValueError("matmul: mismatch of matrix dimensions")
1191            newcontent = []
1192            for t in range(self.T):
1193                if _check_for_none(self, self.content[t]):
1194                    newcontent.append(None)
1195                else:
1196                    newcontent.append(y @ self.content[t])
1197            return Corr(newcontent)
1198        else:
1199            return NotImplemented
1200
1201    def __truediv__(self, y):
1202        if isinstance(y, Corr):
1203            if not ((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T):
1204                raise ValueError("Multiplication of Corr object requires N=N or N=1 and T=T")
1205            newcontent = []
1206            for t in range(self.T):
1207                if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]):
1208                    newcontent.append(None)
1209                else:
1210                    newcontent.append(self.content[t] / y.content[t])
1211            for t in range(self.T):
1212                if _check_for_none(self, newcontent[t]):
1213                    continue
1214                if np.isnan(np.sum(newcontent[t]).value):
1215                    newcontent[t] = None
1216
1217            if all([item is None for item in newcontent]):
1218                raise ValueError("Division returns completely undefined correlator")
1219            return Corr(newcontent)
1220
1221        elif isinstance(y, (Obs, CObs)):
1222            if isinstance(y, Obs):
1223                if y.value == 0:
1224                    raise ValueError('Division by zero will return undefined correlator')
1225            if isinstance(y, CObs):
1226                if y.is_zero():
1227                    raise ValueError('Division by zero will return undefined correlator')
1228
1229            newcontent = []
1230            for t in range(self.T):
1231                if _check_for_none(self, self.content[t]):
1232                    newcontent.append(None)
1233                else:
1234                    newcontent.append(self.content[t] / y)
1235            return Corr(newcontent, prange=self.prange)
1236
1237        elif isinstance(y, (int, float)):
1238            if y == 0:
1239                raise ValueError('Division by zero will return undefined correlator')
1240            newcontent = []
1241            for t in range(self.T):
1242                if _check_for_none(self, self.content[t]):
1243                    newcontent.append(None)
1244                else:
1245                    newcontent.append(self.content[t] / y)
1246            return Corr(newcontent, prange=self.prange)
1247        elif isinstance(y, np.ndarray):
1248            if y.shape == (self.T,):
1249                return Corr(list((np.array(self.content).T / y).T))
1250            else:
1251                raise ValueError("operands could not be broadcast together")
1252        else:
1253            raise TypeError('Corr / wrong type')
1254
1255    def __neg__(self):
1256        newcontent = [None if _check_for_none(self, item) else -1. * item for item in self.content]
1257        return Corr(newcontent, prange=self.prange)
1258
1259    def __sub__(self, y):
1260        return self + (-y)
1261
1262    def __pow__(self, y):
1263        if isinstance(y, (Obs, int, float, CObs)):
1264            newcontent = [None if _check_for_none(self, item) else item**y for item in self.content]
1265            return Corr(newcontent, prange=self.prange)
1266        else:
1267            raise TypeError('Type of exponent not supported')
1268
1269    def __abs__(self):
1270        newcontent = [None if _check_for_none(self, item) else np.abs(item) for item in self.content]
1271        return Corr(newcontent, prange=self.prange)
1272
1273    # The numpy functions:
1274    def sqrt(self):
1275        return self ** 0.5
1276
1277    def log(self):
1278        newcontent = [None if _check_for_none(self, item) else np.log(item) for item in self.content]
1279        return Corr(newcontent, prange=self.prange)
1280
1281    def exp(self):
1282        newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content]
1283        return Corr(newcontent, prange=self.prange)
1284
1285    def _apply_func_to_corr(self, func):
1286        newcontent = [None if _check_for_none(self, item) else func(item) for item in self.content]
1287        for t in range(self.T):
1288            if _check_for_none(self, newcontent[t]):
1289                continue
1290            tmp_sum = np.sum(newcontent[t])
1291            if hasattr(tmp_sum, "value"):
1292                if np.isnan(tmp_sum.value):
1293                    newcontent[t] = None
1294        if all([item is None for item in newcontent]):
1295            raise ValueError('Operation returns undefined correlator')
1296        return Corr(newcontent)
1297
1298    def sin(self):
1299        return self._apply_func_to_corr(np.sin)
1300
1301    def cos(self):
1302        return self._apply_func_to_corr(np.cos)
1303
1304    def tan(self):
1305        return self._apply_func_to_corr(np.tan)
1306
1307    def sinh(self):
1308        return self._apply_func_to_corr(np.sinh)
1309
1310    def cosh(self):
1311        return self._apply_func_to_corr(np.cosh)
1312
1313    def tanh(self):
1314        return self._apply_func_to_corr(np.tanh)
1315
1316    def arcsin(self):
1317        return self._apply_func_to_corr(np.arcsin)
1318
1319    def arccos(self):
1320        return self._apply_func_to_corr(np.arccos)
1321
1322    def arctan(self):
1323        return self._apply_func_to_corr(np.arctan)
1324
1325    def arcsinh(self):
1326        return self._apply_func_to_corr(np.arcsinh)
1327
1328    def arccosh(self):
1329        return self._apply_func_to_corr(np.arccosh)
1330
1331    def arctanh(self):
1332        return self._apply_func_to_corr(np.arctanh)
1333
1334    # Right hand side operations (require tweak in main module to work)
1335    def __radd__(self, y):
1336        return self + y
1337
1338    def __rsub__(self, y):
1339        return -self + y
1340
1341    def __rmul__(self, y):
1342        return self * y
1343
1344    def __rtruediv__(self, y):
1345        return (self / y) ** (-1)
1346
1347    @property
1348    def real(self):
1349        def return_real(obs_OR_cobs):
1350            if isinstance(obs_OR_cobs.flatten()[0], CObs):
1351                return np.vectorize(lambda x: x.real)(obs_OR_cobs)
1352            else:
1353                return obs_OR_cobs
1354
1355        return self._apply_func_to_corr(return_real)
1356
1357    @property
1358    def imag(self):
1359        def return_imag(obs_OR_cobs):
1360            if isinstance(obs_OR_cobs.flatten()[0], CObs):
1361                return np.vectorize(lambda x: x.imag)(obs_OR_cobs)
1362            else:
1363                return obs_OR_cobs * 0  # So it stays the right type
1364
1365        return self._apply_func_to_corr(return_imag)
1366
1367    def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
1368        r''' Project large correlation matrix to lowest states
1369
1370        This method can be used to reduce the size of an (N x N) correlation matrix
1371        to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise
1372        is still small.
1373
1374        Parameters
1375        ----------
1376        Ntrunc: int
1377            Rank of the target matrix.
1378        tproj: int
1379            Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method.
1380            The default value is 3.
1381        t0proj: int
1382            Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly
1383            discouraged for O(a) improved theories, since the correctness of the procedure
1384            cannot be granted in this case. The default value is 2.
1385        basematrix : Corr
1386            Correlation matrix that is used to determine the eigenvectors of the
1387            lowest states based on a GEVP. basematrix is taken to be the Corr itself if
1388            is is not specified.
1389
1390        Notes
1391        -----
1392        We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving
1393        the GEVP $$C(t) v_n(t, t_0) = \lambda_n(t, t_0) C(t_0) v_n(t, t_0)$$ where $t \equiv t_\mathrm{proj}$
1394        and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the
1395        resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via
1396        $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large
1397        correlation matrix and to remove some noise that is added by irrelevant operators.
1398        This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated
1399        bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$.
1400        '''
1401
1402        if self.N == 1:
1403            raise ValueError('Method cannot be applied to one-dimensional correlators.')
1404        if basematrix is None:
1405            basematrix = self
1406        if Ntrunc >= basematrix.N:
1407            raise ValueError(f'Cannot truncate using Ntrunc >= {basematrix.N}')
1408        if basematrix.N != self.N:
1409            raise ValueError('basematrix and targetmatrix have to be of the same size.')
1410
1411        evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc]
1412
1413        tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object)
1414        rmat = []
1415        for t in range(basematrix.T):
1416            if self.content[t] is None:
1417                rmat.append(None)
1418            else:
1419                for i in range(Ntrunc):
1420                    for j in range(Ntrunc):
1421                        tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j]
1422                rmat.append(np.copy(tmpmat))
1423
1424        return Corr(rmat)
1425
1426
1427def _sort_vectors(vec_set_in, ts):
1428    """Helper function used to find a set of Eigenvectors consistent over all timeslices"""
1429
1430    if isinstance(vec_set_in[ts][0][0], Obs):
1431        vec_set = [anp.vectorize(float)(vi) if vi is not None else vi for vi in vec_set_in]
1432    else:
1433        vec_set = vec_set_in
1434    reference_sorting = np.array(vec_set[ts])
1435    N = reference_sorting.shape[0]
1436    sorted_vec_set = []
1437    for t in range(len(vec_set)):
1438        if vec_set[t] is None:
1439            sorted_vec_set.append(None)
1440        elif not t == ts:
1441            perms = [list(o) for o in permutations([i for i in range(N)], N)]
1442            best_score = 0
1443            for perm in perms:
1444                current_score = 1
1445                for k in range(N):
1446                    new_sorting = reference_sorting.copy()
1447                    new_sorting[perm[k], :] = vec_set[t][k]
1448                    current_score *= abs(np.linalg.det(new_sorting))
1449                if current_score > best_score:
1450                    best_score = current_score
1451                    best_perm = perm
1452            sorted_vec_set.append([vec_set_in[t][k] for k in best_perm])
1453        else:
1454            sorted_vec_set.append(vec_set_in[t])
1455
1456    return sorted_vec_set
1457
1458
1459def _check_for_none(corr, entry):
1460    """Checks if entry for correlator corr is None"""
1461    return len(list(filter(None, np.asarray(entry).flatten()))) < corr.N ** 2
1462
1463
1464def _GEVP_solver(Gt, G0, method='eigh', chol_inv=None):
1465    r"""Helper function for solving the GEVP and sorting the eigenvectors.
1466
1467    Solves $G(t)v_i=\lambda_i G(t_0)v_i$ and returns the eigenvectors v_i
1468
1469    The helper function assumes that both provided matrices are symmetric and
1470    only processes the lower triangular part of both matrices. In case the matrices
1471    are not symmetric the upper triangular parts are effectively discarded.
1472
1473    Parameters
1474    ----------
1475    Gt : array
1476        The correlator at time t for the left hand side of the GEVP
1477    G0 : array
1478        The correlator at time t0 for the right hand side of the GEVP
1479    Method used to solve the GEVP.
1480       - "eigh": Use scipy.linalg.eigh to solve the GEVP.
1481       - "cholesky": Use manually implemented solution via the Cholesky decomposition.
1482    chol_inv : array, optional
1483        Inverse of the Cholesky decomposition of G0. May be provided to
1484        speed up the computation in the case of method=='cholesky'
1485
1486    """
1487    if isinstance(G0[0][0], Obs):
1488        vector_obs = True
1489    else:
1490        vector_obs = False
1491
1492    if method == 'cholesky':
1493        if vector_obs:
1494            cholesky = linalg.cholesky
1495            inv = linalg.inv
1496            eigv = linalg.eigv
1497            matmul = linalg.matmul
1498        else:
1499            cholesky = np.linalg.cholesky
1500            inv = np.linalg.inv
1501
1502            def eigv(x, **kwargs):
1503                return np.linalg.eigh(x)[1]
1504
1505            def matmul(*operands):
1506                return np.linalg.multi_dot(operands)
1507        N = Gt.shape[0]
1508        output = [[] for j in range(N)]
1509        if chol_inv is None:
1510            chol = cholesky(G0)  # This will automatically report if the matrix is not pos-def
1511            chol_inv = inv(chol)
1512
1513        try:
1514            new_matrix = matmul(chol_inv, Gt, chol_inv.T)
1515            ev = eigv(new_matrix)
1516            ev = matmul(chol_inv.T, ev)
1517            output = np.flip(ev, axis=1).T
1518        except (np.linalg.LinAlgError, TypeError, ValueError):  # The above code can fail because of linalg-errors or because the entry of the corr is None
1519            for s in range(N):
1520                output[s] = None
1521        return output
1522    elif method == 'eigh':
1523        return scipy.linalg.eigh(Gt, G0, lower=True)[1].T[::-1]