# pyerrors.correlators

View Source

       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]

class Corr: View Source

      15class Corr:
      16    r"""The class for a correlator (time dependent sequence of pe.Obs).
      17
      18    Everything, this class does, can be achieved using lists or arrays of Obs.
      19    But it is simply more convenient to have a dedicated object for correlators.
      20    One often wants to add or multiply correlators of the same length at every timeslice and it is inconvenient
      21    to iterate over all timeslices for every operation. This is especially true, when dealing with matrices.
      22
      23    The correlator can have two types of content: An Obs at every timeslice OR a matrix at every timeslice.
      24    Other dependency (eg. spatial) are not supported.
      25
      26    The Corr class can also deal with missing measurements or paddings for fixed boundary conditions.
      27    The missing entries are represented via the `None` object.
      28
      29    Initialization
      30    --------------
      31    A simple correlator can be initialized with a list or a one-dimensional array of `Obs` or `Cobs`
      32    ```python
      33    corr11 = pe.Corr([obs1, obs2])
      34    corr11 = pe.Corr(np.array([obs1, obs2]))
      35    ```
      36    A matrix-valued correlator can either be initialized via a two-dimensional array of `Corr` objects
      37    ```python
      38    matrix_corr = pe.Corr(np.array([[corr11, corr12], [corr21, corr22]]))
      39    ```
      40    or alternatively via a three-dimensional array of `Obs` or `CObs` of shape (T, N, N) where T is
      41    the temporal extent of the correlator and N is the dimension of the matrix.
      42    """
      43
      44    __slots__ = ["content", "N", "T", "tag", "prange"]
      45
      46    def __init__(self, data_input, padding=[0, 0], prange=None):
      47        """ Initialize a Corr object.
      48
      49        Parameters
      50        ----------
      51        data_input : list or array
      52            list of Obs or list of arrays of Obs or array of Corrs (see class docstring for details).
      53        padding : list, optional
      54            List with two entries where the first labels the padding
      55            at the front of the correlator and the second the padding
      56            at the back.
      57        prange : list, optional
      58            List containing the first and last timeslice of the plateau
      59            region identified for this correlator.
      60        """
      61
      62        if isinstance(data_input, np.ndarray):
      63            if data_input.ndim == 1:
      64                data_input = list(data_input)
      65            elif data_input.ndim == 2:
      66                if not data_input.shape[0] == data_input.shape[1]:
      67                    raise ValueError("Array needs to be square.")
      68                if not all([isinstance(item, Corr) for item in data_input.flatten()]):
      69                    raise ValueError("If the input is an array, its elements must be of type pe.Corr.")
      70                if not all([item.N == 1 for item in data_input.flatten()]):
      71                    raise ValueError("Can only construct matrix correlator from single valued correlators.")
      72                if not len(set([item.T for item in data_input.flatten()])) == 1:
      73                    raise ValueError("All input Correlators must be defined over the same timeslices.")
      74
      75                T = data_input[0, 0].T
      76                N = data_input.shape[0]
      77                input_as_list = []
      78                for t in range(T):
      79                    if any([(item.content[t] is None) for item in data_input.flatten()]):
      80                        if not all([(item.content[t] is None) for item in data_input.flatten()]):
      81                            warnings.warn("Input ill-defined at different timeslices. Conversion leads to data loss.!", RuntimeWarning)
      82                        input_as_list.append(None)
      83                    else:
      84                        array_at_timeslace = np.empty([N, N], dtype="object")
      85                        for i in range(N):
      86                            for j in range(N):
      87                                array_at_timeslace[i, j] = data_input[i, j][t]
      88                        input_as_list.append(array_at_timeslace)
      89                data_input = input_as_list
      90            elif data_input.ndim == 3:
      91                if not data_input.shape[1] == data_input.shape[2]:
      92                    raise ValueError("Array needs to be square.")
      93                data_input = list(data_input)
      94            else:
      95                raise ValueError("Arrays with ndim>3 not supported.")
      96
      97        if isinstance(data_input, list):
      98
      99            if all([isinstance(item, (Obs, CObs)) or item is None for item in data_input]):
     100                _assert_equal_properties([o for o in data_input if o is not None])
     101                self.content = [np.asarray([item]) if item is not None else None for item in data_input]
     102                self.N = 1
     103            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]):
     104                self.content = data_input
     105                noNull = [a for a in self.content if not (a is None)]  # To check if the matrices are correct for all undefined elements
     106                self.N = noNull[0].shape[0]
     107                if self.N > 1 and noNull[0].shape[0] != noNull[0].shape[1]:
     108                    raise ValueError("Smearing matrices are not NxN.")
     109                if (not all([item.shape == noNull[0].shape for item in noNull])):
     110                    raise ValueError("Items in data_input are not of identical shape." + str(noNull))
     111            else:
     112                raise TypeError("'data_input' contains item of wrong type.")
     113        else:
     114            raise TypeError("Data input was not given as list or correct array.")
     115
     116        self.tag = None
     117
     118        # An undefined timeslice is represented by the None object
     119        self.content = [None] * padding[0] + self.content + [None] * padding[1]
     120        self.T = len(self.content)
     121        self.prange = prange
     122
     123    def __getitem__(self, idx):
     124        """Return the content of timeslice idx"""
     125        if self.content[idx] is None:
     126            return None
     127        elif len(self.content[idx]) == 1:
     128            return self.content[idx][0]
     129        else:
     130            return self.content[idx]
     131
     132    @property
     133    def reweighted(self):
     134        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]])
     135        if np.all(bool_array == 1):
     136            return True
     137        elif np.all(bool_array == 0):
     138            return False
     139        else:
     140            raise Exception("Reweighting status of correlator corrupted.")
     141
     142    def gamma_method(self, **kwargs):
     143        """Apply the gamma method to the content of the Corr."""
     144        for item in self.content:
     145            if not (item is None):
     146                if self.N == 1:
     147                    item[0].gamma_method(**kwargs)
     148                else:
     149                    for i in range(self.N):
     150                        for j in range(self.N):
     151                            item[i, j].gamma_method(**kwargs)
     152
     153    gm = gamma_method
     154
     155    def projected(self, vector_l=None, vector_r=None, normalize=False):
     156        """We need to project the Correlator with a Vector to get a single value at each timeslice.
     157
     158        The method can use one or two vectors.
     159        If two are specified it returns v1@G@v2 (the order might be very important.)
     160        By default it will return the lowest source, which usually means unsmeared-unsmeared (0,0), but it does not have to
     161        """
     162        if self.N == 1:
     163            raise Exception("Trying to project a Corr, that already has N=1.")
     164
     165        if vector_l is None:
     166            vector_l, vector_r = np.asarray([1.] + (self.N - 1) * [0.]), np.asarray([1.] + (self.N - 1) * [0.])
     167        elif (vector_r is None):
     168            vector_r = vector_l
     169        if isinstance(vector_l, list) and not isinstance(vector_r, list):
     170            if len(vector_l) != self.T:
     171                raise Exception("Length of vector list must be equal to T")
     172            vector_r = [vector_r] * self.T
     173        if isinstance(vector_r, list) and not isinstance(vector_l, list):
     174            if len(vector_r) != self.T:
     175                raise Exception("Length of vector list must be equal to T")
     176            vector_l = [vector_l] * self.T
     177
     178        if not isinstance(vector_l, list):
     179            if not vector_l.shape == vector_r.shape == (self.N,):
     180                raise Exception("Vectors are of wrong shape!")
     181            if normalize:
     182                vector_l, vector_r = vector_l / np.sqrt((vector_l @ vector_l)), vector_r / np.sqrt(vector_r @ vector_r)
     183            newcontent = [None if _check_for_none(self, item) else np.asarray([vector_l.T @ item @ vector_r]) for item in self.content]
     184
     185        else:
     186            # There are no checks here yet. There are so many possible scenarios, where this can go wrong.
     187            if normalize:
     188                for t in range(self.T):
     189                    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])
     190
     191            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)]
     192        return Corr(newcontent)
     193
     194    def item(self, i, j):
     195        """Picks the element [i,j] from every matrix and returns a correlator containing one Obs per timeslice.
     196
     197        Parameters
     198        ----------
     199        i : int
     200            First index to be picked.
     201        j : int
     202            Second index to be picked.
     203        """
     204        if self.N == 1:
     205            raise Exception("Trying to pick item from projected Corr")
     206        newcontent = [None if (item is None) else item[i, j] for item in self.content]
     207        return Corr(newcontent)
     208
     209    def plottable(self):
     210        """Outputs the correlator in a plotable format.
     211
     212        Outputs three lists containing the timeslice index, the value on each
     213        timeslice and the error on each timeslice.
     214        """
     215        if self.N != 1:
     216            raise Exception("Can only make Corr[N=1] plottable")
     217        x_list = [x for x in range(self.T) if not self.content[x] is None]
     218        y_list = [y[0].value for y in self.content if y is not None]
     219        y_err_list = [y[0].dvalue for y in self.content if y is not None]
     220
     221        return x_list, y_list, y_err_list
     222
     223    def symmetric(self):
     224        """ Symmetrize the correlator around x0=0."""
     225        if self.N != 1:
     226            raise Exception('symmetric cannot be safely applied to multi-dimensional correlators.')
     227        if self.T % 2 != 0:
     228            raise Exception("Can not symmetrize odd T")
     229
     230        if self.content[0] is not None:
     231            if np.argmax(np.abs([o[0].value if o is not None else 0 for o in self.content])) != 0:
     232                warnings.warn("Correlator does not seem to be symmetric around x0=0.", RuntimeWarning)
     233
     234        newcontent = [self.content[0]]
     235        for t in range(1, self.T):
     236            if (self.content[t] is None) or (self.content[self.T - t] is None):
     237                newcontent.append(None)
     238            else:
     239                newcontent.append(0.5 * (self.content[t] + self.content[self.T - t]))
     240        if (all([x is None for x in newcontent])):
     241            raise Exception("Corr could not be symmetrized: No redundant values")
     242        return Corr(newcontent, prange=self.prange)
     243
     244    def anti_symmetric(self):
     245        """Anti-symmetrize the correlator around x0=0."""
     246        if self.N != 1:
     247            raise TypeError('anti_symmetric cannot be safely applied to multi-dimensional correlators.')
     248        if self.T % 2 != 0:
     249            raise Exception("Can not symmetrize odd T")
     250
     251        test = 1 * self
     252        test.gamma_method()
     253        if not all([o.is_zero_within_error(3) for o in test.content[0]]):
     254            warnings.warn("Correlator does not seem to be anti-symmetric around x0=0.", RuntimeWarning)
     255
     256        newcontent = [self.content[0]]
     257        for t in range(1, self.T):
     258            if (self.content[t] is None) or (self.content[self.T - t] is None):
     259                newcontent.append(None)
     260            else:
     261                newcontent.append(0.5 * (self.content[t] - self.content[self.T - t]))
     262        if (all([x is None for x in newcontent])):
     263            raise Exception("Corr could not be symmetrized: No redundant values")
     264        return Corr(newcontent, prange=self.prange)
     265
     266    def is_matrix_symmetric(self):
     267        """Checks whether a correlator matrices is symmetric on every timeslice."""
     268        if self.N == 1:
     269            raise TypeError("Only works for correlator matrices.")
     270        for t in range(self.T):
     271            if self[t] is None:
     272                continue
     273            for i in range(self.N):
     274                for j in range(i + 1, self.N):
     275                    if self[t][i, j] is self[t][j, i]:
     276                        continue
     277                    if hash(self[t][i, j]) != hash(self[t][j, i]):
     278                        return False
     279        return True
     280
     281    def trace(self):
     282        """Calculates the per-timeslice trace of a correlator matrix."""
     283        if self.N == 1:
     284            raise ValueError("Only works for correlator matrices.")
     285        newcontent = []
     286        for t in range(self.T):
     287            if _check_for_none(self, self.content[t]):
     288                newcontent.append(None)
     289            else:
     290                newcontent.append(np.trace(self.content[t]))
     291        return Corr(newcontent)
     292
     293    def matrix_symmetric(self):
     294        """Symmetrizes the correlator matrices on every timeslice."""
     295        if self.N == 1:
     296            raise Exception("Trying to symmetrize a correlator matrix, that already has N=1.")
     297        if self.is_matrix_symmetric():
     298            return 1.0 * self
     299        else:
     300            transposed = [None if _check_for_none(self, G) else G.T for G in self.content]
     301            return 0.5 * (Corr(transposed) + self)
     302
     303    def GEVP(self, t0, ts=None, sort="Eigenvalue", vector_obs=False, **kwargs):
     304        r'''Solve the generalized eigenvalue problem on the correlator matrix and returns the corresponding eigenvectors.
     305
     306        The eigenvectors are sorted according to the descending eigenvalues, the zeroth eigenvector(s) correspond to the
     307        largest eigenvalue(s). The eigenvector(s) for the individual states can be accessed via slicing
     308        ```python
     309        C.GEVP(t0=2)[0]  # Ground state vector(s)
     310        C.GEVP(t0=2)[:3]  # Vectors for the lowest three states
     311        ```
     312
     313        Parameters
     314        ----------
     315        t0 : int
     316            The time t0 for the right hand side of the GEVP according to $G(t)v_i=\lambda_i G(t_0)v_i$
     317        ts : int
     318            fixed time $G(t_s)v_i=\lambda_i G(t_0)v_i$ if sort=None.
     319            If sort="Eigenvector" it gives a reference point for the sorting method.
     320        sort : string
     321            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.
     322            - "Eigenvalue": The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice. (default)
     323            - "Eigenvector": Use the method described in arXiv:2004.10472 to find the set of v(t) belonging to the state.
     324              The reference state is identified by its eigenvalue at $t=t_s$.
     325            - None: The GEVP is solved only at ts, no sorting is necessary
     326        vector_obs : bool
     327            If True, uncertainties are propagated in the eigenvector computation (default False).
     328
     329        Other Parameters
     330        ----------------
     331        state : int
     332           Returns only the vector(s) for a specified state. The lowest state is zero.
     333        method : str
     334           Method used to solve the GEVP.
     335           - "eigh": Use scipy.linalg.eigh to solve the GEVP. (default for vector_obs=False)
     336           - "cholesky": Use manually implemented solution via the Cholesky decomposition. Automatically chosen if vector_obs==True.
     337        '''
     338
     339        if self.N == 1:
     340            raise Exception("GEVP methods only works on correlator matrices and not single correlators.")
     341        if ts is not None:
     342            if (ts <= t0):
     343                raise Exception("ts has to be larger than t0.")
     344
     345        if "sorted_list" in kwargs:
     346            warnings.warn("Argument 'sorted_list' is deprecated, use 'sort' instead.", DeprecationWarning)
     347            sort = kwargs.get("sorted_list")
     348
     349        if self.is_matrix_symmetric():
     350            symmetric_corr = self
     351        else:
     352            symmetric_corr = self.matrix_symmetric()
     353
     354        def _get_mat_at_t(t, vector_obs=vector_obs):
     355            if vector_obs:
     356                return symmetric_corr[t]
     357            else:
     358                return np.vectorize(lambda x: x.value)(symmetric_corr[t])
     359        G0 = _get_mat_at_t(t0)
     360
     361        method = kwargs.get('method', 'eigh')
     362        if vector_obs:
     363            chol = linalg.cholesky(G0)
     364            chol_inv = linalg.inv(chol)
     365            method = 'cholesky'
     366        else:
     367            chol = np.linalg.cholesky(_get_mat_at_t(t0, vector_obs=False))  # Check if matrix G0 is positive-semidefinite.
     368            if method == 'cholesky':
     369                chol_inv = np.linalg.inv(chol)
     370            else:
     371                chol_inv = None
     372
     373        if sort is None:
     374            if (ts is None):
     375                raise Exception("ts is required if sort=None.")
     376            if (self.content[t0] is None) or (self.content[ts] is None):
     377                raise Exception("Corr not defined at t0/ts.")
     378            Gt = _get_mat_at_t(ts)
     379            reordered_vecs = _GEVP_solver(Gt, G0, method=method, chol_inv=chol_inv)
     380            if kwargs.get('auto_gamma', False) and vector_obs:
     381                [[o.gm() for o in ev if isinstance(o, Obs)] for ev in reordered_vecs]
     382
     383        elif sort in ["Eigenvalue", "Eigenvector"]:
     384            if sort == "Eigenvalue" and ts is not None:
     385                warnings.warn("ts has no effect when sorting by eigenvalue is chosen.", RuntimeWarning)
     386            all_vecs = [None] * (t0 + 1)
     387            for t in range(t0 + 1, self.T):
     388                try:
     389                    Gt = _get_mat_at_t(t)
     390                    all_vecs.append(_GEVP_solver(Gt, G0, method=method, chol_inv=chol_inv))
     391                except Exception:
     392                    all_vecs.append(None)
     393            if sort == "Eigenvector":
     394                if ts is None:
     395                    raise Exception("ts is required for the Eigenvector sorting method.")
     396                all_vecs = _sort_vectors(all_vecs, ts)
     397
     398            reordered_vecs = [[v[s] if v is not None else None for v in all_vecs] for s in range(self.N)]
     399            if kwargs.get('auto_gamma', False) and vector_obs:
     400                [[[o.gm() for o in evn] for evn in ev if evn is not None] for ev in reordered_vecs]
     401        else:
     402            raise Exception("Unknown value for 'sort'. Choose 'Eigenvalue', 'Eigenvector' or None.")
     403
     404        if "state" in kwargs:
     405            return reordered_vecs[kwargs.get("state")]
     406        else:
     407            return reordered_vecs
     408
     409    def Eigenvalue(self, t0, ts=None, state=0, sort="Eigenvalue", **kwargs):
     410        """Determines the eigenvalue of the GEVP by solving and projecting the correlator
     411
     412        Parameters
     413        ----------
     414        state : int
     415            The state one is interested in ordered by energy. The lowest state is zero.
     416
     417        All other parameters are identical to the ones of Corr.GEVP.
     418        """
     419        vec = self.GEVP(t0, ts=ts, sort=sort, **kwargs)[state]
     420        return self.projected(vec)
     421
     422    def Hankel(self, N, periodic=False):
     423        """Constructs an NxN Hankel matrix
     424
     425        C(t) c(t+1) ... c(t+n-1)
     426        C(t+1) c(t+2) ... c(t+n)
     427        .................
     428        C(t+(n-1)) c(t+n) ... c(t+2(n-1))
     429
     430        Parameters
     431        ----------
     432        N : int
     433            Dimension of the Hankel matrix
     434        periodic : bool, optional
     435            determines whether the matrix is extended periodically
     436        """
     437
     438        if self.N != 1:
     439            raise Exception("Multi-operator Prony not implemented!")
     440
     441        array = np.empty([N, N], dtype="object")
     442        new_content = []
     443        for t in range(self.T):
     444            new_content.append(array.copy())
     445
     446        def wrap(i):
     447            while i >= self.T:
     448                i -= self.T
     449            return i
     450
     451        for t in range(self.T):
     452            for i in range(N):
     453                for j in range(N):
     454                    if periodic:
     455                        new_content[t][i, j] = self.content[wrap(t + i + j)][0]
     456                    elif (t + i + j) >= self.T:
     457                        new_content[t] = None
     458                    else:
     459                        new_content[t][i, j] = self.content[t + i + j][0]
     460
     461        return Corr(new_content)
     462
     463    def roll(self, dt):
     464        """Periodically shift the correlator by dt timeslices
     465
     466        Parameters
     467        ----------
     468        dt : int
     469            number of timeslices
     470        """
     471        return Corr(list(np.roll(np.array(self.content, dtype=object), dt, axis=0)))
     472
     473    def reverse(self):
     474        """Reverse the time ordering of the Corr"""
     475        return Corr(self.content[:: -1])
     476
     477    def thin(self, spacing=2, offset=0):
     478        """Thin out a correlator to suppress correlations
     479
     480        Parameters
     481        ----------
     482        spacing : int
     483            Keep only every 'spacing'th entry of the correlator
     484        offset : int
     485            Offset the equal spacing
     486        """
     487        new_content = []
     488        for t in range(self.T):
     489            if (offset + t) % spacing != 0:
     490                new_content.append(None)
     491            else:
     492                new_content.append(self.content[t])
     493        return Corr(new_content)
     494
     495    def correlate(self, partner):
     496        """Correlate the correlator with another correlator or Obs
     497
     498        Parameters
     499        ----------
     500        partner : Obs or Corr
     501            partner to correlate the correlator with.
     502            Can either be an Obs which is correlated with all entries of the
     503            correlator or a Corr of same length.
     504        """
     505        if self.N != 1:
     506            raise Exception("Only one-dimensional correlators can be safely correlated.")
     507        new_content = []
     508        for x0, t_slice in enumerate(self.content):
     509            if _check_for_none(self, t_slice):
     510                new_content.append(None)
     511            else:
     512                if isinstance(partner, Corr):
     513                    if _check_for_none(partner, partner.content[x0]):
     514                        new_content.append(None)
     515                    else:
     516                        new_content.append(np.array([correlate(o, partner.content[x0][0]) for o in t_slice]))
     517                elif isinstance(partner, Obs):  # Should this include CObs?
     518                    new_content.append(np.array([correlate(o, partner) for o in t_slice]))
     519                else:
     520                    raise Exception("Can only correlate with an Obs or a Corr.")
     521
     522        return Corr(new_content)
     523
     524    def reweight(self, weight, **kwargs):
     525        """Reweight the correlator.
     526
     527        Parameters
     528        ----------
     529        weight : Obs
     530            Reweighting factor. An Observable that has to be defined on a superset of the
     531            configurations in obs[i].idl for all i.
     532        all_configs : bool
     533            if True, the reweighted observables are normalized by the average of
     534            the reweighting factor on all configurations in weight.idl and not
     535            on the configurations in obs[i].idl.
     536        """
     537        if self.N != 1:
     538            raise Exception("Reweighting only implemented for one-dimensional correlators.")
     539        new_content = []
     540        for t_slice in self.content:
     541            if _check_for_none(self, t_slice):
     542                new_content.append(None)
     543            else:
     544                new_content.append(np.array(reweight(weight, t_slice, **kwargs)))
     545        return Corr(new_content)
     546
     547    def T_symmetry(self, partner, parity=+1):
     548        """Return the time symmetry average of the correlator and its partner
     549
     550        Parameters
     551        ----------
     552        partner : Corr
     553            Time symmetry partner of the Corr
     554        parity : int
     555            Parity quantum number of the correlator, can be +1 or -1
     556        """
     557        if self.N != 1:
     558            raise Exception("T_symmetry only implemented for one-dimensional correlators.")
     559        if not isinstance(partner, Corr):
     560            raise Exception("T partner has to be a Corr object.")
     561        if parity not in [+1, -1]:
     562            raise Exception("Parity has to be +1 or -1.")
     563        T_partner = parity * partner.reverse()
     564
     565        t_slices = []
     566        test = (self - T_partner)
     567        test.gamma_method()
     568        for x0, t_slice in enumerate(test.content):
     569            if t_slice is not None:
     570                if not t_slice[0].is_zero_within_error(5):
     571                    t_slices.append(x0)
     572        if t_slices:
     573            warnings.warn("T symmetry partners do not agree within 5 sigma on time slices " + str(t_slices) + ".", RuntimeWarning)
     574
     575        return (self + T_partner) / 2
     576
     577    def deriv(self, variant="symmetric"):
     578        """Return the first derivative of the correlator with respect to x0.
     579
     580        Parameters
     581        ----------
     582        variant : str
     583            decides which definition of the finite differences derivative is used.
     584            Available choice: symmetric, forward, backward, improved, log, default: symmetric
     585        """
     586        if self.N != 1:
     587            raise Exception("deriv only implemented for one-dimensional correlators.")
     588        if variant == "symmetric":
     589            newcontent = []
     590            for t in range(1, self.T - 1):
     591                if (self.content[t - 1] is None) or (self.content[t + 1] is None):
     592                    newcontent.append(None)
     593                else:
     594                    newcontent.append(0.5 * (self.content[t + 1] - self.content[t - 1]))
     595            if (all([x is None for x in newcontent])):
     596                raise Exception('Derivative is undefined at all timeslices')
     597            return Corr(newcontent, padding=[1, 1])
     598        elif variant == "forward":
     599            newcontent = []
     600            for t in range(self.T - 1):
     601                if (self.content[t] is None) or (self.content[t + 1] is None):
     602                    newcontent.append(None)
     603                else:
     604                    newcontent.append(self.content[t + 1] - self.content[t])
     605            if (all([x is None for x in newcontent])):
     606                raise Exception("Derivative is undefined at all timeslices")
     607            return Corr(newcontent, padding=[0, 1])
     608        elif variant == "backward":
     609            newcontent = []
     610            for t in range(1, self.T):
     611                if (self.content[t - 1] is None) or (self.content[t] is None):
     612                    newcontent.append(None)
     613                else:
     614                    newcontent.append(self.content[t] - self.content[t - 1])
     615            if (all([x is None for x in newcontent])):
     616                raise Exception("Derivative is undefined at all timeslices")
     617            return Corr(newcontent, padding=[1, 0])
     618        elif variant == "improved":
     619            newcontent = []
     620            for t in range(2, self.T - 2):
     621                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):
     622                    newcontent.append(None)
     623                else:
     624                    newcontent.append((1 / 12) * (self.content[t - 2] - 8 * self.content[t - 1] + 8 * self.content[t + 1] - self.content[t + 2]))
     625            if (all([x is None for x in newcontent])):
     626                raise Exception('Derivative is undefined at all timeslices')
     627            return Corr(newcontent, padding=[2, 2])
     628        elif variant == 'log':
     629            newcontent = []
     630            for t in range(self.T):
     631                if (self.content[t] is None) or (self.content[t] <= 0):
     632                    newcontent.append(None)
     633                else:
     634                    newcontent.append(np.log(self.content[t]))
     635            if (all([x is None for x in newcontent])):
     636                raise Exception("Log is undefined at all timeslices")
     637            logcorr = Corr(newcontent)
     638            return self * logcorr.deriv('symmetric')
     639        else:
     640            raise Exception("Unknown variant.")
     641
     642    def second_deriv(self, variant="symmetric"):
     643        r"""Return the second derivative of the correlator with respect to x0.
     644
     645        Parameters
     646        ----------
     647        variant : str
     648            decides which definition of the finite differences derivative is used.
     649            Available choice:
     650                - symmetric (default)
     651                    $$\tilde{\partial}^2_0 f(x_0) = f(x_0+1)-2f(x_0)+f(x_0-1)$$
     652                - big_symmetric
     653                    $$\partial^2_0 f(x_0) = \frac{f(x_0+2)-2f(x_0)+f(x_0-2)}{4}$$
     654                - improved
     655                    $$\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}$$
     656                - log
     657                    $$f(x) = \tilde{\partial}^2_0 log(f(x_0))+(\tilde{\partial}_0 log(f(x_0)))^2$$
     658        """
     659        if self.N != 1:
     660            raise Exception("second_deriv only implemented for one-dimensional correlators.")
     661        if variant == "symmetric":
     662            newcontent = []
     663            for t in range(1, self.T - 1):
     664                if (self.content[t - 1] is None) or (self.content[t + 1] is None):
     665                    newcontent.append(None)
     666                else:
     667                    newcontent.append((self.content[t + 1] - 2 * self.content[t] + self.content[t - 1]))
     668            if (all([x is None for x in newcontent])):
     669                raise Exception("Derivative is undefined at all timeslices")
     670            return Corr(newcontent, padding=[1, 1])
     671        elif variant == "big_symmetric":
     672            newcontent = []
     673            for t in range(2, self.T - 2):
     674                if (self.content[t - 2] is None) or (self.content[t + 2] is None):
     675                    newcontent.append(None)
     676                else:
     677                    newcontent.append((self.content[t + 2] - 2 * self.content[t] + self.content[t - 2]) / 4)
     678            if (all([x is None for x in newcontent])):
     679                raise Exception("Derivative is undefined at all timeslices")
     680            return Corr(newcontent, padding=[2, 2])
     681        elif variant == "improved":
     682            newcontent = []
     683            for t in range(2, self.T - 2):
     684                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):
     685                    newcontent.append(None)
     686                else:
     687                    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]))
     688            if (all([x is None for x in newcontent])):
     689                raise Exception("Derivative is undefined at all timeslices")
     690            return Corr(newcontent, padding=[2, 2])
     691        elif variant == 'log':
     692            newcontent = []
     693            for t in range(self.T):
     694                if (self.content[t] is None) or (self.content[t] <= 0):
     695                    newcontent.append(None)
     696                else:
     697                    newcontent.append(np.log(self.content[t]))
     698            if (all([x is None for x in newcontent])):
     699                raise Exception("Log is undefined at all timeslices")
     700            logcorr = Corr(newcontent)
     701            return self * (logcorr.second_deriv('symmetric') + (logcorr.deriv('symmetric'))**2)
     702        else:
     703            raise Exception("Unknown variant.")
     704
     705    def m_eff(self, variant='log', guess=1.0):
     706        """Returns the effective mass of the correlator as correlator object
     707
     708        Parameters
     709        ----------
     710        variant : str
     711            log : uses the standard effective mass log(C(t) / C(t+1))
     712            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.
     713            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.
     714            See, e.g., arXiv:1205.5380
     715            arccosh : Uses the explicit form of the symmetrized correlator (not recommended)
     716            logsym: uses the symmetric effective mass log(C(t-1) / C(t+1))/2
     717        guess : float
     718            guess for the root finder, only relevant for the root variant
     719        """
     720        if self.N != 1:
     721            raise Exception('Correlator must be projected before getting m_eff')
     722        if variant == 'log':
     723            newcontent = []
     724            for t in range(self.T - 1):
     725                if ((self.content[t] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0):
     726                    newcontent.append(None)
     727                elif self.content[t][0].value / self.content[t + 1][0].value < 0:
     728                    newcontent.append(None)
     729                else:
     730                    newcontent.append(self.content[t] / self.content[t + 1])
     731            if (all([x is None for x in newcontent])):
     732                raise Exception('m_eff is undefined at all timeslices')
     733
     734            return np.log(Corr(newcontent, padding=[0, 1]))
     735
     736        elif variant == 'logsym':
     737            newcontent = []
     738            for t in range(1, self.T - 1):
     739                if ((self.content[t - 1] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0):
     740                    newcontent.append(None)
     741                elif self.content[t - 1][0].value / self.content[t + 1][0].value < 0:
     742                    newcontent.append(None)
     743                else:
     744                    newcontent.append(self.content[t - 1] / self.content[t + 1])
     745            if (all([x is None for x in newcontent])):
     746                raise Exception('m_eff is undefined at all timeslices')
     747
     748            return np.log(Corr(newcontent, padding=[1, 1])) / 2
     749
     750        elif variant in ['periodic', 'cosh', 'sinh']:
     751            if variant in ['periodic', 'cosh']:
     752                func = anp.cosh
     753            else:
     754                func = anp.sinh
     755
     756            def root_function(x, d):
     757                return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d
     758
     759            newcontent = []
     760            for t in range(self.T - 1):
     761                if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t + 1][0].value == 0):
     762                    newcontent.append(None)
     763                # Fill the two timeslices in the middle of the lattice with their predecessors
     764                elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]:
     765                    newcontent.append(newcontent[-1])
     766                elif self.content[t][0].value / self.content[t + 1][0].value < 0:
     767                    newcontent.append(None)
     768                else:
     769                    newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess)))
     770            if (all([x is None for x in newcontent])):
     771                raise Exception('m_eff is undefined at all timeslices')
     772
     773            return Corr(newcontent, padding=[0, 1])
     774
     775        elif variant == 'arccosh':
     776            newcontent = []
     777            for t in range(1, self.T - 1):
     778                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):
     779                    newcontent.append(None)
     780                else:
     781                    newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t]))
     782            if (all([x is None for x in newcontent])):
     783                raise Exception("m_eff is undefined at all timeslices")
     784            return np.arccosh(Corr(newcontent, padding=[1, 1]))
     785
     786        else:
     787            raise Exception('Unknown variant.')
     788
     789    def fit(self, function, fitrange=None, silent=False, **kwargs):
     790        r'''Fits function to the data
     791
     792        Parameters
     793        ----------
     794        function : obj
     795            function to fit to the data. See fits.least_squares for details.
     796        fitrange : list
     797            Two element list containing the timeslices on which the fit is supposed to start and stop.
     798            Caution: This range is inclusive as opposed to standard python indexing.
     799            `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6.
     800            If not specified, self.prange or all timeslices are used.
     801        silent : bool
     802            Decides whether output is printed to the standard output.
     803        '''
     804        if self.N != 1:
     805            raise Exception("Correlator must be projected before fitting")
     806
     807        if fitrange is None:
     808            if self.prange:
     809                fitrange = self.prange
     810            else:
     811                fitrange = [0, self.T - 1]
     812        else:
     813            if not isinstance(fitrange, list):
     814                raise Exception("fitrange has to be a list with two elements")
     815            if len(fitrange) != 2:
     816                raise Exception("fitrange has to have exactly two elements [fit_start, fit_stop]")
     817
     818        xs = np.array([x for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None])
     819        ys = np.array([self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None])
     820        result = least_squares(xs, ys, function, silent=silent, **kwargs)
     821        return result
     822
     823    def plateau(self, plateau_range=None, method="fit", auto_gamma=False):
     824        """ Extract a plateau value from a Corr object
     825
     826        Parameters
     827        ----------
     828        plateau_range : list
     829            list with two entries, indicating the first and the last timeslice
     830            of the plateau region.
     831        method : str
     832            method to extract the plateau.
     833                'fit' fits a constant to the plateau region
     834                'avg', 'average' or 'mean' just average over the given timeslices.
     835        auto_gamma : bool
     836            apply gamma_method with default parameters to the Corr. Defaults to None
     837        """
     838        if not plateau_range:
     839            if self.prange:
     840                plateau_range = self.prange
     841            else:
     842                raise Exception("no plateau range provided")
     843        if self.N != 1:
     844            raise Exception("Correlator must be projected before getting a plateau.")
     845        if (all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])):
     846            raise Exception("plateau is undefined at all timeslices in plateaurange.")
     847        if auto_gamma:
     848            self.gamma_method()
     849        if method == "fit":
     850            def const_func(a, t):
     851                return a[0]
     852            return self.fit(const_func, plateau_range)[0]
     853        elif method in ["avg", "average", "mean"]:
     854            returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None])
     855            return returnvalue
     856
     857        else:
     858            raise Exception("Unsupported plateau method: " + method)
     859
     860    def set_prange(self, prange):
     861        """Sets the attribute prange of the Corr object."""
     862        if not len(prange) == 2:
     863            raise Exception("prange must be a list or array with two values")
     864        if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))):
     865            raise Exception("Start and end point must be integers")
     866        if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] < prange[1]):
     867            raise Exception("Start and end point must define a range in the interval 0,T")
     868
     869        self.prange = prange
     870        return
     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 Exception("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 Exception("'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 Exception("'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 Exception("'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 Exception("'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 Exception("Correlator needs to be projected first.")
    1007
    1008        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]))
    1009        x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) 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 Exception("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    def __add__(self, y):
    1096        if isinstance(y, Corr):
    1097            if ((self.N != y.N) or (self.T != y.T)):
    1098                raise Exception("Addition of Corrs with different shape")
    1099            newcontent = []
    1100            for t in range(self.T):
    1101                if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]):
    1102                    newcontent.append(None)
    1103                else:
    1104                    newcontent.append(self.content[t] + y.content[t])
    1105            return Corr(newcontent)
    1106
    1107        elif isinstance(y, (Obs, int, float, CObs, complex)):
    1108            newcontent = []
    1109            for t in range(self.T):
    1110                if _check_for_none(self, self.content[t]):
    1111                    newcontent.append(None)
    1112                else:
    1113                    newcontent.append(self.content[t] + y)
    1114            return Corr(newcontent, prange=self.prange)
    1115        elif isinstance(y, np.ndarray):
    1116            if y.shape == (self.T,):
    1117                return Corr(list((np.array(self.content).T + y).T))
    1118            else:
    1119                raise ValueError("operands could not be broadcast together")
    1120        else:
    1121            raise TypeError("Corr + wrong type")
    1122
    1123    def __mul__(self, y):
    1124        if isinstance(y, Corr):
    1125            if not ((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T):
    1126                raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T")
    1127            newcontent = []
    1128            for t in range(self.T):
    1129                if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]):
    1130                    newcontent.append(None)
    1131                else:
    1132                    newcontent.append(self.content[t] * y.content[t])
    1133            return Corr(newcontent)
    1134
    1135        elif isinstance(y, (Obs, int, float, CObs, complex)):
    1136            newcontent = []
    1137            for t in range(self.T):
    1138                if _check_for_none(self, self.content[t]):
    1139                    newcontent.append(None)
    1140                else:
    1141                    newcontent.append(self.content[t] * y)
    1142            return Corr(newcontent, prange=self.prange)
    1143        elif isinstance(y, np.ndarray):
    1144            if y.shape == (self.T,):
    1145                return Corr(list((np.array(self.content).T * y).T))
    1146            else:
    1147                raise ValueError("operands could not be broadcast together")
    1148        else:
    1149            raise TypeError("Corr * wrong type")
    1150
    1151    def __matmul__(self, y):
    1152        if isinstance(y, np.ndarray):
    1153            if y.ndim != 2 or y.shape[0] != y.shape[1]:
    1154                raise ValueError("Can only multiply correlators by square matrices.")
    1155            if not self.N == y.shape[0]:
    1156                raise ValueError("matmul: mismatch of matrix dimensions")
    1157            newcontent = []
    1158            for t in range(self.T):
    1159                if _check_for_none(self, self.content[t]):
    1160                    newcontent.append(None)
    1161                else:
    1162                    newcontent.append(self.content[t] @ y)
    1163            return Corr(newcontent)
    1164        elif isinstance(y, Corr):
    1165            if not self.N == y.N:
    1166                raise ValueError("matmul: mismatch of matrix dimensions")
    1167            newcontent = []
    1168            for t in range(self.T):
    1169                if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]):
    1170                    newcontent.append(None)
    1171                else:
    1172                    newcontent.append(self.content[t] @ y.content[t])
    1173            return Corr(newcontent)
    1174
    1175        else:
    1176            return NotImplemented
    1177
    1178    def __rmatmul__(self, y):
    1179        if isinstance(y, np.ndarray):
    1180            if y.ndim != 2 or y.shape[0] != y.shape[1]:
    1181                raise ValueError("Can only multiply correlators by square matrices.")
    1182            if not self.N == y.shape[0]:
    1183                raise ValueError("matmul: mismatch of matrix dimensions")
    1184            newcontent = []
    1185            for t in range(self.T):
    1186                if _check_for_none(self, self.content[t]):
    1187                    newcontent.append(None)
    1188                else:
    1189                    newcontent.append(y @ self.content[t])
    1190            return Corr(newcontent)
    1191        else:
    1192            return NotImplemented
    1193
    1194    def __truediv__(self, y):
    1195        if isinstance(y, Corr):
    1196            if not ((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T):
    1197                raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T")
    1198            newcontent = []
    1199            for t in range(self.T):
    1200                if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]):
    1201                    newcontent.append(None)
    1202                else:
    1203                    newcontent.append(self.content[t] / y.content[t])
    1204            for t in range(self.T):
    1205                if _check_for_none(self, newcontent[t]):
    1206                    continue
    1207                if np.isnan(np.sum(newcontent[t]).value):
    1208                    newcontent[t] = None
    1209
    1210            if all([item is None for item in newcontent]):
    1211                raise Exception("Division returns completely undefined correlator")
    1212            return Corr(newcontent)
    1213
    1214        elif isinstance(y, (Obs, CObs)):
    1215            if isinstance(y, Obs):
    1216                if y.value == 0:
    1217                    raise Exception('Division by zero will return undefined correlator')
    1218            if isinstance(y, CObs):
    1219                if y.is_zero():
    1220                    raise Exception('Division by zero will return undefined correlator')
    1221
    1222            newcontent = []
    1223            for t in range(self.T):
    1224                if _check_for_none(self, self.content[t]):
    1225                    newcontent.append(None)
    1226                else:
    1227                    newcontent.append(self.content[t] / y)
    1228            return Corr(newcontent, prange=self.prange)
    1229
    1230        elif isinstance(y, (int, float)):
    1231            if y == 0:
    1232                raise Exception('Division by zero will return undefined correlator')
    1233            newcontent = []
    1234            for t in range(self.T):
    1235                if _check_for_none(self, self.content[t]):
    1236                    newcontent.append(None)
    1237                else:
    1238                    newcontent.append(self.content[t] / y)
    1239            return Corr(newcontent, prange=self.prange)
    1240        elif isinstance(y, np.ndarray):
    1241            if y.shape == (self.T,):
    1242                return Corr(list((np.array(self.content).T / y).T))
    1243            else:
    1244                raise ValueError("operands could not be broadcast together")
    1245        else:
    1246            raise TypeError('Corr / wrong type')
    1247
    1248    def __neg__(self):
    1249        newcontent = [None if _check_for_none(self, item) else -1. * item for item in self.content]
    1250        return Corr(newcontent, prange=self.prange)
    1251
    1252    def __sub__(self, y):
    1253        return self + (-y)
    1254
    1255    def __pow__(self, y):
    1256        if isinstance(y, (Obs, int, float, CObs)):
    1257            newcontent = [None if _check_for_none(self, item) else item**y for item in self.content]
    1258            return Corr(newcontent, prange=self.prange)
    1259        else:
    1260            raise TypeError('Type of exponent not supported')
    1261
    1262    def __abs__(self):
    1263        newcontent = [None if _check_for_none(self, item) else np.abs(item) for item in self.content]
    1264        return Corr(newcontent, prange=self.prange)
    1265
    1266    # The numpy functions:
    1267    def sqrt(self):
    1268        return self ** 0.5
    1269
    1270    def log(self):
    1271        newcontent = [None if _check_for_none(self, item) else np.log(item) for item in self.content]
    1272        return Corr(newcontent, prange=self.prange)
    1273
    1274    def exp(self):
    1275        newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content]
    1276        return Corr(newcontent, prange=self.prange)
    1277
    1278    def _apply_func_to_corr(self, func):
    1279        newcontent = [None if _check_for_none(self, item) else func(item) for item in self.content]
    1280        for t in range(self.T):
    1281            if _check_for_none(self, newcontent[t]):
    1282                continue
    1283            tmp_sum = np.sum(newcontent[t])
    1284            if hasattr(tmp_sum, "value"):
    1285                if np.isnan(tmp_sum.value):
    1286                    newcontent[t] = None
    1287        if all([item is None for item in newcontent]):
    1288            raise Exception('Operation returns undefined correlator')
    1289        return Corr(newcontent)
    1290
    1291    def sin(self):
    1292        return self._apply_func_to_corr(np.sin)
    1293
    1294    def cos(self):
    1295        return self._apply_func_to_corr(np.cos)
    1296
    1297    def tan(self):
    1298        return self._apply_func_to_corr(np.tan)
    1299
    1300    def sinh(self):
    1301        return self._apply_func_to_corr(np.sinh)
    1302
    1303    def cosh(self):
    1304        return self._apply_func_to_corr(np.cosh)
    1305
    1306    def tanh(self):
    1307        return self._apply_func_to_corr(np.tanh)
    1308
    1309    def arcsin(self):
    1310        return self._apply_func_to_corr(np.arcsin)
    1311
    1312    def arccos(self):
    1313        return self._apply_func_to_corr(np.arccos)
    1314
    1315    def arctan(self):
    1316        return self._apply_func_to_corr(np.arctan)
    1317
    1318    def arcsinh(self):
    1319        return self._apply_func_to_corr(np.arcsinh)
    1320
    1321    def arccosh(self):
    1322        return self._apply_func_to_corr(np.arccosh)
    1323
    1324    def arctanh(self):
    1325        return self._apply_func_to_corr(np.arctanh)
    1326
    1327    # Right hand side operations (require tweak in main module to work)
    1328    def __radd__(self, y):
    1329        return self + y
    1330
    1331    def __rsub__(self, y):
    1332        return -self + y
    1333
    1334    def __rmul__(self, y):
    1335        return self * y
    1336
    1337    def __rtruediv__(self, y):
    1338        return (self / y) ** (-1)
    1339
    1340    @property
    1341    def real(self):
    1342        def return_real(obs_OR_cobs):
    1343            if isinstance(obs_OR_cobs.flatten()[0], CObs):
    1344                return np.vectorize(lambda x: x.real)(obs_OR_cobs)
    1345            else:
    1346                return obs_OR_cobs
    1347
    1348        return self._apply_func_to_corr(return_real)
    1349
    1350    @property
    1351    def imag(self):
    1352        def return_imag(obs_OR_cobs):
    1353            if isinstance(obs_OR_cobs.flatten()[0], CObs):
    1354                return np.vectorize(lambda x: x.imag)(obs_OR_cobs)
    1355            else:
    1356                return obs_OR_cobs * 0  # So it stays the right type
    1357
    1358        return self._apply_func_to_corr(return_imag)
    1359
    1360    def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
    1361        r''' Project large correlation matrix to lowest states
    1362
    1363        This method can be used to reduce the size of an (N x N) correlation matrix
    1364        to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise
    1365        is still small.
    1366
    1367        Parameters
    1368        ----------
    1369        Ntrunc: int
    1370            Rank of the target matrix.
    1371        tproj: int
    1372            Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method.
    1373            The default value is 3.
    1374        t0proj: int
    1375            Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly
    1376            discouraged for O(a) improved theories, since the correctness of the procedure
    1377            cannot be granted in this case. The default value is 2.
    1378        basematrix : Corr
    1379            Correlation matrix that is used to determine the eigenvectors of the
    1380            lowest states based on a GEVP. basematrix is taken to be the Corr itself if
    1381            is is not specified.
    1382
    1383        Notes
    1384        -----
    1385        We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving
    1386        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}$
    1387        and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the
    1388        resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via
    1389        $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large
    1390        correlation matrix and to remove some noise that is added by irrelevant operators.
    1391        This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated
    1392        bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$.
    1393        '''
    1394
    1395        if self.N == 1:
    1396            raise Exception('Method cannot be applied to one-dimensional correlators.')
    1397        if basematrix is None:
    1398            basematrix = self
    1399        if Ntrunc >= basematrix.N:
    1400            raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N))
    1401        if basematrix.N != self.N:
    1402            raise Exception('basematrix and targetmatrix have to be of the same size.')
    1403
    1404        evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc]
    1405
    1406        tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object)
    1407        rmat = []
    1408        for t in range(basematrix.T):
    1409            for i in range(Ntrunc):
    1410                for j in range(Ntrunc):
    1411                    tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j]
    1412            rmat.append(np.copy(tmpmat))
    1413
    1414        newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)]
    1415        return Corr(newcontent)

