# pyerrors.correlators API documentation

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

class Corr: View Source

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

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

    29    def __init__(self, data_input, padding=[0, 0], prange=None):
    30        """ Initialize a Corr object.
    31
    32        Parameters
    33        ----------
    34        data_input : list or array
    35            list of Obs or list of arrays of Obs or array of Corrs
    36        padding : list, optional
    37            List with two entries where the first labels the padding
    38            at the front of the correlator and the second the padding
    39            at the back.
    40        prange : list, optional
    41            List containing the first and last timeslice of the plateau
    42            region indentified for this correlator.
    43        """
    44
    45        if isinstance(data_input, np.ndarray):
    46
    47            # This only works, if the array fulfills the conditions below
    48            if not len(data_input.shape) == 2 and data_input.shape[0] == data_input.shape[1]:
    49                raise Exception("Incompatible array shape")
    50            if not all([isinstance(item, Corr) for item in data_input.flatten()]):
    51                raise Exception("If the input is an array, its elements must be of type pe.Corr")
    52            if not all([item.N == 1 for item in data_input.flatten()]):
    53                raise Exception("Can only construct matrix correlator from single valued correlators")
    54            if not len(set([item.T for item in data_input.flatten()])) == 1:
    55                raise Exception("All input Correlators must be defined over the same timeslices.")
    56
    57            T = data_input[0, 0].T
    58            N = data_input.shape[0]
    59            input_as_list = []
    60            for t in range(T):
    61                if any([(item.content[t] is None) for item in data_input.flatten()]):
    62                    if not all([(item.content[t] is None) for item in data_input.flatten()]):
    63                        warnings.warn("Input ill-defined at different timeslices. Conversion leads to data loss!", RuntimeWarning)
    64                    input_as_list.append(None)
    65                else:
    66                    array_at_timeslace = np.empty([N, N], dtype="object")
    67                    for i in range(N):
    68                        for j in range(N):
    69                            array_at_timeslace[i, j] = data_input[i, j][t]
    70                    input_as_list.append(array_at_timeslace)
    71            data_input = input_as_list
    72
    73        if isinstance(data_input, list):
    74
    75            if all([isinstance(item, (Obs, CObs)) or item is None for item in data_input]):
    76                _assert_equal_properties([o for o in data_input if o is not None])
    77                self.content = [np.asarray([item]) if item is not None else None for item in data_input]
    78                self.N = 1
    79
    80            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]):
    81                self.content = data_input
    82                noNull = [a for a in self.content if not (a is None)]  # To check if the matrices are correct for all undefined elements
    83                self.N = noNull[0].shape[0]
    84                if self.N > 1 and noNull[0].shape[0] != noNull[0].shape[1]:
    85                    raise Exception("Smearing matrices are not NxN")
    86                if (not all([item.shape == noNull[0].shape for item in noNull])):
    87                    raise Exception("Items in data_input are not of identical shape." + str(noNull))
    88            else:
    89                raise Exception("data_input contains item of wrong type")
    90        else:
    91            raise Exception("Data input was not given as list or correct array")
    92
    93        self.tag = None
    94
    95        # An undefined timeslice is represented by the None object
    96        self.content = [None] * padding[0] + self.content + [None] * padding[1]
    97        self.T = len(self.content)
    98        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
- 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 indentified for this correlator.

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

    119    def gamma_method(self, **kwargs):
    120        """Apply the gamma method to the content of the Corr."""
    121        for item in self.content:
    122            if not (item is None):
    123                if self.N == 1:
    124                    item[0].gamma_method(**kwargs)
    125                else:
    126                    for i in range(self.N):
    127                        for j in range(self.N):
    128                            item[i, j].gamma_method(**kwargs)

