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) or self.content[t][0].value / self.content[t + 1][0].value < 0:
 731                    newcontent.append(None)
 732                else:
 733                    newcontent.append(self.content[t] / self.content[t + 1])
 734            if (all([x is None for x in newcontent])):
 735                raise ValueError('m_eff is undefined at all timeslices')
 736
 737            return np.log(Corr(newcontent, padding=[0, 1]))
 738
 739        elif variant == 'logsym':
 740            newcontent = []
 741            for t in range(1, self.T - 1):
 742                if ((self.content[t - 1] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0) or self.content[t - 1][0].value / self.content[t + 1][0].value < 0:
 743                    newcontent.append(None)
 744                else:
 745                    newcontent.append(self.content[t - 1] / self.content[t + 1])
 746            if (all([x is None for x in newcontent])):
 747                raise ValueError('m_eff is undefined at all timeslices')
 748
 749            return np.log(Corr(newcontent, padding=[1, 1])) / 2
 750
 751        elif variant in ['periodic', 'cosh', 'sinh']:
 752            if variant in ['periodic', 'cosh']:
 753                func = anp.cosh
 754            else:
 755                func = anp.sinh
 756
 757            def root_function(x, d):
 758                return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d
 759
 760            newcontent = []
 761            for t in range(self.T - 1):
 762                if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t + 1][0].value == 0):
 763                    newcontent.append(None)
 764                # Fill the two timeslices in the middle of the lattice with their predecessors
 765                elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]:
 766                    newcontent.append(newcontent[-1])
 767                elif self.content[t][0].value / self.content[t + 1][0].value < 0:
 768                    newcontent.append(None)
 769                else:
 770                    newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess)))
 771            if (all([x is None for x in newcontent])):
 772                raise ValueError('m_eff is undefined at all timeslices')
 773
 774            return Corr(newcontent, padding=[0, 1])
 775
 776        elif variant == 'arccosh':
 777            newcontent = []
 778            for t in range(1, self.T - 1):
 779                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):
 780                    newcontent.append(None)
 781                else:
 782                    newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t]))
 783            if (all([x is None for x in newcontent])):
 784                raise ValueError("m_eff is undefined at all timeslices")
 785            return np.arccosh(Corr(newcontent, padding=[1, 1]))
 786
 787        else:
 788            raise ValueError('Unknown variant.')
 789
 790    def fit(self, function, fitrange=None, silent=False, **kwargs):
 791        r'''Fits function to the data
 792
 793        Parameters
 794        ----------
 795        function : obj
 796            function to fit to the data. See fits.least_squares for details.
 797        fitrange : list
 798            Two element list containing the timeslices on which the fit is supposed to start and stop.
 799            Caution: This range is inclusive as opposed to standard python indexing.
 800            `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6.
 801            If not specified, self.prange or all timeslices are used.
 802        silent : bool
 803            Decides whether output is printed to the standard output.
 804        '''
 805        if self.N != 1:
 806            raise ValueError("Correlator must be projected before fitting")
 807
 808        if fitrange is None:
 809            if self.prange:
 810                fitrange = self.prange
 811            else:
 812                fitrange = [0, self.T - 1]
 813        else:
 814            if not isinstance(fitrange, list):
 815                raise TypeError("fitrange has to be a list with two elements")
 816            if len(fitrange) != 2:
 817                raise ValueError("fitrange has to have exactly two elements [fit_start, fit_stop]")
 818
 819        xs = np.array([x for x in range(fitrange[0], fitrange[1] + 1) if self.content[x] is not None])
 820        ys = np.array([self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if self.content[x] is not None])
 821        result = least_squares(xs, ys, function, silent=silent, **kwargs)
 822        return result
 823
 824    def plateau(self, plateau_range=None, method="fit", auto_gamma=False):
 825        """ Extract a plateau value from a Corr object
 826
 827        Parameters
 828        ----------
 829        plateau_range : list
 830            list with two entries, indicating the first and the last timeslice
 831            of the plateau region.
 832        method : str
 833            method to extract the plateau.
 834                'fit' fits a constant to the plateau region
 835                'avg', 'average' or 'mean' just average over the given timeslices.
 836        auto_gamma : bool
 837            apply gamma_method with default parameters to the Corr. Defaults to None
 838        """
 839        if not plateau_range:
 840            if self.prange:
 841                plateau_range = self.prange
 842            else:
 843                raise ValueError("no plateau range provided")
 844        if self.N != 1:
 845            raise ValueError("Correlator must be projected before getting a plateau.")
 846        if (all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])):
 847            raise ValueError("plateau is undefined at all timeslices in plateaurange.")
 848        if auto_gamma:
 849            self.gamma_method()
 850        if method == "fit":
 851            def const_func(a, t):
 852                return a[0]
 853            return self.fit(const_func, plateau_range)[0]
 854        elif method in ["avg", "average", "mean"]:
 855            returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None])
 856            return returnvalue
 857
 858        else:
 859            raise ValueError("Unsupported plateau method: " + method)
 860
 861    def set_prange(self, prange):
 862        """Sets the attribute prange of the Corr object."""
 863        if not len(prange) == 2:
 864            raise ValueError("prange must be a list or array with two values")
 865        if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))):
 866            raise TypeError("Start and end point must be integers")
 867        if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] <= prange[1]):
 868            raise ValueError("Start and end point must define a range in the interval 0,T")
 869
 870        self.prange = prange
 871
 872    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):
 873        """Plots the correlator using the tag of the correlator as label if available.
 874
 875        Parameters
 876        ----------
 877        x_range : list
 878            list of two values, determining the range of the x-axis e.g. [4, 8].
 879        comp : Corr or list of Corr
 880            Correlator or list of correlators which are plotted for comparison.
 881            The tags of these correlators are used as labels if available.
 882        logscale : bool
 883            Sets y-axis to logscale.
 884        plateau : Obs
 885            Plateau value to be visualized in the figure.
 886        fit_res : Fit_result
 887            Fit_result object to be visualized.
 888        fit_key : str
 889            Key for the fit function in Fit_result.fit_function (for combined fits).
 890        ylabel : str
 891            Label for the y-axis.
 892        save : str
 893            path to file in which the figure should be saved.
 894        auto_gamma : bool
 895            Apply the gamma method with standard parameters to all correlators and plateau values before plotting.
 896        hide_sigma : float
 897            Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors.
 898        references : list
 899            List of floating point values that are displayed as horizontal lines for reference.
 900        title : string
 901            Optional title of the figure.
 902        """
 903        if self.N != 1:
 904            raise ValueError("Correlator must be projected before plotting")
 905
 906        if auto_gamma:
 907            self.gamma_method()
 908
 909        if x_range is None:
 910            x_range = [0, self.T - 1]
 911
 912        fig = plt.figure()
 913        ax1 = fig.add_subplot(111)
 914
 915        x, y, y_err = self.plottable()
 916        if hide_sigma:
 917            hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
 918        else:
 919            hide_from = None
 920        ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag)
 921        if logscale:
 922            ax1.set_yscale('log')
 923        else:
 924            if y_range is None:
 925                try:
 926                    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)])
 927                    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)])
 928                    ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)])
 929                except Exception:
 930                    pass
 931            else:
 932                ax1.set_ylim(y_range)
 933        if comp:
 934            if isinstance(comp, (Corr, list)):
 935                for corr in comp if isinstance(comp, list) else [comp]:
 936                    if auto_gamma:
 937                        corr.gamma_method()
 938                    x, y, y_err = corr.plottable()
 939                    if hide_sigma:
 940                        hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
 941                    else:
 942                        hide_from = None
 943                    ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor'])
 944            else:
 945                raise TypeError("'comp' must be a correlator or a list of correlators.")
 946
 947        if plateau:
 948            if isinstance(plateau, Obs):
 949                if auto_gamma:
 950                    plateau.gamma_method()
 951                ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau))
 952                ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-')
 953            else:
 954                raise TypeError("'plateau' must be an Obs")
 955
 956        if references:
 957            if isinstance(references, list):
 958                for ref in references:
 959                    ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--')
 960            else:
 961                raise TypeError("'references' must be a list of floating pint values.")
 962
 963        if self.prange:
 964            ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',', color="black", zorder=0)
 965            ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',', color="black", zorder=0)
 966
 967        if fit_res:
 968            x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05)
 969            if isinstance(fit_res.fit_function, dict):
 970                if fit_key:
 971                    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)
 972                else:
 973                    raise ValueError("Please provide a 'fit_key' for visualizing combined fits.")
 974            else:
 975                ax1.plot(x_samples, fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples), ls='-', marker=',', lw=2)
 976
 977        ax1.set_xlabel(r'$x_0 / a$')
 978        if ylabel:
 979            ax1.set_ylabel(ylabel)
 980        ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5])
 981
 982        _handles, labels = ax1.get_legend_handles_labels()
 983        if labels:
 984            ax1.legend()
 985
 986        if title:
 987            plt.title(title)
 988
 989        plt.draw()
 990
 991        if save:
 992            if isinstance(save, str):
 993                fig.savefig(save, bbox_inches='tight')
 994            else:
 995                raise TypeError("'save' has to be a string.")
 996
 997    def spaghetti_plot(self, logscale=True):
 998        """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations.
 999
1000        Parameters
1001        ----------
1002        logscale : bool
1003            Determines whether the scale of the y-axis is logarithmic or standard.
1004        """
1005        if self.N != 1:
1006            raise ValueError("Correlator needs to be projected first.")
1007
1008        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]))
1009        x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content, strict=True) if o is not None]
1010
1011        for name in mc_names:
1012            data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T
1013
1014            fig = plt.figure()
1015            ax = fig.add_subplot(111)
1016            for dat in data:
1017                ax.plot(x0_vals, dat, ls='-', marker='')
1018
1019            if logscale is True:
1020                ax.set_yscale('log')
1021
1022            ax.set_xlabel(r'$x_0 / a$')
1023            plt.title(name)
1024            plt.draw()
1025
1026    def dump(self, filename, datatype="json.gz", **kwargs):
1027        """Dumps the Corr into a file of chosen type
1028        Parameters
1029        ----------
1030        filename : str
1031            Name of the file to be saved.
1032        datatype : str
1033            Format of the exported file. Supported formats include
1034            "json.gz" and "pickle"
1035        path : str
1036            specifies a custom path for the file (default '.')
1037        """
1038        if datatype == "json.gz":
1039            from .input.json import dump_to_json
1040            if 'path' in kwargs:
1041                file_name = kwargs.get('path') + '/' + filename
1042            else:
1043                file_name = filename
1044            dump_to_json(self, file_name)
1045        elif datatype == "pickle":
1046            dump_object(self, filename, **kwargs)
1047        else:
1048            raise ValueError("Unknown datatype " + str(datatype))
1049
1050    def print(self, print_range=None):
1051        print(self.__repr__(print_range))
1052
1053    def __repr__(self, print_range=None):
1054        if print_range is None:
1055            print_range = [0, None]
1056
1057        content_string = ""
1058        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
1059
1060        if self.tag is not None:
1061            content_string += "Description: " + self.tag + "\n"
1062        if self.N != 1:
1063            return content_string
1064
1065        if print_range[1]:
1066            print_range[1] += 1
1067        content_string += 'x0/a\tCorr(x0/a)\n------------------\n'
1068        for i, sub_corr in enumerate(self.content[print_range[0]:print_range[1]]):
1069            if sub_corr is None:
1070                content_string += str(i + print_range[0]) + '\n'
1071            else:
1072                content_string += str(i + print_range[0])
1073                for element in sub_corr:
1074                    content_string += f"\t{element:+2}"
1075                content_string += '\n'
1076        return content_string
1077
1078    def __str__(self):
1079        return self.__repr__()
1080
1081    # We define the basic operations, that can be performed with correlators.
1082    # While */+- get defined here, they only work for Corr*Obs and not Obs*Corr.
1083    # This is because Obs*Corr checks Obs.__mul__ first and does not catch an exception.
1084    # One could try and tell Obs to check if the y in __mul__ is a Corr and
1085
1086    __array_priority__ = 10000
1087
1088    def __eq__(self, y):
1089        if isinstance(y, Corr):
1090            comp = np.asarray(y.content, dtype=object)
1091        else:
1092            comp = np.asarray(y)
1093        return np.asarray(self.content, dtype=object) == comp
1094
1095    __hash__ = None
1096
1097    def __add__(self, y):
1098        if isinstance(y, Corr):
1099            if ((self.N != y.N) or (self.T != y.T)):
1100                raise ValueError("Addition of Corrs with different shape")
1101            newcontent = []
1102            for t in range(self.T):
1103                if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]):
1104                    newcontent.append(None)
1105                else:
1106                    newcontent.append(self.content[t] + y.content[t])
1107            return Corr(newcontent)
1108
1109        elif isinstance(y, (Obs, int, float, CObs, complex)):
1110            newcontent = []
1111            for t in range(self.T):
1112                if _check_for_none(self, self.content[t]):
1113                    newcontent.append(None)
1114                else:
1115                    newcontent.append(self.content[t] + y)
1116            return Corr(newcontent, prange=self.prange)
1117        elif isinstance(y, np.ndarray):
1118            if y.shape == (self.T,):
1119                return Corr(list((np.array(self.content).T + y).T))
1120            else:
1121                raise ValueError("operands could not be broadcast together")
1122        else:
1123            raise TypeError("Corr + wrong type")
1124
1125    def __mul__(self, y):
1126        if isinstance(y, Corr):
1127            if not ((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T):
1128                raise ValueError("Multiplication of Corr object requires N=N or N=1 and T=T")
1129            newcontent = []
1130            for t in range(self.T):
1131                if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]):
1132                    newcontent.append(None)
1133                else:
1134                    newcontent.append(self.content[t] * y.content[t])
1135            return Corr(newcontent)
1136
1137        elif isinstance(y, (Obs, int, float, CObs, complex)):
1138            newcontent = []
1139            for t in range(self.T):
1140                if _check_for_none(self, self.content[t]):
1141                    newcontent.append(None)
1142                else:
1143                    newcontent.append(self.content[t] * y)
1144            return Corr(newcontent, prange=self.prange)
1145        elif isinstance(y, np.ndarray):
1146            if y.shape == (self.T,):
1147                return Corr(list((np.array(self.content).T * y).T))
1148            else:
1149                raise ValueError("operands could not be broadcast together")
1150        else:
1151            raise TypeError("Corr * wrong type")
1152
1153    def __matmul__(self, y):
1154        if isinstance(y, np.ndarray):
1155            if y.ndim != 2 or y.shape[0] != y.shape[1]:
1156                raise ValueError("Can only multiply correlators by square matrices.")
1157            if not self.N == y.shape[0]:
1158                raise ValueError("matmul: mismatch of matrix dimensions")
1159            newcontent = []
1160            for t in range(self.T):
1161                if _check_for_none(self, self.content[t]):
1162                    newcontent.append(None)
1163                else:
1164                    newcontent.append(self.content[t] @ y)
1165            return Corr(newcontent)
1166        elif isinstance(y, Corr):
1167            if not self.N == y.N:
1168                raise ValueError("matmul: mismatch of matrix dimensions")
1169            newcontent = []
1170            for t in range(self.T):
1171                if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]):
1172                    newcontent.append(None)
1173                else:
1174                    newcontent.append(self.content[t] @ y.content[t])
1175            return Corr(newcontent)
1176
1177        else:
1178            return NotImplemented
1179
1180    def __rmatmul__(self, y):
1181        if isinstance(y, np.ndarray):
1182            if y.ndim != 2 or y.shape[0] != y.shape[1]:
1183                raise ValueError("Can only multiply correlators by square matrices.")
1184            if not self.N == y.shape[0]:
1185                raise ValueError("matmul: mismatch of matrix dimensions")
1186            newcontent = []
1187            for t in range(self.T):
1188                if _check_for_none(self, self.content[t]):
1189                    newcontent.append(None)
1190                else:
1191                    newcontent.append(y @ self.content[t])
1192            return Corr(newcontent)
1193        else:
1194            return NotImplemented
1195
1196    def __truediv__(self, y):
1197        if isinstance(y, Corr):
1198            if not ((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T):
1199                raise ValueError("Multiplication of Corr object requires N=N or N=1 and T=T")
1200            newcontent = []
1201            for t in range(self.T):
1202                if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]):
1203                    newcontent.append(None)
1204                else:
1205                    newcontent.append(self.content[t] / y.content[t])
1206            for t in range(self.T):
1207                if _check_for_none(self, newcontent[t]):
1208                    continue
1209                if np.isnan(np.sum(newcontent[t]).value):
1210                    newcontent[t] = None
1211
1212            if all([item is None for item in newcontent]):
1213                raise ValueError("Division returns completely undefined correlator")
1214            return Corr(newcontent)
1215
1216        elif isinstance(y, (Obs, CObs)):
1217            if isinstance(y, Obs):
1218                if y.value == 0:
1219                    raise ValueError('Division by zero will return undefined correlator')
1220            if isinstance(y, CObs):
1221                if y.is_zero():
1222                    raise ValueError('Division by zero will return undefined correlator')
1223
1224            newcontent = []
1225            for t in range(self.T):
1226                if _check_for_none(self, self.content[t]):
1227                    newcontent.append(None)
1228                else:
1229                    newcontent.append(self.content[t] / y)
1230            return Corr(newcontent, prange=self.prange)
1231
1232        elif isinstance(y, (int, float)):
1233            if y == 0:
1234                raise ValueError('Division by zero will return undefined correlator')
1235            newcontent = []
1236            for t in range(self.T):
1237                if _check_for_none(self, self.content[t]):
1238                    newcontent.append(None)
1239                else:
1240                    newcontent.append(self.content[t] / y)
1241            return Corr(newcontent, prange=self.prange)
1242        elif isinstance(y, np.ndarray):
1243            if y.shape == (self.T,):
1244                return Corr(list((np.array(self.content).T / y).T))
1245            else:
1246                raise ValueError("operands could not be broadcast together")
1247        else:
1248            raise TypeError('Corr / wrong type')
1249
1250    def __neg__(self):
1251        newcontent = [None if _check_for_none(self, item) else -1. * item for item in self.content]
1252        return Corr(newcontent, prange=self.prange)
1253
1254    def __sub__(self, y):
1255        return self + (-y)
1256
1257    def __pow__(self, y):
1258        if isinstance(y, (Obs, int, float, CObs)):
1259            newcontent = [None if _check_for_none(self, item) else item**y for item in self.content]
1260            return Corr(newcontent, prange=self.prange)
1261        else:
1262            raise TypeError('Type of exponent not supported')
1263
1264    def __abs__(self):
1265        newcontent = [None if _check_for_none(self, item) else np.abs(item) for item in self.content]
1266        return Corr(newcontent, prange=self.prange)
1267
1268    # The numpy functions:
1269    def sqrt(self):
1270        return self ** 0.5
1271
1272    def log(self):
1273        newcontent = [None if _check_for_none(self, item) else np.log(item) for item in self.content]
1274        return Corr(newcontent, prange=self.prange)
1275
1276    def exp(self):
1277        newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content]
1278        return Corr(newcontent, prange=self.prange)
1279
1280    def _apply_func_to_corr(self, func):
1281        newcontent = [None if _check_for_none(self, item) else func(item) for item in self.content]
1282        for t in range(self.T):
1283            if _check_for_none(self, newcontent[t]):
1284                continue
1285            tmp_sum = np.sum(newcontent[t])
1286            if hasattr(tmp_sum, "value"):
1287                if np.isnan(tmp_sum.value):
1288                    newcontent[t] = None
1289        if all([item is None for item in newcontent]):
1290            raise ValueError('Operation returns undefined correlator')
1291        return Corr(newcontent)
1292
1293    def sin(self):
1294        return self._apply_func_to_corr(np.sin)
1295
1296    def cos(self):
1297        return self._apply_func_to_corr(np.cos)
1298
1299    def tan(self):
1300        return self._apply_func_to_corr(np.tan)
1301
1302    def sinh(self):
1303        return self._apply_func_to_corr(np.sinh)
1304
1305    def cosh(self):
1306        return self._apply_func_to_corr(np.cosh)
1307
1308    def tanh(self):
1309        return self._apply_func_to_corr(np.tanh)
1310
1311    def arcsin(self):
1312        return self._apply_func_to_corr(np.arcsin)
1313
1314    def arccos(self):
1315        return self._apply_func_to_corr(np.arccos)
1316
1317    def arctan(self):
1318        return self._apply_func_to_corr(np.arctan)
1319
1320    def arcsinh(self):
1321        return self._apply_func_to_corr(np.arcsinh)
1322
1323    def arccosh(self):
1324        return self._apply_func_to_corr(np.arccosh)
1325
1326    def arctanh(self):
1327        return self._apply_func_to_corr(np.arctanh)
1328
1329    # Right hand side operations (require tweak in main module to work)
1330    def __radd__(self, y):
1331        return self + y
1332
1333    def __rsub__(self, y):
1334        return -self + y
1335
1336    def __rmul__(self, y):
1337        return self * y
1338
1339    def __rtruediv__(self, y):
1340        return (self / y) ** (-1)
1341
1342    @property
1343    def real(self):
1344        def return_real(obs_OR_cobs):
1345            if isinstance(obs_OR_cobs.flatten()[0], CObs):
1346                return np.vectorize(lambda x: x.real)(obs_OR_cobs)
1347            else:
1348                return obs_OR_cobs
1349
1350        return self._apply_func_to_corr(return_real)
1351
1352    @property
1353    def imag(self):
1354        def return_imag(obs_OR_cobs):
1355            if isinstance(obs_OR_cobs.flatten()[0], CObs):
1356                return np.vectorize(lambda x: x.imag)(obs_OR_cobs)
1357            else:
1358                return obs_OR_cobs * 0  # So it stays the right type
1359
1360        return self._apply_func_to_corr(return_imag)
1361
1362    def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
1363        r''' Project large correlation matrix to lowest states
1364
1365        This method can be used to reduce the size of an (N x N) correlation matrix
1366        to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise
1367        is still small.
1368
1369        Parameters
1370        ----------
1371        Ntrunc: int
1372            Rank of the target matrix.
1373        tproj: int
1374            Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method.
1375            The default value is 3.
1376        t0proj: int
1377            Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly
1378            discouraged for O(a) improved theories, since the correctness of the procedure
1379            cannot be granted in this case. The default value is 2.
1380        basematrix : Corr
1381            Correlation matrix that is used to determine the eigenvectors of the
1382            lowest states based on a GEVP. basematrix is taken to be the Corr itself if
1383            is is not specified.
1384
1385        Notes
1386        -----
1387        We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving
1388        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}$
1389        and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the
1390        resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via
1391        $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large
1392        correlation matrix and to remove some noise that is added by irrelevant operators.
1393        This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated
1394        bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$.
1395        '''
1396
1397        if self.N == 1:
1398            raise ValueError('Method cannot be applied to one-dimensional correlators.')
1399        if basematrix is None:
1400            basematrix = self
1401        if Ntrunc >= basematrix.N:
1402            raise ValueError(f'Cannot truncate using Ntrunc >= {basematrix.N}')
1403        if basematrix.N != self.N:
1404            raise ValueError('basematrix and targetmatrix have to be of the same size.')
1405
1406        evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc]
1407
1408        tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object)
1409        rmat = []
1410        for t in range(basematrix.T):
1411            if self.content[t] is None:
1412                rmat.append(None)
1413            else:
1414                for i in range(Ntrunc):
1415                    for j in range(Ntrunc):
1416                        tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j]
1417                rmat.append(np.copy(tmpmat))
1418
1419        return Corr(rmat)
1420
1421
1422def _sort_vectors(vec_set_in, ts):
1423    """Helper function used to find a set of Eigenvectors consistent over all timeslices"""
1424
1425    if isinstance(vec_set_in[ts][0][0], Obs):
1426        vec_set = [anp.vectorize(float)(vi) if vi is not None else vi for vi in vec_set_in]
1427    else:
1428        vec_set = vec_set_in
1429    reference_sorting = np.array(vec_set[ts])
1430    N = reference_sorting.shape[0]
1431    sorted_vec_set = []
1432    for t in range(len(vec_set)):
1433        if vec_set[t] is None:
1434            sorted_vec_set.append(None)
1435        elif not t == ts:
1436            perms = [list(o) for o in permutations([i for i in range(N)], N)]
1437            best_score = 0
1438            for perm in perms:
1439                current_score = 1
1440                for k in range(N):
1441                    new_sorting = reference_sorting.copy()
1442                    new_sorting[perm[k], :] = vec_set[t][k]
1443                    current_score *= abs(np.linalg.det(new_sorting))
1444                if current_score > best_score:
1445                    best_score = current_score
1446                    best_perm = perm
1447            sorted_vec_set.append([vec_set_in[t][k] for k in best_perm])
1448        else:
1449            sorted_vec_set.append(vec_set_in[t])
1450
1451    return sorted_vec_set
1452
1453
1454def _check_for_none(corr, entry):
1455    """Checks if entry for correlator corr is None"""
1456    return len(list(filter(None, np.asarray(entry).flatten()))) < corr.N ** 2
1457
1458
1459def _GEVP_solver(Gt, G0, method='eigh', chol_inv=None):
1460    r"""Helper function for solving the GEVP and sorting the eigenvectors.
1461
1462    Solves $G(t)v_i=\lambda_i G(t_0)v_i$ and returns the eigenvectors v_i
1463
1464    The helper function assumes that both provided matrices are symmetric and
1465    only processes the lower triangular part of both matrices. In case the matrices
1466    are not symmetric the upper triangular parts are effectively discarded.
1467
1468    Parameters
1469    ----------
1470    Gt : array
1471        The correlator at time t for the left hand side of the GEVP
1472    G0 : array
1473        The correlator at time t0 for the right hand side of the GEVP
1474    Method used to solve the GEVP.
1475       - "eigh": Use scipy.linalg.eigh to solve the GEVP.
1476       - "cholesky": Use manually implemented solution via the Cholesky decomposition.
1477    chol_inv : array, optional
1478        Inverse of the Cholesky decomposition of G0. May be provided to
1479        speed up the computation in the case of method=='cholesky'
1480
1481    """
1482    if isinstance(G0[0][0], Obs):
1483        vector_obs = True
1484    else:
1485        vector_obs = False
1486
1487    if method == 'cholesky':
1488        if vector_obs:
1489            cholesky = linalg.cholesky
1490            inv = linalg.inv
1491            eigv = linalg.eigv
1492            matmul = linalg.matmul
1493        else:
1494            cholesky = np.linalg.cholesky
1495            inv = np.linalg.inv
1496
1497            def eigv(x, **kwargs):
1498                return np.linalg.eigh(x)[1]
1499
1500            def matmul(*operands):
1501                return np.linalg.multi_dot(operands)
1502        N = Gt.shape[0]
1503        output = [[] for j in range(N)]
1504        if chol_inv is None:
1505            chol = cholesky(G0)  # This will automatically report if the matrix is not pos-def
1506            chol_inv = inv(chol)
1507
1508        try:
1509            new_matrix = matmul(chol_inv, Gt, chol_inv.T)
1510            ev = eigv(new_matrix)
1511            ev = matmul(chol_inv.T, ev)
1512            output = np.flip(ev, axis=1).T
1513        except (np.linalg.LinAlgError, TypeError, ValueError):  # The above code can fail because of linalg-errors or because the entry of the corr is None
1514            for s in range(N):
1515                output[s] = None
1516        return output
1517    elif method == 'eigh':
1518        return scipy.linalg.eigh(Gt, G0, lower=True)[1].T[::-1]