The class for a correlator (time dependent sequence of pe.Obs).

Everything, this class does, can be achieved using lists or arrays of Obs. But it is simply more convenient to have a dedicated object for correlators. One often wants to add or multiply correlators of the same length at every timeslice and it is inconvenient to iterate over all timeslices for every operation. This is especially true, when dealing with matrices.

The correlator can have two types of content: An Obs at every timeslice OR a matrix at every timeslice. Other dependency (eg. spatial) are not supported.

The Corr class can also deal with missing measurements or paddings for fixed boundary conditions. The missing entries are represented via the `None` object.

###### Initialization

A simple correlator can be initialized with a list or a one-dimensional array of `Obs` or `Cobs`

    corr11 = pe.Corr([obs1, obs2])
    corr11 = pe.Corr(np.array([obs1, obs2]))

A matrix-valued correlator can either be initialized via a two-dimensional array of `Corr` objects

    matrix_corr = pe.Corr(np.array([[corr11, corr12], [corr21, corr22]]))

or alternatively via a three-dimensional array of `Obs` or `CObs` of shape (T, N, N) where T is the temporal extent of the correlator and N is the dimension of the matrix.

Corr(data_input, padding=[0, 0], prange=None) View Source

     46    def __init__(self, data_input, padding=[0, 0], prange=None):
     47        """ Initialize a Corr object.
     48
     49        Parameters
     50        ----------
     51        data_input : list or array
     52            list of Obs or list of arrays of Obs or array of Corrs (see class docstring for details).
     53        padding : list, optional
     54            List with two entries where the first labels the padding
     55            at the front of the correlator and the second the padding
     56            at the back.
     57        prange : list, optional
     58            List containing the first and last timeslice of the plateau
     59            region identified for this correlator.
     60        """
     61
     62        if isinstance(data_input, np.ndarray):
     63            if data_input.ndim == 1:
     64                data_input = list(data_input)
     65            elif data_input.ndim == 2:
     66                if not data_input.shape[0] == data_input.shape[1]:
     67                    raise ValueError("Array needs to be square.")
     68                if not all([isinstance(item, Corr) for item in data_input.flatten()]):
     69                    raise ValueError("If the input is an array, its elements must be of type pe.Corr.")
     70                if not all([item.N == 1 for item in data_input.flatten()]):
     71                    raise ValueError("Can only construct matrix correlator from single valued correlators.")
     72                if not len(set([item.T for item in data_input.flatten()])) == 1:
     73                    raise ValueError("All input Correlators must be defined over the same timeslices.")
     74
     75                T = data_input[0, 0].T
     76                N = data_input.shape[0]
     77                input_as_list = []
     78                for t in range(T):
     79                    if any([(item.content[t] is None) for item in data_input.flatten()]):
     80                        if not all([(item.content[t] is None) for item in data_input.flatten()]):
     81                            warnings.warn("Input ill-defined at different timeslices. Conversion leads to data loss.!", RuntimeWarning)
     82                        input_as_list.append(None)
     83                    else:
     84                        array_at_timeslace = np.empty([N, N], dtype="object")
     85                        for i in range(N):
     86                            for j in range(N):
     87                                array_at_timeslace[i, j] = data_input[i, j][t]
     88                        input_as_list.append(array_at_timeslace)
     89                data_input = input_as_list
     90            elif data_input.ndim == 3:
     91                if not data_input.shape[1] == data_input.shape[2]:
     92                    raise ValueError("Array needs to be square.")
     93                data_input = list(data_input)
     94            else:
     95                raise ValueError("Arrays with ndim>3 not supported.")
     96
     97        if isinstance(data_input, list):
     98
     99            if all([isinstance(item, (Obs, CObs)) or item is None for item in data_input]):
    100                _assert_equal_properties([o for o in data_input if o is not None])
    101                self.content = [np.asarray([item]) if item is not None else None for item in data_input]
    102                self.N = 1
    103            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]):
    104                self.content = data_input
    105                noNull = [a for a in self.content if not (a is None)]  # To check if the matrices are correct for all undefined elements
    106                self.N = noNull[0].shape[0]
    107                if self.N > 1 and noNull[0].shape[0] != noNull[0].shape[1]:
    108                    raise ValueError("Smearing matrices are not NxN.")
    109                if (not all([item.shape == noNull[0].shape for item in noNull])):
    110                    raise ValueError("Items in data_input are not of identical shape." + str(noNull))
    111            else:
    112                raise TypeError("'data_input' contains item of wrong type.")
    113        else:
    114            raise TypeError("Data input was not given as list or correct array.")
    115
    116        self.tag = None
    117
    118        # An undefined timeslice is represented by the None object
    119        self.content = [None] * padding[0] + self.content + [None] * padding[1]
    120        self.T = len(self.content)
    121        self.prange = prange