Apply the gamma method to the content of the Corr.

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

    119    def gamma_method(self, **kwargs):
    120        """Apply the gamma method to the content of the Corr."""
    121        for item in self.content:
    122            if not (item is None):
    123                if self.N == 1:
    124                    item[0].gamma_method(**kwargs)
    125                else:
    126                    for i in range(self.N):
    127                        for j in range(self.N):
    128                            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

    132    def projected(self, vector_l=None, vector_r=None, normalize=False):
    133        """We need to project the Correlator with a Vector to get a single value at each timeslice.
    134
    135        The method can use one or two vectors.
    136        If two are specified it returns v1@G@v2 (the order might be very important.)
    137        By default it will return the lowest source, which usually means unsmeared-unsmeared (0,0), but it does not have to
    138        """
    139        if self.N == 1:
    140            raise Exception("Trying to project a Corr, that already has N=1.")
    141
    142        if vector_l is None:
    143            vector_l, vector_r = np.asarray([1.] + (self.N - 1) * [0.]), np.asarray([1.] + (self.N - 1) * [0.])
    144        elif (vector_r is None):
    145            vector_r = vector_l
    146        if isinstance(vector_l, list) and not isinstance(vector_r, list):
    147            if len(vector_l) != self.T:
    148                raise Exception("Length of vector list must be equal to T")
    149            vector_r = [vector_r] * self.T
    150        if isinstance(vector_r, list) and not isinstance(vector_l, list):
    151            if len(vector_r) != self.T:
    152                raise Exception("Length of vector list must be equal to T")
    153            vector_l = [vector_l] * self.T
    154
    155        if not isinstance(vector_l, list):
    156            if not vector_l.shape == vector_r.shape == (self.N,):
    157                raise Exception("Vectors are of wrong shape!")
    158            if normalize:
    159                vector_l, vector_r = vector_l / np.sqrt((vector_l @ vector_l)), vector_r / np.sqrt(vector_r @ vector_r)
    160            newcontent = [None if _check_for_none(self, item) else np.asarray([vector_l.T @ item @ vector_r]) for item in self.content]
    161
    162        else:
    163            # There are no checks here yet. There are so many possible scenarios, where this can go wrong.
    164            if normalize:
    165                for t in range(self.T):
    166                    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])
    167
    168            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)]
    169        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

    171    def item(self, i, j):
    172        """Picks the element [i,j] from every matrix and returns a correlator containing one Obs per timeslice.
    173
    174        Parameters
    175        ----------
    176        i : int
    177            First index to be picked.
    178        j : int
    179            Second index to be picked.
    180        """
    181        if self.N == 1:
    182            raise Exception("Trying to pick item from projected Corr")
    183        newcontent = [None if (item is None) else item[i, j] for item in self.content]
    184        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

    186    def plottable(self):
    187        """Outputs the correlator in a plotable format.
    188
    189        Outputs three lists containing the timeslice index, the value on each
    190        timeslice and the error on each timeslice.
    191        """
    192        if self.N != 1:
    193            raise Exception("Can only make Corr[N=1] plottable")
    194        x_list = [x for x in range(self.T) if not self.content[x] is None]
    195        y_list = [y[0].value for y in self.content if y is not None]
    196        y_err_list = [y[0].dvalue for y in self.content if y is not None]
    197
    198        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

    200    def symmetric(self):
    201        """ Symmetrize the correlator around x0=0."""
    202        if self.N != 1:
    203            raise Exception('symmetric cannot be safely applied to multi-dimensional correlators.')
    204        if self.T % 2 != 0:
    205            raise Exception("Can not symmetrize odd T")
    206
    207        if self.content[0] is not None:
    208            if np.argmax(np.abs([o[0].value if o is not None else 0 for o in self.content])) != 0:
    209                warnings.warn("Correlator does not seem to be symmetric around x0=0.", RuntimeWarning)
    210
    211        newcontent = [self.content[0]]
    212        for t in range(1, self.T):
    213            if (self.content[t] is None) or (self.content[self.T - t] is None):
    214                newcontent.append(None)
    215            else:
    216                newcontent.append(0.5 * (self.content[t] + self.content[self.T - t]))
    217        if (all([x is None for x in newcontent])):
    218            raise Exception("Corr could not be symmetrized: No redundant values")
    219        return Corr(newcontent, prange=self.prange)

Symmetrize the correlator around x0=0.

