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

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

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

    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)

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

    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)

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

    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

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

    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 np.argmax(np.abs(self.content)) != 0:
    207            warnings.warn("Correlator does not seem to be symmetric around x0=0.", RuntimeWarning)
    208
    209        newcontent = [self.content[0]]
    210        for t in range(1, self.T):
    211            if (self.content[t] is None) or (self.content[self.T - t] is None):
    212                newcontent.append(None)
    213            else:
    214                newcontent.append(0.5 * (self.content[t] + self.content[self.T - t]))
    215        if (all([x is None for x in newcontent])):
    216            raise Exception("Corr could not be symmetrized: No redundant values")
    217        return Corr(newcontent, prange=self.prange)

Symmetrize the correlator around x0=0.

def anti_symmetric(self): View Source

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

Anti-symmetrize the correlator around x0=0.

def is_matrix_symmetric(self): View Source

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

Checks whether a correlator matrices is symmetric on every timeslice.

def matrix_symmetric(self): View Source

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

Symmetrizes the correlator matrices on every timeslice.

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

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

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

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

    406    def roll(self, dt):
    407        """Periodically shift the correlator by dt timeslices
    408
    409        Parameters
    410        ----------
    411        dt : int
    412            number of timeslices
    413        """
    414        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

    416    def reverse(self):
    417        """Reverse the time ordering of the Corr"""
    418        return Corr(self.content[:: -1])

Reverse the time ordering of the Corr

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

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

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

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

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

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

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

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

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

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

    785    def set_prange(self, prange):
    786        """Sets the attribute prange of the Corr object."""
    787        if not len(prange) == 2:
    788            raise Exception("prange must be a list or array with two values")
    789        if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))):
    790            raise Exception("Start and end point must be integers")
    791        if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] < prange[1]):
    792            raise Exception("Start and end point must define a range in the interval 0,T")
    793
    794        self.prange = prange
    795        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, ylabel=None, save=None, auto_gamma=False, hide_sigma=None, references=None, title=None): View Source

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

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.
- 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

    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 [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()

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

    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))

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

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

def sqrt(self): View Source

    1134    def sqrt(self):
    1135        return self ** 0.5

def log(self): View Source

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

def exp(self): View Source

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

def sin(self): View Source

    1156    def sin(self):
    1157        return self._apply_func_to_corr(np.sin)

def cos(self): View Source

    1159    def cos(self):
    1160        return self._apply_func_to_corr(np.cos)

def tan(self): View Source

    1162    def tan(self):
    1163        return self._apply_func_to_corr(np.tan)

def sinh(self): View Source

    1165    def sinh(self):
    1166        return self._apply_func_to_corr(np.sinh)

def cosh(self): View Source

    1168    def cosh(self):
    1169        return self._apply_func_to_corr(np.cosh)

def tanh(self): View Source

    1171    def tanh(self):
    1172        return self._apply_func_to_corr(np.tanh)

def arcsin(self): View Source

    1174    def arcsin(self):
    1175        return self._apply_func_to_corr(np.arcsin)

def arccos(self): View Source

    1177    def arccos(self):
    1178        return self._apply_func_to_corr(np.arccos)

def arctan(self): View Source

    1180    def arctan(self):
    1181        return self._apply_func_to_corr(np.arctan)

def arcsinh(self): View Source

    1183    def arcsinh(self):
    1184        return self._apply_func_to_corr(np.arcsinh)

def arccosh(self): View Source

    1186    def arccosh(self):
    1187        return self._apply_func_to_corr(np.arccosh)

def arctanh(self): View Source

    1189    def arctanh(self):
    1190        return self._apply_func_to_corr(np.arctanh)

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

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