Initialize a Corr object.

###### Parameters

- data_input (list or array): list of Obs or list of arrays of Obs or array of Corrs (see class docstring for details).
- padding (list, optional): List with two entries where the first labels the padding at the front of the correlator and the second the padding at the back.
- prange (list, optional): List containing the first and last timeslice of the plateau region identified for this correlator.

tag
content
T
prange
reweighted View Source

    132    @property
    133    def reweighted(self):
    134        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]])
    135        if np.all(bool_array == 1):
    136            return True
    137        elif np.all(bool_array == 0):
    138            return False
    139        else:
    140            raise Exception("Reweighting status of correlator corrupted.")

def gamma_method(self, **kwargs): View Source

    142    def gamma_method(self, **kwargs):
    143        """Apply the gamma method to the content of the Corr."""
    144        for item in self.content:
    145            if not (item is None):
    146                if self.N == 1:
    147                    item[0].gamma_method(**kwargs)
    148                else:
    149                    for i in range(self.N):
    150                        for j in range(self.N):
    151                            item[i, j].gamma_method(**kwargs)

Apply the gamma method to the content of the Corr.

def gm(self, **kwargs): View Source

    142    def gamma_method(self, **kwargs):
    143        """Apply the gamma method to the content of the Corr."""
    144        for item in self.content:
    145            if not (item is None):
    146                if self.N == 1:
    147                    item[0].gamma_method(**kwargs)
    148                else:
    149                    for i in range(self.N):
    150                        for j in range(self.N):
    151                            item[i, j].gamma_method(**kwargs)