def anti_symmetric(self): View Source

    221    def anti_symmetric(self):
    222        """Anti-symmetrize the correlator around x0=0."""
    223        if self.N != 1:
    224            raise Exception('anti_symmetric cannot be safely applied to multi-dimensional correlators.')
    225        if self.T % 2 != 0:
    226            raise Exception("Can not symmetrize odd T")
    227
    228        test = 1 * self
    229        test.gamma_method()
    230        if not all([o.is_zero_within_error(3) for o in test.content[0]]):
    231            warnings.warn("Correlator does not seem to be anti-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)

Anti-symmetrize the correlator around x0=0.

def is_matrix_symmetric(self): View Source

    243    def is_matrix_symmetric(self):
    244        """Checks whether a correlator matrices is symmetric on every timeslice."""
    245        if self.N == 1:
    246            raise Exception("Only works for correlator matrices.")
    247        for t in range(self.T):
    248            if self[t] is None:
    249                continue
    250            for i in range(self.N):
    251                for j in range(i + 1, self.N):
    252                    if self[t][i, j] is self[t][j, i]:
    253                        continue
    254                    if hash(self[t][i, j]) != hash(self[t][j, i]):
    255                        return False
    256        return True

Checks whether a correlator matrices is symmetric on every timeslice.

def matrix_symmetric(self): View Source

    258    def matrix_symmetric(self):
    259        """Symmetrizes the correlator matrices on every timeslice."""
    260        if self.N == 1:
    261            raise Exception("Trying to symmetrize a correlator matrix, that already has N=1.")
    262        if self.is_matrix_symmetric():
    263            return 1.0 * self
    264        else:
    265            transposed = [None if _check_for_none(self, G) else G.T for G in self.content]
    266            return 0.5 * (Corr(transposed) + self)

Symmetrizes the correlator matrices on every timeslice.

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

    268    def GEVP(self, t0, ts=None, sort="Eigenvalue", **kwargs):
    269        r'''Solve the generalized eigenvalue problem on the correlator matrix and returns the corresponding eigenvectors.
    270
    271        The eigenvectors are sorted according to the descending eigenvalues, the zeroth eigenvector(s) correspond to the
    272        largest eigenvalue(s). The eigenvector(s) for the individual states can be accessed via slicing
    273        ```python
    274        C.GEVP(t0=2)[0]  # Ground state vector(s)
    275        C.GEVP(t0=2)[:3]  # Vectors for the lowest three states
    276        ```
    277
    278        Parameters
    279        ----------
    280        t0 : int
    281            The time t0 for the right hand side of the GEVP according to $G(t)v_i=\lambda_i G(t_0)v_i$
    282        ts : int
    283            fixed time $G(t_s)v_i=\lambda_i G(t_0)v_i$ if sort=None.
    284            If sort="Eigenvector" it gives a reference point for the sorting method.
    285        sort : string
    286            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.
    287            - "Eigenvalue": The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice.
    288            - "Eigenvector": Use the method described in arXiv:2004.10472 to find the set of v(t) belonging to the state.
    289              The reference state is identified by its eigenvalue at $t=t_s$.
    290
    291        Other Parameters
    292        ----------------
    293        state : int
    294           Returns only the vector(s) for a specified state. The lowest state is zero.
    295        '''
    296
    297        if self.N == 1:
    298            raise Exception("GEVP methods only works on correlator matrices and not single correlators.")
    299        if ts is not None:
    300            if (ts <= t0):
    301                raise Exception("ts has to be larger than t0.")
    302
    303        if "sorted_list" in kwargs:
    304            warnings.warn("Argument 'sorted_list' is deprecated, use 'sort' instead.", DeprecationWarning)
    305            sort = kwargs.get("sorted_list")
    306
    307        if self.is_matrix_symmetric():
    308            symmetric_corr = self
    309        else:
    310            symmetric_corr = self.matrix_symmetric()
    311
    312        G0 = np.vectorize(lambda x: x.value)(symmetric_corr[t0])
    313        np.linalg.cholesky(G0)  # Check if matrix G0 is positive-semidefinite.
    314
    315        if sort is None:
    316            if (ts is None):
    317                raise Exception("ts is required if sort=None.")
    318            if (self.content[t0] is None) or (self.content[ts] is None):
    319                raise Exception("Corr not defined at t0/ts.")
    320            Gt = np.vectorize(lambda x: x.value)(symmetric_corr[ts])
    321            reordered_vecs = _GEVP_solver(Gt, G0)
    322
    323        elif sort in ["Eigenvalue", "Eigenvector"]:
    324            if sort == "Eigenvalue" and ts is not None:
    325                warnings.warn("ts has no effect when sorting by eigenvalue is chosen.", RuntimeWarning)
    326            all_vecs = [None] * (t0 + 1)
    327            for t in range(t0 + 1, self.T):
    328                try:
    329                    Gt = np.vectorize(lambda x: x.value)(symmetric_corr[t])
    330                    all_vecs.append(_GEVP_solver(Gt, G0))
    331                except Exception:
    332                    all_vecs.append(None)
    333            if sort == "Eigenvector":
    334                if ts is None:
    335                    raise Exception("ts is required for the Eigenvector sorting method.")
    336                all_vecs = _sort_vectors(all_vecs, ts)
    337
    338            reordered_vecs = [[v[s] if v is not None else None for v in all_vecs] for s in range(self.N)]
    339        else:
    340            raise Exception("Unkown value for 'sort'.")
    341
    342        if "state" in kwargs:
    343            return reordered_vecs[kwargs.get("state")]
    344        else:
    345            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.
  - "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$.

###### Other Parameters

- state (int): Returns only the vector(s) for a specified state. The lowest state is zero.

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

    347    def Eigenvalue(self, t0, ts=None, state=0, sort="Eigenvalue"):
    348        """Determines the eigenvalue of the GEVP by solving and projecting the correlator
    349
    350        Parameters
    351        ----------
    352        state : int
    353            The state one is interested in ordered by energy. The lowest state is zero.
    354
    355        All other parameters are identical to the ones of Corr.GEVP.
    356        """
    357        vec = self.GEVP(t0, ts=ts, sort=sort)[state]
    358        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

    360    def Hankel(self, N, periodic=False):
    361        """Constructs an NxN Hankel matrix
    362
    363        C(t) c(t+1) ... c(t+n-1)
    364        C(t+1) c(t+2) ... c(t+n)
    365        .................
    366        C(t+(n-1)) c(t+n) ... c(t+2(n-1))
    367
    368        Parameters
    369        ----------
    370        N : int
    371            Dimension of the Hankel matrix
    372        periodic : bool, optional
    373            determines whether the matrix is extended periodically
    374        """
    375
    376        if self.N != 1:
    377            raise Exception("Multi-operator Prony not implemented!")
    378
    379        array = np.empty([N, N], dtype="object")
    380        new_content = []
    381        for t in range(self.T):
    382            new_content.append(array.copy())
    383
    384        def wrap(i):
    385            while i >= self.T:
    386                i -= self.T
    387            return i
    388
    389        for t in range(self.T):
    390            for i in range(N):
    391                for j in range(N):
    392                    if periodic:
    393                        new_content[t][i, j] = self.content[wrap(t + i + j)][0]
    394                    elif (t + i + j) >= self.T:
    395                        new_content[t] = None
    396                    else:
    397                        new_content[t][i, j] = self.content[t + i + j][0]
    398
    399        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

    401    def roll(self, dt):
    402        """Periodically shift the correlator by dt timeslices
    403
    404        Parameters
    405        ----------
    406        dt : int
    407            number of timeslices
    408        """
    409        return Corr(list(np.roll(np.array(self.content, dtype=object), dt)))

Periodically shift the correlator by dt timeslices

###### Parameters

- dt (int): number of timeslices

def reverse(self): View Source

    411    def reverse(self):
    412        """Reverse the time ordering of the Corr"""
    413        return Corr(self.content[:: -1])

Reverse the time ordering of the Corr

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

    415    def thin(self, spacing=2, offset=0):
    416        """Thin out a correlator to suppress correlations
    417
    418        Parameters
    419        ----------
    420        spacing : int
    421            Keep only every 'spacing'th entry of the correlator
    422        offset : int
    423            Offset the equal spacing
    424        """
    425        new_content = []
    426        for t in range(self.T):
    427            if (offset + t) % spacing != 0:
    428                new_content.append(None)
    429            else:
    430                new_content.append(self.content[t])
    431        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

    433    def correlate(self, partner):
    434        """Correlate the correlator with another correlator or Obs
    435
    436        Parameters
    437        ----------
    438        partner : Obs or Corr
    439            partner to correlate the correlator with.
    440            Can either be an Obs which is correlated with all entries of the
    441            correlator or a Corr of same length.
    442        """
    443        if self.N != 1:
    444            raise Exception("Only one-dimensional correlators can be safely correlated.")
    445        new_content = []
    446        for x0, t_slice in enumerate(self.content):
    447            if _check_for_none(self, t_slice):
    448                new_content.append(None)
    449            else:
    450                if isinstance(partner, Corr):
    451                    if _check_for_none(partner, partner.content[x0]):
    452                        new_content.append(None)
    453                    else:
    454                        new_content.append(np.array([correlate(o, partner.content[x0][0]) for o in t_slice]))
    455                elif isinstance(partner, Obs):  # Should this include CObs?
    456                    new_content.append(np.array([correlate(o, partner) for o in t_slice]))
    457                else:
    458                    raise Exception("Can only correlate with an Obs or a Corr.")
    459
    460        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

    462    def reweight(self, weight, **kwargs):
    463        """Reweight the correlator.
    464
    465        Parameters
    466        ----------
    467        weight : Obs
    468            Reweighting factor. An Observable that has to be defined on a superset of the
    469            configurations in obs[i].idl for all i.
    470        all_configs : bool
    471            if True, the reweighted observables are normalized by the average of
    472            the reweighting factor on all configurations in weight.idl and not
    473            on the configurations in obs[i].idl.
    474        """
    475        if self.N != 1:
    476            raise Exception("Reweighting only implemented for one-dimensional correlators.")
    477        new_content = []
    478        for t_slice in self.content:
    479            if _check_for_none(self, t_slice):
    480                new_content.append(None)
    481            else:
    482                new_content.append(np.array(reweight(weight, t_slice, **kwargs)))
    483        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

    485    def T_symmetry(self, partner, parity=+1):
    486        """Return the time symmetry average of the correlator and its partner
    487
    488        Parameters
    489        ----------
    490        partner : Corr
    491            Time symmetry partner of the Corr
    492        partity : int
    493            Parity quantum number of the correlator, can be +1 or -1
    494        """
    495        if self.N != 1:
    496            raise Exception("T_symmetry only implemented for one-dimensional correlators.")
    497        if not isinstance(partner, Corr):
    498            raise Exception("T partner has to be a Corr object.")
    499        if parity not in [+1, -1]:
    500            raise Exception("Parity has to be +1 or -1.")
    501        T_partner = parity * partner.reverse()
    502
    503        t_slices = []
    504        test = (self - T_partner)
    505        test.gamma_method()
    506        for x0, t_slice in enumerate(test.content):
    507            if t_slice is not None:
    508                if not t_slice[0].is_zero_within_error(5):
    509                    t_slices.append(x0)
    510        if t_slices:
    511            warnings.warn("T symmetry partners do not agree within 5 sigma on time slices " + str(t_slices) + ".", RuntimeWarning)
    512
    513        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
- partity (int): Parity quantum number of the correlator, can be +1 or -1

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

    515    def deriv(self, variant="symmetric"):
    516        """Return the first derivative of the correlator with respect to x0.
    517
    518        Parameters
    519        ----------
    520        variant : str
    521            decides which definition of the finite differences derivative is used.
    522            Available choice: symmetric, forward, backward, improved, log, default: symmetric
    523        """
    524        if self.N != 1:
    525            raise Exception("deriv only implemented for one-dimensional correlators.")
    526        if variant == "symmetric":
    527            newcontent = []
    528            for t in range(1, self.T - 1):
    529                if (self.content[t - 1] is None) or (self.content[t + 1] is None):
    530                    newcontent.append(None)
    531                else:
    532                    newcontent.append(0.5 * (self.content[t + 1] - self.content[t - 1]))
    533            if (all([x is None for x in newcontent])):
    534                raise Exception('Derivative is undefined at all timeslices')
    535            return Corr(newcontent, padding=[1, 1])
    536        elif variant == "forward":
    537            newcontent = []
    538            for t in range(self.T - 1):
    539                if (self.content[t] is None) or (self.content[t + 1] is None):
    540                    newcontent.append(None)
    541                else:
    542                    newcontent.append(self.content[t + 1] - self.content[t])
    543            if (all([x is None for x in newcontent])):
    544                raise Exception("Derivative is undefined at all timeslices")
    545            return Corr(newcontent, padding=[0, 1])
    546        elif variant == "backward":
    547            newcontent = []
    548            for t in range(1, self.T):
    549                if (self.content[t - 1] is None) or (self.content[t] is None):
    550                    newcontent.append(None)
    551                else:
    552                    newcontent.append(self.content[t] - self.content[t - 1])
    553            if (all([x is None for x in newcontent])):
    554                raise Exception("Derivative is undefined at all timeslices")
    555            return Corr(newcontent, padding=[1, 0])
    556        elif variant == "improved":
    557            newcontent = []
    558            for t in range(2, self.T - 2):
    559                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):
    560                    newcontent.append(None)
    561                else:
    562                    newcontent.append((1 / 12) * (self.content[t - 2] - 8 * self.content[t - 1] + 8 * self.content[t + 1] - self.content[t + 2]))
    563            if (all([x is None for x in newcontent])):
    564                raise Exception('Derivative is undefined at all timeslices')
    565            return Corr(newcontent, padding=[2, 2])
    566        elif variant == 'log':
    567            newcontent = []
    568            for t in range(self.T):
    569                if (self.content[t] is None) or (self.content[t] <= 0):
    570                    newcontent.append(None)
    571                else:
    572                    newcontent.append(np.log(self.content[t]))
    573            if (all([x is None for x in newcontent])):
    574                raise Exception("Log is undefined at all timeslices")
    575            logcorr = Corr(newcontent)
    576            return self * logcorr.deriv('symmetric')
    577        else:
    578            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

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

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

    625    def m_eff(self, variant='log', guess=1.0):
    626        """Returns the effective mass of the correlator as correlator object
    627
    628        Parameters
    629        ----------
    630        variant : str
    631            log : uses the standard effective mass log(C(t) / C(t+1))
    632            cosh, periodic : Use periodicitiy of the correlator by solving C(t) / C(t+1) = cosh(m * (t - T/2)) / cosh(m * (t + 1 - T/2)) for m.
    633            sinh : Use anti-periodicitiy of the correlator by solving C(t) / C(t+1) = sinh(m * (t - T/2)) / sinh(m * (t + 1 - T/2)) for m.
    634            See, e.g., arXiv:1205.5380
    635            arccosh : Uses the explicit form of the symmetrized correlator (not recommended)
    636            logsym: uses the symmetric effective mass log(C(t-1) / C(t+1))/2
    637        guess : float
    638            guess for the root finder, only relevant for the root variant
    639        """
    640        if self.N != 1:
    641            raise Exception('Correlator must be projected before getting m_eff')
    642        if variant == 'log':
    643            newcontent = []
    644            for t in range(self.T - 1):
    645                if ((self.content[t] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0):
    646                    newcontent.append(None)
    647                elif self.content[t][0].value / self.content[t + 1][0].value < 0:
    648                    newcontent.append(None)
    649                else:
    650                    newcontent.append(self.content[t] / self.content[t + 1])
    651            if (all([x is None for x in newcontent])):
    652                raise Exception('m_eff is undefined at all timeslices')
    653
    654            return np.log(Corr(newcontent, padding=[0, 1]))
    655
    656        elif variant == 'logsym':
    657            newcontent = []
    658            for t in range(1, self.T - 1):
    659                if ((self.content[t - 1] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0):
    660                    newcontent.append(None)
    661                elif self.content[t - 1][0].value / self.content[t + 1][0].value < 0:
    662                    newcontent.append(None)
    663                else:
    664                    newcontent.append(self.content[t - 1] / self.content[t + 1])
    665            if (all([x is None for x in newcontent])):
    666                raise Exception('m_eff is undefined at all timeslices')
    667
    668            return np.log(Corr(newcontent, padding=[1, 1])) / 2
    669
    670        elif variant in ['periodic', 'cosh', 'sinh']:
    671            if variant in ['periodic', 'cosh']:
    672                func = anp.cosh
    673            else:
    674                func = anp.sinh
    675
    676            def root_function(x, d):
    677                return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d
    678
    679            newcontent = []
    680            for t in range(self.T - 1):
    681                if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t + 1][0].value == 0):
    682                    newcontent.append(None)
    683                # Fill the two timeslices in the middle of the lattice with their predecessors
    684                elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]:
    685                    newcontent.append(newcontent[-1])
    686                elif self.content[t][0].value / self.content[t + 1][0].value < 0:
    687                    newcontent.append(None)
    688                else:
    689                    newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess)))
    690            if (all([x is None for x in newcontent])):
    691                raise Exception('m_eff is undefined at all timeslices')
    692
    693            return Corr(newcontent, padding=[0, 1])
    694
    695        elif variant == 'arccosh':
    696            newcontent = []
    697            for t in range(1, self.T - 1):
    698                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):
    699                    newcontent.append(None)
    700                else:
    701                    newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t]))
    702            if (all([x is None for x in newcontent])):
    703                raise Exception("m_eff is undefined at all timeslices")
    704            return np.arccosh(Corr(newcontent, padding=[1, 1]))
    705
    706        else:
    707            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 periodicitiy 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-periodicitiy 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

    709    def fit(self, function, fitrange=None, silent=False, **kwargs):
    710        r'''Fits function to the data
    711
    712        Parameters
    713        ----------
    714        function : obj
    715            function to fit to the data. See fits.least_squares for details.
    716        fitrange : list
    717            Two element list containing the timeslices on which the fit is supposed to start and stop.
    718            Caution: This range is inclusive as opposed to standard python indexing.
    719            `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6.
    720            If not specified, self.prange or all timeslices are used.
    721        silent : bool
    722            Decides whether output is printed to the standard output.
    723        '''
    724        if self.N != 1:
    725            raise Exception("Correlator must be projected before fitting")
    726
    727        if fitrange is None:
    728            if self.prange:
    729                fitrange = self.prange
    730            else:
    731                fitrange = [0, self.T - 1]
    732        else:
    733            if not isinstance(fitrange, list):
    734                raise Exception("fitrange has to be a list with two elements")
    735            if len(fitrange) != 2:
    736                raise Exception("fitrange has to have exactly two elements [fit_start, fit_stop]")
    737
    738        xs = np.array([x for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None])
    739        ys = np.array([self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None])
    740        result = least_squares(xs, ys, function, silent=silent, **kwargs)
    741        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

    743    def plateau(self, plateau_range=None, method="fit", auto_gamma=False):
    744        """ Extract a plateau value from a Corr object
    745
    746        Parameters
    747        ----------
    748        plateau_range : list
    749            list with two entries, indicating the first and the last timeslice
    750            of the plateau region.
    751        method : str
    752            method to extract the plateau.
    753                'fit' fits a constant to the plateau region
    754                'avg', 'average' or 'mean' just average over the given timeslices.
    755        auto_gamma : bool
    756            apply gamma_method with default parameters to the Corr. Defaults to None
    757        """
    758        if not plateau_range:
    759            if self.prange:
    760                plateau_range = self.prange
    761            else:
    762                raise Exception("no plateau range provided")
    763        if self.N != 1:
    764            raise Exception("Correlator must be projected before getting a plateau.")
    765        if (all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])):
    766            raise Exception("plateau is undefined at all timeslices in plateaurange.")
    767        if auto_gamma:
    768            self.gamma_method()
    769        if method == "fit":
    770            def const_func(a, t):
    771                return a[0]
    772            return self.fit(const_func, plateau_range)[0]
    773        elif method in ["avg", "average", "mean"]:
    774            returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None])
    775            return returnvalue
    776
    777        else:
    778            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

    780    def set_prange(self, prange):
    781        """Sets the attribute prange of the Corr object."""
    782        if not len(prange) == 2:
    783            raise Exception("prange must be a list or array with two values")
    784        if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))):
    785            raise Exception("Start and end point must be integers")
    786        if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] < prange[1]):
    787            raise Exception("Start and end point must define a range in the interval 0,T")
    788
    789        self.prange = prange
    790        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

    792    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):
    793        """Plots the correlator using the tag of the correlator as label if available.
    794
    795        Parameters
    796        ----------
    797        x_range : list
    798            list of two values, determining the range of the x-axis e.g. [4, 8].
    799        comp : Corr or list of Corr
    800            Correlator or list of correlators which are plotted for comparison.
    801            The tags of these correlators are used as labels if available.
    802        logscale : bool
    803            Sets y-axis to logscale.
    804        plateau : Obs
    805            Plateau value to be visualized in the figure.
    806        fit_res : Fit_result
    807            Fit_result object to be visualized.
    808        fit_key : str
    809            Key for the fit function in Fit_result.fit_function (for combined fits).
    810        ylabel : str
    811            Label for the y-axis.
    812        save : str
    813            path to file in which the figure should be saved.
    814        auto_gamma : bool
    815            Apply the gamma method with standard parameters to all correlators and plateau values before plotting.
    816        hide_sigma : float
    817            Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors.
    818        references : list
    819            List of floating point values that are displayed as horizontal lines for reference.
    820        title : string
    821            Optional title of the figure.
    822        """
    823        if self.N != 1:
    824            raise Exception("Correlator must be projected before plotting")
    825
    826        if auto_gamma:
    827            self.gamma_method()
    828
    829        if x_range is None:
    830            x_range = [0, self.T - 1]
    831
    832        fig = plt.figure()
    833        ax1 = fig.add_subplot(111)
    834
    835        x, y, y_err = self.plottable()
    836        if hide_sigma:
    837            hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
    838        else:
    839            hide_from = None
    840        ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag)
    841        if logscale:
    842            ax1.set_yscale('log')
    843        else:
    844            if y_range is None:
    845                try:
    846                    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)])
    847                    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)])
    848                    ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)])
    849                except Exception:
    850                    pass
    851            else:
    852                ax1.set_ylim(y_range)
    853        if comp:
    854            if isinstance(comp, (Corr, list)):
    855                for corr in comp if isinstance(comp, list) else [comp]:
    856                    if auto_gamma:
    857                        corr.gamma_method()
    858                    x, y, y_err = corr.plottable()
    859                    if hide_sigma:
    860                        hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
    861                    else:
    862                        hide_from = None
    863                    ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor'])
    864            else:
    865                raise Exception("'comp' must be a correlator or a list of correlators.")
    866
    867        if plateau:
    868            if isinstance(plateau, Obs):
    869                if auto_gamma:
    870                    plateau.gamma_method()
    871                ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau))
    872                ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-')
    873            else:
    874                raise Exception("'plateau' must be an Obs")
    875
    876        if references:
    877            if isinstance(references, list):
    878                for ref in references:
    879                    ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--')
    880            else:
    881                raise Exception("'references' must be a list of floating pint values.")
    882
    883        if self.prange:
    884            ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',', color="black", zorder=0)
    885            ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',', color="black", zorder=0)
    886
    887        if fit_res:
    888            x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05)
    889            if isinstance(fit_res.fit_function, dict):
    890                if fit_key:
    891                    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)
    892                else:
    893                    raise ValueError("Please provide a 'fit_key' for visualizing combined fits.")
    894            else:
    895                ax1.plot(x_samples, fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples), ls='-', marker=',', lw=2)
    896
    897        ax1.set_xlabel(r'$x_0 / a$')
    898        if ylabel:
    899            ax1.set_ylabel(ylabel)
    900        ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5])
    901
    902        handles, labels = ax1.get_legend_handles_labels()
    903        if labels:
    904            ax1.legend()
    905
    906        if title:
    907            plt.title(title)
    908
    909        plt.draw()
    910
    911        if save:
    912            if isinstance(save, str):
    913                fig.savefig(save, bbox_inches='tight')
    914            else:
    915                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

    917    def spaghetti_plot(self, logscale=True):
    918        """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations.
    919
    920        Parameters
    921        ----------
    922        logscale : bool
    923            Determines whether the scale of the y-axis is logarithmic or standard.
    924        """
    925        if self.N != 1:
    926            raise Exception("Correlator needs to be projected first.")
    927
    928        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]))
    929        x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None]
    930
    931        for name in mc_names:
    932            data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T
    933
    934            fig = plt.figure()
    935            ax = fig.add_subplot(111)
    936            for dat in data:
    937                ax.plot(x0_vals, dat, ls='-', marker='')
    938
    939            if logscale is True:
    940                ax.set_yscale('log')
    941
    942            ax.set_xlabel(r'$x_0 / a$')
    943            plt.title(name)
    944            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

    946    def dump(self, filename, datatype="json.gz", **kwargs):
    947        """Dumps the Corr into a file of chosen type
    948        Parameters
    949        ----------
    950        filename : str
    951            Name of the file to be saved.
    952        datatype : str
    953            Format of the exported file. Supported formats include
    954            "json.gz" and "pickle"
    955        path : str
    956            specifies a custom path for the file (default '.')
    957        """
    958        if datatype == "json.gz":
    959            from .input.json import dump_to_json
    960            if 'path' in kwargs:
    961                file_name = kwargs.get('path') + '/' + filename
    962            else:
    963                file_name = filename
    964            dump_to_json(self, file_name)
    965        elif datatype == "pickle":
    966            dump_object(self, filename, **kwargs)
    967        else:
    968            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

    970    def print(self, print_range=None):
    971        print(self.__repr__(print_range))

