pyerrors.correlators

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