Apply the gamma method to the content of the Corr.

def projected(self, vector_l=None, vector_r=None, normalize=False): View Source

    155    def projected(self, vector_l=None, vector_r=None, normalize=False):
    156        """We need to project the Correlator with a Vector to get a single value at each timeslice.
    157
    158        The method can use one or two vectors.
    159        If two are specified it returns v1@G@v2 (the order might be very important.)
    160        By default it will return the lowest source, which usually means unsmeared-unsmeared (0,0), but it does not have to
    161        """
    162        if self.N == 1:
    163            raise Exception("Trying to project a Corr, that already has N=1.")
    164
    165        if vector_l is None:
    166            vector_l, vector_r = np.asarray([1.] + (self.N - 1) * [0.]), np.asarray([1.] + (self.N - 1) * [0.])
    167        elif (vector_r is None):
    168            vector_r = vector_l
    169        if isinstance(vector_l, list) and not isinstance(vector_r, list):
    170            if len(vector_l) != self.T:
    171                raise Exception("Length of vector list must be equal to T")
    172            vector_r = [vector_r] * self.T
    173        if isinstance(vector_r, list) and not isinstance(vector_l, list):
    174            if len(vector_r) != self.T:
    175                raise Exception("Length of vector list must be equal to T")
    176            vector_l = [vector_l] * self.T
    177
    178        if not isinstance(vector_l, list):
    179            if not vector_l.shape == vector_r.shape == (self.N,):
    180                raise Exception("Vectors are of wrong shape!")
    181            if normalize:
    182                vector_l, vector_r = vector_l / np.sqrt((vector_l @ vector_l)), vector_r / np.sqrt(vector_r @ vector_r)
    183            newcontent = [None if _check_for_none(self, item) else np.asarray([vector_l.T @ item @ vector_r]) for item in self.content]
    184
    185        else:
    186            # There are no checks here yet. There are so many possible scenarios, where this can go wrong.
    187            if normalize:
    188                for t in range(self.T):
    189                    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])
    190
    191            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)]
    192        return Corr(newcontent)