def sqrt(self): View Source

    1137    def sqrt(self):
    1138        return self ** 0.5

def log(self): View Source

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

def exp(self): View Source

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

def sin(self): View Source

    1161    def sin(self):
    1162        return self._apply_func_to_corr(np.sin)

def cos(self): View Source

    1164    def cos(self):
    1165        return self._apply_func_to_corr(np.cos)

def tan(self): View Source

    1167    def tan(self):
    1168        return self._apply_func_to_corr(np.tan)

def sinh(self): View Source

    1170    def sinh(self):
    1171        return self._apply_func_to_corr(np.sinh)

def cosh(self): View Source

    1173    def cosh(self):
    1174        return self._apply_func_to_corr(np.cosh)

def tanh(self): View Source

    1176    def tanh(self):
    1177        return self._apply_func_to_corr(np.tanh)

def arcsin(self): View Source

    1179    def arcsin(self):
    1180        return self._apply_func_to_corr(np.arcsin)

def arccos(self): View Source

    1182    def arccos(self):
    1183        return self._apply_func_to_corr(np.arccos)

def arctan(self): View Source

    1185    def arctan(self):
    1186        return self._apply_func_to_corr(np.arctan)

def arcsinh(self): View Source

    1188    def arcsinh(self):
    1189        return self._apply_func_to_corr(np.arcsinh)