We need to project the Correlator with a Vector to get a single value at each timeslice.

The method can use one or two vectors. If two are specified it returns v1@G@v2 (the order might be very important.) By default it will return the lowest source, which usually means unsmeared-unsmeared (0,0), but it does not have to

def item(self, i, j): View Source

    194    def item(self, i, j):
    195        """Picks the element [i,j] from every matrix and returns a correlator containing one Obs per timeslice.
    196
    197        Parameters
    198        ----------
    199        i : int
    200            First index to be picked.
    201        j : int
    202            Second index to be picked.
    203        """
    204        if self.N == 1:
    205            raise Exception("Trying to pick item from projected Corr")
    206        newcontent = [None if (item is None) else item[i, j] for item in self.content]
    207        return Corr(newcontent)

Picks the element [i,j] from every matrix and returns a correlator containing one Obs per timeslice.

###### Parameters

- i (int): First index to be picked.
- j (int): Second index to be picked.

def plottable(self): View Source

    209    def plottable(self):
    210        """Outputs the correlator in a plotable format.
    211
    212        Outputs three lists containing the timeslice index, the value on each
    213        timeslice and the error on each timeslice.
    214        """
    215        if self.N != 1:
    216            raise Exception("Can only make Corr[N=1] plottable")
    217        x_list = [x for x in range(self.T) if not self.content[x] is None]
    218        y_list = [y[0].value for y in self.content if y is not None]
    219        y_err_list = [y[0].dvalue for y in self.content if y is not None]
    220
    221        return x_list, y_list, y_err_list

Outputs the correlator in a plotable format.

Outputs three lists containing the timeslice index, the value on each timeslice and the error on each timeslice.

def symmetric(self): View Source

    223    def symmetric(self):
    224        """ Symmetrize the correlator around x0=0."""
    225        if self.N != 1:
    226            raise Exception('symmetric cannot be safely applied to multi-dimensional correlators.')
    227        if self.T % 2 != 0:
    228            raise Exception("Can not symmetrize odd T")
    229
    230        if self.content[0] is not None:
    231            if np.argmax(np.abs([o[0].value if o is not None else 0 for o in self.content])) != 0:
    232                warnings.warn("Correlator does not seem to be symmetric around x0=0.", RuntimeWarning)
    233
    234        newcontent = [self.content[0]]
    235        for t in range(1, self.T):
    236            if (self.content[t] is None) or (self.content[self.T - t] is None):
    237                newcontent.append(None)
    238            else:
    239                newcontent.append(0.5 * (self.content[t] + self.content[self.T - t]))
    240        if (all([x is None for x in newcontent])):
    241            raise Exception("Corr could not be symmetrized: No redundant values")
    242        return Corr(newcontent, prange=self.prange)

Symmetrize the correlator around x0=0.

def anti_symmetric(self): View Source

    244    def anti_symmetric(self):
    245        """Anti-symmetrize the correlator around x0=0."""
    246        if self.N != 1:
    247            raise TypeError('anti_symmetric cannot be safely applied to multi-dimensional correlators.')
    248        if self.T % 2 != 0:
    249            raise Exception("Can not symmetrize odd T")
    250
    251        test = 1 * self
    252        test.gamma_method()
    253        if not all([o.is_zero_within_error(3) for o in test.content[0]]):
    254            warnings.warn("Correlator does not seem to be anti-symmetric around x0=0.", RuntimeWarning)
    255
    256        newcontent = [self.content[0]]
    257        for t in range(1, self.T):
    258            if (self.content[t] is None) or (self.content[self.T - t] is None):
    259                newcontent.append(None)
    260            else:
    261                newcontent.append(0.5 * (self.content[t] - self.content[self.T - t]))
    262        if (all([x is None for x in newcontent])):
    263            raise Exception("Corr could not be symmetrized: No redundant values")
    264        return Corr(newcontent, prange=self.prange)

Anti-symmetrize the correlator around x0=0.

def is_matrix_symmetric(self): View Source

    266    def is_matrix_symmetric(self):
    267        """Checks whether a correlator matrices is symmetric on every timeslice."""
    268        if self.N == 1:
    269            raise TypeError("Only works for correlator matrices.")
    270        for t in range(self.T):
    271            if self[t] is None:
    272                continue
    273            for i in range(self.N):
    274                for j in range(i + 1, self.N):
    275                    if self[t][i, j] is self[t][j, i]:
    276                        continue
    277                    if hash(self[t][i, j]) != hash(self[t][j, i]):
    278                        return False
    279        return True

Checks whether a correlator matrices is symmetric on every timeslice.

def trace(self): View Source

    281    def trace(self):
    282        """Calculates the per-timeslice trace of a correlator matrix."""
    283        if self.N == 1:
    284            raise ValueError("Only works for correlator matrices.")
    285        newcontent = []
    286        for t in range(self.T):
    287            if _check_for_none(self, self.content[t]):
    288                newcontent.append(None)
    289            else:
    290                newcontent.append(np.trace(self.content[t]))
    291        return Corr(newcontent)

Calculates the per-timeslice trace of a correlator matrix.

def matrix_symmetric(self): View Source

    293    def matrix_symmetric(self):
    294        """Symmetrizes the correlator matrices on every timeslice."""
    295        if self.N == 1:
    296            raise Exception("Trying to symmetrize a correlator matrix, that already has N=1.")
    297        if self.is_matrix_symmetric():
    298            return 1.0 * self
    299        else:
    300            transposed = [None if _check_for_none(self, G) else G.T for G in self.content]
    301            return 0.5 * (Corr(transposed) + self)

Symmetrizes the correlator matrices on every timeslice.

def GEVP(self, t0, ts=None, sort='Eigenvalue', vector_obs=False, **kwargs): View Source

    303    def GEVP(self, t0, ts=None, sort="Eigenvalue", vector_obs=False, **kwargs):
    304        r'''Solve the generalized eigenvalue problem on the correlator matrix and returns the corresponding eigenvectors.
    305
    306        The eigenvectors are sorted according to the descending eigenvalues, the zeroth eigenvector(s) correspond to the
    307        largest eigenvalue(s). The eigenvector(s) for the individual states can be accessed via slicing
    308        ```python
    309        C.GEVP(t0=2)[0]  # Ground state vector(s)
    310        C.GEVP(t0=2)[:3]  # Vectors for the lowest three states
    311        ```
    312
    313        Parameters
    314        ----------
    315        t0 : int
    316            The time t0 for the right hand side of the GEVP according to $G(t)v_i=\lambda_i G(t_0)v_i$
    317        ts : int
    318            fixed time $G(t_s)v_i=\lambda_i G(t_0)v_i$ if sort=None.
    319            If sort="Eigenvector" it gives a reference point for the sorting method.
    320        sort : string
    321            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.
    322            - "Eigenvalue": The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice. (default)
    323            - "Eigenvector": Use the method described in arXiv:2004.10472 to find the set of v(t) belonging to the state.
    324              The reference state is identified by its eigenvalue at $t=t_s$.
    325            - None: The GEVP is solved only at ts, no sorting is necessary
    326        vector_obs : bool
    327            If True, uncertainties are propagated in the eigenvector computation (default False).
    328
    329        Other Parameters
    330        ----------------
    331        state : int
    332           Returns only the vector(s) for a specified state. The lowest state is zero.
    333        method : str
    334           Method used to solve the GEVP.
    335           - "eigh": Use scipy.linalg.eigh to solve the GEVP. (default for vector_obs=False)
    336           - "cholesky": Use manually implemented solution via the Cholesky decomposition. Automatically chosen if vector_obs==True.
    337        '''
    338
    339        if self.N == 1:
    340            raise Exception("GEVP methods only works on correlator matrices and not single correlators.")
    341        if ts is not None:
    342            if (ts <= t0):
    343                raise Exception("ts has to be larger than t0.")
    344
    345        if "sorted_list" in kwargs:
    346            warnings.warn("Argument 'sorted_list' is deprecated, use 'sort' instead.", DeprecationWarning)
    347            sort = kwargs.get("sorted_list")
    348
    349        if self.is_matrix_symmetric():
    350            symmetric_corr = self
    351        else:
    352            symmetric_corr = self.matrix_symmetric()
    353
    354        def _get_mat_at_t(t, vector_obs=vector_obs):
    355            if vector_obs:
    356                return symmetric_corr[t]
    357            else:
    358                return np.vectorize(lambda x: x.value)(symmetric_corr[t])
    359        G0 = _get_mat_at_t(t0)
    360
    361        method = kwargs.get('method', 'eigh')
    362        if vector_obs:
    363            chol = linalg.cholesky(G0)
    364            chol_inv = linalg.inv(chol)
    365            method = 'cholesky'
    366        else:
    367            chol = np.linalg.cholesky(_get_mat_at_t(t0, vector_obs=False))  # Check if matrix G0 is positive-semidefinite.
    368            if method == 'cholesky':
    369                chol_inv = np.linalg.inv(chol)
    370            else:
    371                chol_inv = None
    372
    373        if sort is None:
    374            if (ts is None):
    375                raise Exception("ts is required if sort=None.")
    376            if (self.content[t0] is None) or (self.content[ts] is None):
    377                raise Exception("Corr not defined at t0/ts.")
    378            Gt = _get_mat_at_t(ts)
    379            reordered_vecs = _GEVP_solver(Gt, G0, method=method, chol_inv=chol_inv)
    380            if kwargs.get('auto_gamma', False) and vector_obs:
    381                [[o.gm() for o in ev if isinstance(o, Obs)] for ev in reordered_vecs]
    382
    383        elif sort in ["Eigenvalue", "Eigenvector"]:
    384            if sort == "Eigenvalue" and ts is not None:
    385                warnings.warn("ts has no effect when sorting by eigenvalue is chosen.", RuntimeWarning)
    386            all_vecs = [None] * (t0 + 1)
    387            for t in range(t0 + 1, self.T):
    388                try:
    389                    Gt = _get_mat_at_t(t)
    390                    all_vecs.append(_GEVP_solver(Gt, G0, method=method, chol_inv=chol_inv))
    391                except Exception:
    392                    all_vecs.append(None)
    393            if sort == "Eigenvector":
    394                if ts is None:
    395                    raise Exception("ts is required for the Eigenvector sorting method.")
    396                all_vecs = _sort_vectors(all_vecs, ts)
    397
    398            reordered_vecs = [[v[s] if v is not None else None for v in all_vecs] for s in range(self.N)]
    399            if kwargs.get('auto_gamma', False) and vector_obs:
    400                [[[o.gm() for o in evn] for evn in ev if evn is not None] for ev in reordered_vecs]
    401        else:
    402            raise Exception("Unknown value for 'sort'. Choose 'Eigenvalue', 'Eigenvector' or None.")
    403
    404        if "state" in kwargs:
    405            return reordered_vecs[kwargs.get("state")]
    406        else:
    407            return reordered_vecs

Solve the generalized eigenvalue problem on the correlator matrix and returns the corresponding eigenvectors.

The eigenvectors are sorted according to the descending eigenvalues, the zeroth eigenvector(s) correspond to the largest eigenvalue(s). The eigenvector(s) for the individual states can be accessed via slicing

    C.GEVP(t0=2)[0]  # Ground state vector(s)
    C.GEVP(t0=2)[:3]  # Vectors for the lowest three states

###### Parameters

- t0 (int): The time t0 for the right hand side of the GEVP according to $G(t)v_i=\lambda_i G(t_0)v_i$
- ts (int): fixed time $G(t_s)v_i=\lambda_i G(t_0)v_i$ if sort=None. If sort="Eigenvector" it gives a reference point for the sorting method.
- sort (string): 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.
  - "Eigenvalue": The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice. (default)
  - "Eigenvector": Use the method described in arXiv:2004.10472 to find the set of v(t) belonging to the state. The reference state is identified by its eigenvalue at $t=t_s$.
  - None: The GEVP is solved only at ts, no sorting is necessary
- vector_obs (bool): If True, uncertainties are propagated in the eigenvector computation (default False).

###### Other Parameters

- state (int): Returns only the vector(s) for a specified state. The lowest state is zero.
- method (str): Method used to solve the GEVP.
  - "eigh": Use scipy.linalg.eigh to solve the GEVP. (default for vector_obs=False)
  - "cholesky": Use manually implemented solution via the Cholesky decomposition. Automatically chosen if vector_obs==True.

def Eigenvalue(self, t0, ts=None, state=0, sort='Eigenvalue', **kwargs): View Source

    409    def Eigenvalue(self, t0, ts=None, state=0, sort="Eigenvalue", **kwargs):
    410        """Determines the eigenvalue of the GEVP by solving and projecting the correlator
    411
    412        Parameters
    413        ----------
    414        state : int
    415            The state one is interested in ordered by energy. The lowest state is zero.
    416
    417        All other parameters are identical to the ones of Corr.GEVP.
    418        """
    419        vec = self.GEVP(t0, ts=ts, sort=sort, **kwargs)[state]
    420        return self.projected(vec)

Determines the eigenvalue of the GEVP by solving and projecting the correlator

###### Parameters

- state (int): The state one is interested in ordered by energy. The lowest state is zero.
- All other parameters are identical to the ones of Corr.GEVP.

def Hankel(self, N, periodic=False): View Source

    422    def Hankel(self, N, periodic=False):
    423        """Constructs an NxN Hankel matrix
    424
    425        C(t) c(t+1) ... c(t+n-1)
    426        C(t+1) c(t+2) ... c(t+n)
    427        .................
    428        C(t+(n-1)) c(t+n) ... c(t+2(n-1))
    429
    430        Parameters
    431        ----------
    432        N : int
    433            Dimension of the Hankel matrix
    434        periodic : bool, optional
    435            determines whether the matrix is extended periodically
    436        """
    437
    438        if self.N != 1:
    439            raise Exception("Multi-operator Prony not implemented!")
    440
    441        array = np.empty([N, N], dtype="object")
    442        new_content = []
    443        for t in range(self.T):
    444            new_content.append(array.copy())
    445
    446        def wrap(i):
    447            while i >= self.T:
    448                i -= self.T
    449            return i
    450
    451        for t in range(self.T):
    452            for i in range(N):
    453                for j in range(N):
    454                    if periodic:
    455                        new_content[t][i, j] = self.content[wrap(t + i + j)][0]
    456                    elif (t + i + j) >= self.T:
    457                        new_content[t] = None
    458                    else:
    459                        new_content[t][i, j] = self.content[t + i + j][0]
    460
    461        return Corr(new_content)

Constructs an NxN Hankel matrix

C(t) c(t+1) ... c(t+n-1) C(t+1) c(t+2) ... c(t+n) ................. C(t+(n-1)) c(t+n) ... c(t+2(n-1))

###### Parameters

- N (int): Dimension of the Hankel matrix
- periodic (bool, optional): determines whether the matrix is extended periodically

def roll(self, dt): View Source

    463    def roll(self, dt):
    464        """Periodically shift the correlator by dt timeslices
    465
    466        Parameters
    467        ----------
    468        dt : int
    469            number of timeslices
    470        """
    471        return Corr(list(np.roll(np.array(self.content, dtype=object), dt, axis=0)))

Periodically shift the correlator by dt timeslices

###### Parameters

- dt (int): number of timeslices

def reverse(self): View Source

    473    def reverse(self):
    474        """Reverse the time ordering of the Corr"""
    475        return Corr(self.content[:: -1])

Reverse the time ordering of the Corr

def thin(self, spacing=2, offset=0): View Source

    477    def thin(self, spacing=2, offset=0):
    478        """Thin out a correlator to suppress correlations
    479
    480        Parameters
    481        ----------
    482        spacing : int
    483            Keep only every 'spacing'th entry of the correlator
    484        offset : int
    485            Offset the equal spacing
    486        """
    487        new_content = []
    488        for t in range(self.T):
    489            if (offset + t) % spacing != 0:
    490                new_content.append(None)
    491            else:
    492                new_content.append(self.content[t])
    493        return Corr(new_content)

Thin out a correlator to suppress correlations

###### Parameters

- spacing (int): Keep only every 'spacing'th entry of the correlator
- offset (int): Offset the equal spacing

def correlate(self, partner): View Source

    495    def correlate(self, partner):
    496        """Correlate the correlator with another correlator or Obs
    497
    498        Parameters
    499        ----------
    500        partner : Obs or Corr
    501            partner to correlate the correlator with.
    502            Can either be an Obs which is correlated with all entries of the
    503            correlator or a Corr of same length.
    504        """
    505        if self.N != 1:
    506            raise Exception("Only one-dimensional correlators can be safely correlated.")
    507        new_content = []
    508        for x0, t_slice in enumerate(self.content):
    509            if _check_for_none(self, t_slice):
    510                new_content.append(None)
    511            else:
    512                if isinstance(partner, Corr):
    513                    if _check_for_none(partner, partner.content[x0]):
    514                        new_content.append(None)
    515                    else:
    516                        new_content.append(np.array([correlate(o, partner.content[x0][0]) for o in t_slice]))
    517                elif isinstance(partner, Obs):  # Should this include CObs?
    518                    new_content.append(np.array([correlate(o, partner) for o in t_slice]))
    519                else:
    520                    raise Exception("Can only correlate with an Obs or a Corr.")
    521
    522        return Corr(new_content)