def arccosh(self): View Source

    1191    def arccosh(self):
    1192        return self._apply_func_to_corr(np.arccosh)

def arctanh(self): View Source

    1194    def arctanh(self):
    1195        return self._apply_func_to_corr(np.arctanh)

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

    1230    def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
    1231        r''' Project large correlation matrix to lowest states
    1232
    1233        This method can be used to reduce the size of an (N x N) correlation matrix
    1234        to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise
    1235        is still small.
    1236
    1237        Parameters
    1238        ----------
    1239        Ntrunc: int
    1240            Rank of the target matrix.
    1241        tproj: int
    1242            Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method.
    1243            The default value is 3.
    1244        t0proj: int
    1245            Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly
    1246            discouraged for O(a) improved theories, since the correctness of the procedure
    1247            cannot be granted in this case. The default value is 2.
    1248        basematrix : Corr
    1249            Correlation matrix that is used to determine the eigenvectors of the
    1250            lowest states based on a GEVP. basematrix is taken to be the Corr itself if
    1251            is is not specified.
    1252
    1253        Notes
    1254        -----
    1255        We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving
    1256        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}$
    1257        and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the
    1258        resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via
    1259        $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large
    1260        correlation matrix and to remove some noise that is added by irrelevant operators.
    1261        This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated
    1262        bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$.
    1263        '''
    1264
    1265        if self.N == 1:
    1266            raise Exception('Method cannot be applied to one-dimensional correlators.')
    1267        if basematrix is None:
    1268            basematrix = self
    1269        if Ntrunc >= basematrix.N:
    1270            raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N))
    1271        if basematrix.N != self.N:
    1272            raise Exception('basematrix and targetmatrix have to be of the same size.')
    1273
    1274        evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc]
    1275
    1276        tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object)
    1277        rmat = []
    1278        for t in range(basematrix.T):
    1279            for i in range(Ntrunc):
    1280                for j in range(Ntrunc):
    1281                    tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j]
    1282            rmat.append(np.copy(tmpmat))
    1283
    1284        newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)]
    1285        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)$.