Correlate the correlator with another correlator or Obs

###### Parameters

- partner (Obs or Corr): partner to correlate the correlator with. Can either be an Obs which is correlated with all entries of the correlator or a Corr of same length.

def reweight(self, weight, **kwargs): View Source

    524    def reweight(self, weight, **kwargs):
    525        """Reweight the correlator.
    526
    527        Parameters
    528        ----------
    529        weight : Obs
    530            Reweighting factor. An Observable that has to be defined on a superset of the
    531            configurations in obs[i].idl for all i.
    532        all_configs : bool
    533            if True, the reweighted observables are normalized by the average of
    534            the reweighting factor on all configurations in weight.idl and not
    535            on the configurations in obs[i].idl.
    536        """
    537        if self.N != 1:
    538            raise Exception("Reweighting only implemented for one-dimensional correlators.")
    539        new_content = []
    540        for t_slice in self.content:
    541            if _check_for_none(self, t_slice):
    542                new_content.append(None)
    543            else:
    544                new_content.append(np.array(reweight(weight, t_slice, **kwargs)))
    545        return Corr(new_content)

Reweight the correlator.

###### Parameters

- weight (Obs): Reweighting factor. An Observable that has to be defined on a superset of the configurations in obs[i].idl for all i.
- all_configs (bool): if True, the reweighted observables are normalized by the average of the reweighting factor on all configurations in weight.idl and not on the configurations in obs[i].idl.

def T_symmetry(self, partner, parity=1): View Source

    547    def T_symmetry(self, partner, parity=+1):
    548        """Return the time symmetry average of the correlator and its partner
    549
    550        Parameters
    551        ----------
    552        partner : Corr
    553            Time symmetry partner of the Corr
    554        parity : int
    555            Parity quantum number of the correlator, can be +1 or -1
    556        """
    557        if self.N != 1:
    558            raise Exception("T_symmetry only implemented for one-dimensional correlators.")
    559        if not isinstance(partner, Corr):
    560            raise Exception("T partner has to be a Corr object.")
    561        if parity not in [+1, -1]:
    562            raise Exception("Parity has to be +1 or -1.")
    563        T_partner = parity * partner.reverse()
    564
    565        t_slices = []
    566        test = (self - T_partner)
    567        test.gamma_method()
    568        for x0, t_slice in enumerate(test.content):
    569            if t_slice is not None:
    570                if not t_slice[0].is_zero_within_error(5):
    571                    t_slices.append(x0)
    572        if t_slices:
    573            warnings.warn("T symmetry partners do not agree within 5 sigma on time slices " + str(t_slices) + ".", RuntimeWarning)
    574
    575        return (self + T_partner) / 2

Return the time symmetry average of the correlator and its partner

###### Parameters

- partner (Corr): Time symmetry partner of the Corr
- parity (int): Parity quantum number of the correlator, can be +1 or -1

def deriv(self, variant='symmetric'): View Source

    577    def deriv(self, variant="symmetric"):
    578        """Return the first derivative of the correlator with respect to x0.
    579
    580        Parameters
    581        ----------
    582        variant : str
    583            decides which definition of the finite differences derivative is used.
    584            Available choice: symmetric, forward, backward, improved, log, default: symmetric
    585        """
    586        if self.N != 1:
    587            raise Exception("deriv only implemented for one-dimensional correlators.")
    588        if variant == "symmetric":
    589            newcontent = []
    590            for t in range(1, self.T - 1):
    591                if (self.content[t - 1] is None) or (self.content[t + 1] is None):
    592                    newcontent.append(None)
    593                else:
    594                    newcontent.append(0.5 * (self.content[t + 1] - self.content[t - 1]))
    595            if (all([x is None for x in newcontent])):
    596                raise Exception('Derivative is undefined at all timeslices')
    597            return Corr(newcontent, padding=[1, 1])
    598        elif variant == "forward":
    599            newcontent = []
    600            for t in range(self.T - 1):
    601                if (self.content[t] is None) or (self.content[t + 1] is None):
    602                    newcontent.append(None)
    603                else:
    604                    newcontent.append(self.content[t + 1] - self.content[t])
    605            if (all([x is None for x in newcontent])):
    606                raise Exception("Derivative is undefined at all timeslices")
    607            return Corr(newcontent, padding=[0, 1])
    608        elif variant == "backward":
    609            newcontent = []
    610            for t in range(1, self.T):
    611                if (self.content[t - 1] is None) or (self.content[t] is None):
    612                    newcontent.append(None)
    613                else:
    614                    newcontent.append(self.content[t] - self.content[t - 1])
    615            if (all([x is None for x in newcontent])):
    616                raise Exception("Derivative is undefined at all timeslices")
    617            return Corr(newcontent, padding=[1, 0])
    618        elif variant == "improved":
    619            newcontent = []
    620            for t in range(2, self.T - 2):
    621                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):
    622                    newcontent.append(None)
    623                else:
    624                    newcontent.append((1 / 12) * (self.content[t - 2] - 8 * self.content[t - 1] + 8 * self.content[t + 1] - self.content[t + 2]))
    625            if (all([x is None for x in newcontent])):
    626                raise Exception('Derivative is undefined at all timeslices')
    627            return Corr(newcontent, padding=[2, 2])
    628        elif variant == 'log':
    629            newcontent = []
    630            for t in range(self.T):
    631                if (self.content[t] is None) or (self.content[t] <= 0):
    632                    newcontent.append(None)
    633                else:
    634                    newcontent.append(np.log(self.content[t]))
    635            if (all([x is None for x in newcontent])):
    636                raise Exception("Log is undefined at all timeslices")
    637            logcorr = Corr(newcontent)
    638            return self * logcorr.deriv('symmetric')
    639        else:
    640            raise Exception("Unknown variant.")

Return the first derivative of the correlator with respect to x0.

###### Parameters

- variant (str): decides which definition of the finite differences derivative is used. Available choice: symmetric, forward, backward, improved, log, default: symmetric

def second_deriv(self, variant='symmetric'): View Source

    642    def second_deriv(self, variant="symmetric"):
    643        r"""Return the second derivative of the correlator with respect to x0.
    644
    645        Parameters
    646        ----------
    647        variant : str
    648            decides which definition of the finite differences derivative is used.
    649            Available choice:
    650                - symmetric (default)
    651                    $$\tilde{\partial}^2_0 f(x_0) = f(x_0+1)-2f(x_0)+f(x_0-1)$$
    652                - big_symmetric
    653                    $$\partial^2_0 f(x_0) = \frac{f(x_0+2)-2f(x_0)+f(x_0-2)}{4}$$
    654                - improved
    655                    $$\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}$$
    656                - log
    657                    $$f(x) = \tilde{\partial}^2_0 log(f(x_0))+(\tilde{\partial}_0 log(f(x_0)))^2$$
    658        """
    659        if self.N != 1:
    660            raise Exception("second_deriv only implemented for one-dimensional correlators.")
    661        if variant == "symmetric":
    662            newcontent = []
    663            for t in range(1, self.T - 1):
    664                if (self.content[t - 1] is None) or (self.content[t + 1] is None):
    665                    newcontent.append(None)
    666                else:
    667                    newcontent.append((self.content[t + 1] - 2 * self.content[t] + self.content[t - 1]))
    668            if (all([x is None for x in newcontent])):
    669                raise Exception("Derivative is undefined at all timeslices")
    670            return Corr(newcontent, padding=[1, 1])
    671        elif variant == "big_symmetric":
    672            newcontent = []
    673            for t in range(2, self.T - 2):
    674                if (self.content[t - 2] is None) or (self.content[t + 2] is None):
    675                    newcontent.append(None)
    676                else:
    677                    newcontent.append((self.content[t + 2] - 2 * self.content[t] + self.content[t - 2]) / 4)
    678            if (all([x is None for x in newcontent])):
    679                raise Exception("Derivative is undefined at all timeslices")
    680            return Corr(newcontent, padding=[2, 2])
    681        elif variant == "improved":
    682            newcontent = []
    683            for t in range(2, self.T - 2):
    684                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):
    685                    newcontent.append(None)
    686                else:
    687                    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]))
    688            if (all([x is None for x in newcontent])):
    689                raise Exception("Derivative is undefined at all timeslices")
    690            return Corr(newcontent, padding=[2, 2])
    691        elif variant == 'log':
    692            newcontent = []
    693            for t in range(self.T):
    694                if (self.content[t] is None) or (self.content[t] <= 0):
    695                    newcontent.append(None)
    696                else:
    697                    newcontent.append(np.log(self.content[t]))
    698            if (all([x is None for x in newcontent])):
    699                raise Exception("Log is undefined at all timeslices")
    700            logcorr = Corr(newcontent)
    701            return self * (logcorr.second_deriv('symmetric') + (logcorr.deriv('symmetric'))**2)
    702        else:
    703            raise Exception("Unknown variant.")

Return the second derivative of the correlator with respect to x0.

###### Parameters

- variant (str): decides which definition of the finite differences derivative is used. Available choice: - symmetric (default) $$\tilde{\partial}^2_0 f(x_0) = f(x_0+1)-2f(x_0)+f(x_0-1)$$ - big_symmetric $$\partial^2_0 f(x_0) = \frac{f(x_0+2)-2f(x_0)+f(x_0-2)}{4}$$ - improved $$\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}$$ - log $$f(x) = \tilde{\partial}^2_0 log(f(x_0))+(\tilde{\partial}_0 log(f(x_0)))^2$$

def m_eff(self, variant='log', guess=1.0): View Source

    705    def m_eff(self, variant='log', guess=1.0):
    706        """Returns the effective mass of the correlator as correlator object
    707
    708        Parameters
    709        ----------
    710        variant : str
    711            log : uses the standard effective mass log(C(t) / C(t+1))
    712            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.
    713            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.
    714            See, e.g., arXiv:1205.5380
    715            arccosh : Uses the explicit form of the symmetrized correlator (not recommended)
    716            logsym: uses the symmetric effective mass log(C(t-1) / C(t+1))/2
    717        guess : float
    718            guess for the root finder, only relevant for the root variant
    719        """
    720        if self.N != 1:
    721            raise Exception('Correlator must be projected before getting m_eff')
    722        if variant == 'log':
    723            newcontent = []
    724            for t in range(self.T - 1):
    725                if ((self.content[t] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0):
    726                    newcontent.append(None)
    727                elif self.content[t][0].value / self.content[t + 1][0].value < 0:
    728                    newcontent.append(None)
    729                else:
    730                    newcontent.append(self.content[t] / self.content[t + 1])
    731            if (all([x is None for x in newcontent])):
    732                raise Exception('m_eff is undefined at all timeslices')
    733
    734            return np.log(Corr(newcontent, padding=[0, 1]))
    735
    736        elif variant == 'logsym':
    737            newcontent = []
    738            for t in range(1, self.T - 1):
    739                if ((self.content[t - 1] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0):
    740                    newcontent.append(None)
    741                elif self.content[t - 1][0].value / self.content[t + 1][0].value < 0:
    742                    newcontent.append(None)
    743                else:
    744                    newcontent.append(self.content[t - 1] / self.content[t + 1])
    745            if (all([x is None for x in newcontent])):
    746                raise Exception('m_eff is undefined at all timeslices')
    747
    748            return np.log(Corr(newcontent, padding=[1, 1])) / 2
    749
    750        elif variant in ['periodic', 'cosh', 'sinh']:
    751            if variant in ['periodic', 'cosh']:
    752                func = anp.cosh
    753            else:
    754                func = anp.sinh
    755
    756            def root_function(x, d):
    757                return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d
    758
    759            newcontent = []
    760            for t in range(self.T - 1):
    761                if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t + 1][0].value == 0):
    762                    newcontent.append(None)
    763                # Fill the two timeslices in the middle of the lattice with their predecessors
    764                elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]:
    765                    newcontent.append(newcontent[-1])
    766                elif self.content[t][0].value / self.content[t + 1][0].value < 0:
    767                    newcontent.append(None)
    768                else:
    769                    newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess)))
    770            if (all([x is None for x in newcontent])):
    771                raise Exception('m_eff is undefined at all timeslices')
    772
    773            return Corr(newcontent, padding=[0, 1])
    774
    775        elif variant == 'arccosh':
    776            newcontent = []
    777            for t in range(1, self.T - 1):
    778                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):
    779                    newcontent.append(None)
    780                else:
    781                    newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t]))
    782            if (all([x is None for x in newcontent])):
    783                raise Exception("m_eff is undefined at all timeslices")
    784            return np.arccosh(Corr(newcontent, padding=[1, 1]))
    785
    786        else:
    787            raise Exception('Unknown variant.')

Returns the effective mass of the correlator as correlator object

###### Parameters

- variant (str): log : uses the standard effective mass log(C(t) / C(t+1)) 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. 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. See, e.g., arXiv:1205.5380 arccosh : Uses the explicit form of the symmetrized correlator (not recommended) logsym: uses the symmetric effective mass log(C(t-1) / C(t+1))/2
- guess (float): guess for the root finder, only relevant for the root variant

def fit(self, function, fitrange=None, silent=False, **kwargs): View Source

    789    def fit(self, function, fitrange=None, silent=False, **kwargs):
    790        r'''Fits function to the data
    791
    792        Parameters
    793        ----------
    794        function : obj
    795            function to fit to the data. See fits.least_squares for details.
    796        fitrange : list
    797            Two element list containing the timeslices on which the fit is supposed to start and stop.
    798            Caution: This range is inclusive as opposed to standard python indexing.
    799            `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6.
    800            If not specified, self.prange or all timeslices are used.
    801        silent : bool
    802            Decides whether output is printed to the standard output.
    803        '''
    804        if self.N != 1:
    805            raise Exception("Correlator must be projected before fitting")
    806
    807        if fitrange is None:
    808            if self.prange:
    809                fitrange = self.prange
    810            else:
    811                fitrange = [0, self.T - 1]
    812        else:
    813            if not isinstance(fitrange, list):
    814                raise Exception("fitrange has to be a list with two elements")
    815            if len(fitrange) != 2:
    816                raise Exception("fitrange has to have exactly two elements [fit_start, fit_stop]")
    817
    818        xs = np.array([x for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None])
    819        ys = np.array([self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None])
    820        result = least_squares(xs, ys, function, silent=silent, **kwargs)
    821        return result

Fits function to the data

###### Parameters

- function (obj): function to fit to the data. See fits.least_squares for details.
- fitrange (list): Two element list containing the timeslices on which the fit is supposed to start and stop. Caution: This range is inclusive as opposed to standard python indexing. `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6. If not specified, self.prange or all timeslices are used.
- silent (bool): Decides whether output is printed to the standard output.

def plateau(self, plateau_range=None, method='fit', auto_gamma=False): View Source

    823    def plateau(self, plateau_range=None, method="fit", auto_gamma=False):
    824        """ Extract a plateau value from a Corr object
    825
    826        Parameters
    827        ----------
    828        plateau_range : list
    829            list with two entries, indicating the first and the last timeslice
    830            of the plateau region.
    831        method : str
    832            method to extract the plateau.
    833                'fit' fits a constant to the plateau region
    834                'avg', 'average' or 'mean' just average over the given timeslices.
    835        auto_gamma : bool
    836            apply gamma_method with default parameters to the Corr. Defaults to None
    837        """
    838        if not plateau_range:
    839            if self.prange:
    840                plateau_range = self.prange
    841            else:
    842                raise Exception("no plateau range provided")
    843        if self.N != 1:
    844            raise Exception("Correlator must be projected before getting a plateau.")
    845        if (all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])):
    846            raise Exception("plateau is undefined at all timeslices in plateaurange.")
    847        if auto_gamma:
    848            self.gamma_method()
    849        if method == "fit":
    850            def const_func(a, t):
    851                return a[0]
    852            return self.fit(const_func, plateau_range)[0]
    853        elif method in ["avg", "average", "mean"]:
    854            returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None])
    855            return returnvalue
    856
    857        else:
    858            raise Exception("Unsupported plateau method: " + method)

Extract a plateau value from a Corr object

###### Parameters

- plateau_range (list): list with two entries, indicating the first and the last timeslice of the plateau region.
- method (str): method to extract the plateau. 'fit' fits a constant to the plateau region 'avg', 'average' or 'mean' just average over the given timeslices.
- auto_gamma (bool): apply gamma_method with default parameters to the Corr. Defaults to None

def set_prange(self, prange): View Source

    860    def set_prange(self, prange):
    861        """Sets the attribute prange of the Corr object."""
    862        if not len(prange) == 2:
    863            raise Exception("prange must be a list or array with two values")
    864        if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))):
    865            raise Exception("Start and end point must be integers")
    866        if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] < prange[1]):
    867            raise Exception("Start and end point must define a range in the interval 0,T")
    868
    869        self.prange = prange
    870        return

Sets the attribute prange of the Corr object.

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): View Source

    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 Exception("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 Exception("'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 Exception("'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 Exception("'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 Exception("'save' has to be a string.")

Plots the correlator using the tag of the correlator as label if available.

###### Parameters

- x_range (list): list of two values, determining the range of the x-axis e.g. [4, 8].
- comp (Corr or list of Corr): Correlator or list of correlators which are plotted for comparison. The tags of these correlators are used as labels if available.
- logscale (bool): Sets y-axis to logscale.
- plateau (Obs): Plateau value to be visualized in the figure.
- fit_res (Fit_result): Fit_result object to be visualized.
- fit_key (str): Key for the fit function in Fit_result.fit_function (for combined fits).
- ylabel (str): Label for the y-axis.
- save (str): path to file in which the figure should be saved.
- auto_gamma (bool): Apply the gamma method with standard parameters to all correlators and plateau values before plotting.
- hide_sigma (float): Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors.
- references (list): List of floating point values that are displayed as horizontal lines for reference.
- title (string): Optional title of the figure.

def spaghetti_plot(self, logscale=True): View Source

     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 Exception("Correlator needs to be projected first.")
    1007
    1008        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]))
    1009        x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) 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()

Produces a spaghetti plot of the correlator suited to monitor exceptional configurations.

###### Parameters

- logscale (bool): Determines whether the scale of the y-axis is logarithmic or standard.

def dump(self, filename, datatype='json.gz', **kwargs): View Source

    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 Exception("Unknown datatype " + str(datatype))

Dumps the Corr into a file of chosen type

###### Parameters

- filename (str): Name of the file to be saved.
- datatype (str): Format of the exported file. Supported formats include "json.gz" and "pickle"
- path (str): specifies a custom path for the file (default '.')

def print(self, print_range=None): View Source

    1050    def print(self, print_range=None):
    1051        print(self.__repr__(print_range))

def sqrt(self): View Source

    1267    def sqrt(self):
    1268        return self ** 0.5

def log(self): View Source

    1270    def log(self):
    1271        newcontent = [None if _check_for_none(self, item) else np.log(item) for item in self.content]
    1272        return Corr(newcontent, prange=self.prange)

def exp(self): View Source

    1274    def exp(self):
    1275        newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content]
    1276        return Corr(newcontent, prange=self.prange)

def sin(self): View Source

    1291    def sin(self):
    1292        return self._apply_func_to_corr(np.sin)

def cos(self): View Source

    1294    def cos(self):
    1295        return self._apply_func_to_corr(np.cos)

def tan(self): View Source

    1297    def tan(self):
    1298        return self._apply_func_to_corr(np.tan)

def sinh(self): View Source

    1300    def sinh(self):
    1301        return self._apply_func_to_corr(np.sinh)

def cosh(self): View Source

    1303    def cosh(self):
    1304        return self._apply_func_to_corr(np.cosh)

def tanh(self): View Source

    1306    def tanh(self):
    1307        return self._apply_func_to_corr(np.tanh)

def arcsin(self): View Source

    1309    def arcsin(self):
    1310        return self._apply_func_to_corr(np.arcsin)

def arccos(self): View Source

    1312    def arccos(self):
    1313        return self._apply_func_to_corr(np.arccos)

def arctan(self): View Source

    1315    def arctan(self):
    1316        return self._apply_func_to_corr(np.arctan)

def arcsinh(self): View Source

    1318    def arcsinh(self):
    1319        return self._apply_func_to_corr(np.arcsinh)

def arccosh(self): View Source

    1321    def arccosh(self):
    1322        return self._apply_func_to_corr(np.arccosh)

def arctanh(self): View Source

    1324    def arctanh(self):
    1325        return self._apply_func_to_corr(np.arctanh)

real View Source

    1340    @property
    1341    def real(self):
    1342        def return_real(obs_OR_cobs):
    1343            if isinstance(obs_OR_cobs.flatten()[0], CObs):
    1344                return np.vectorize(lambda x: x.real)(obs_OR_cobs)
    1345            else:
    1346                return obs_OR_cobs
    1347
    1348        return self._apply_func_to_corr(return_real)

imag View Source

    1350    @property
    1351    def imag(self):
    1352        def return_imag(obs_OR_cobs):
    1353            if isinstance(obs_OR_cobs.flatten()[0], CObs):
    1354                return np.vectorize(lambda x: x.imag)(obs_OR_cobs)
    1355            else:
    1356                return obs_OR_cobs * 0  # So it stays the right type
    1357
    1358        return self._apply_func_to_corr(return_imag)

def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None): View Source

    1360    def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
    1361        r''' Project large correlation matrix to lowest states
    1362
    1363        This method can be used to reduce the size of an (N x N) correlation matrix
    1364        to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise
    1365        is still small.
    1366
    1367        Parameters
    1368        ----------
    1369        Ntrunc: int
    1370            Rank of the target matrix.
    1371        tproj: int
    1372            Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method.
    1373            The default value is 3.
    1374        t0proj: int
    1375            Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly
    1376            discouraged for O(a) improved theories, since the correctness of the procedure
    1377            cannot be granted in this case. The default value is 2.
    1378        basematrix : Corr
    1379            Correlation matrix that is used to determine the eigenvectors of the
    1380            lowest states based on a GEVP. basematrix is taken to be the Corr itself if
    1381            is is not specified.
    1382
    1383        Notes
    1384        -----
    1385        We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving
    1386        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}$
    1387        and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the
    1388        resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via
    1389        $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large
    1390        correlation matrix and to remove some noise that is added by irrelevant operators.
    1391        This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated
    1392        bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$.
    1393        '''
    1394
    1395        if self.N == 1:
    1396            raise Exception('Method cannot be applied to one-dimensional correlators.')
    1397        if basematrix is None:
    1398            basematrix = self
    1399        if Ntrunc >= basematrix.N:
    1400            raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N))
    1401        if basematrix.N != self.N:
    1402            raise Exception('basematrix and targetmatrix have to be of the same size.')
    1403
    1404        evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc]
    1405
    1406        tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object)
    1407        rmat = []
    1408        for t in range(basematrix.T):
    1409            for i in range(Ntrunc):
    1410                for j in range(Ntrunc):
    1411                    tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j]
    1412            rmat.append(np.copy(tmpmat))
    1413
    1414        newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)]
    1415        return Corr(newcontent)

Project large correlation matrix to lowest states

This method can be used to reduce the size of an (N x N) correlation matrix to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise is still small.

###### Parameters

- Ntrunc (int): Rank of the target matrix.
- tproj (int): Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method. The default value is 3.
- t0proj (int): Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly discouraged for O(a) improved theories, since the correctness of the procedure cannot be granted in this case. The default value is 2.
- basematrix (Corr): Correlation matrix that is used to determine the eigenvectors of the lowest states based on a GEVP. basematrix is taken to be the Corr itself if is is not specified.

###### Notes

We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving 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}$ and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large correlation matrix and to remove some noise that is added by irrelevant operators. This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$.

N