# pyerrors.correlators

View Source

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

class Corr:
View Source

      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 (item is None) 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 (self.content[t] is None 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.T % 2 != 0:
     201            raise Exception("Can not symmetrize odd T")
     202
     203        if np.argmax(np.abs(self.content)) != 0:
     204            warnings.warn("Correlator does not seem to be symmetric around x0=0.", RuntimeWarning)
     205
     206        newcontent = [self.content[0]]
     207        for t in range(1, self.T):
     208            if (self.content[t] is None) or (self.content[self.T - t] is None):
     209                newcontent.append(None)
     210            else:
     211                newcontent.append(0.5 * (self.content[t] + self.content[self.T - t]))
     212        if(all([x is None for x in newcontent])):
     213            raise Exception("Corr could not be symmetrized: No redundant values")
     214        return Corr(newcontent, prange=self.prange)
     215
     216    def anti_symmetric(self):
     217        """Anti-symmetrize the correlator around x0=0."""
     218        if self.T % 2 != 0:
     219            raise Exception("Can not symmetrize odd T")
     220
     221        test = 1 * self
     222        test.gamma_method()
     223        if not all([o.is_zero_within_error(3) for o in test.content[0]]):
     224            warnings.warn("Correlator does not seem to be anti-symmetric around x0=0.", RuntimeWarning)
     225
     226        newcontent = [self.content[0]]
     227        for t in range(1, self.T):
     228            if (self.content[t] is None) or (self.content[self.T - t] is None):
     229                newcontent.append(None)
     230            else:
     231                newcontent.append(0.5 * (self.content[t] - self.content[self.T - t]))
     232        if(all([x is None for x in newcontent])):
     233            raise Exception("Corr could not be symmetrized: No redundant values")
     234        return Corr(newcontent, prange=self.prange)
     235
     236    def matrix_symmetric(self):
     237        """Symmetrizes the correlator matrices on every timeslice."""
     238        if self.N > 1:
     239            transposed = [None if (G is None) else G.T for G in self.content]
     240            return 0.5 * (Corr(transposed) + self)
     241        if self.N == 1:
     242            raise Exception("Trying to symmetrize a correlator matrix, that already has N=1.")
     243
     244    def GEVP(self, t0, ts=None, state=0, sorted_list="Eigenvalue"):
     245        """Solve the general eigenvalue problem on the current correlator
     246
     247        Parameters
     248        ----------
     249        t0 : int
     250            The time t0 for G(t)v= lambda G(t_0)v
     251        ts : int
     252            fixed time G(t_s)v= lambda G(t_0)v  if return_list=False
     253            If return_list=True and sorting=Eigenvector it gives a reference point for the sorting method.
     254        state : int
     255            The state one is interested in, ordered by energy. The lowest state is zero.
     256        sorted_list : string
     257            if this argument is set, a list of vectors (len=self.T) is returned. If it is left as None, only one vector is returned.
     258             "Eigenvalue"  -  The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice.
     259             "Eigenvector" -  Use the method described in arXiv:2004.10472 [hep-lat] to find the set of v(t) belonging to the state.
     260                              The reference state is identified by its eigenvalue at t=ts
     261        """
     262
     263        if self.N == 1:
     264            raise Exception("GEVP methods only works on correlator matrices and not single correlators.")
     265
     266        symmetric_corr = self.matrix_symmetric()
     267        if sorted_list is None:
     268            if (ts is None):
     269                raise Exception("ts is required if sorted_list=None")
     270            if (self.content[t0] is None) or (self.content[ts] is None):
     271                raise Exception("Corr not defined at t0/ts")
     272            G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double")
     273            for i in range(self.N):
     274                for j in range(self.N):
     275                    G0[i, j] = symmetric_corr[t0][i, j].value
     276                    Gt[i, j] = symmetric_corr[ts][i, j].value
     277
     278            sp_vecs = _GEVP_solver(Gt, G0)
     279            sp_vec = sp_vecs[state]
     280            return sp_vec
     281        elif sorted_list in ["Eigenvalue", "Eigenvector"]:
     282            all_vecs = []
     283            for t in range(self.T):
     284                try:
     285                    G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double")
     286                    for i in range(self.N):
     287                        for j in range(self.N):
     288                            G0[i, j] = symmetric_corr[t0][i, j].value
     289                            Gt[i, j] = symmetric_corr[t][i, j].value
     290
     291                    sp_vecs = _GEVP_solver(Gt, G0)
     292                    if sorted_list == "Eigenvalue":
     293                        sp_vec = sp_vecs[state]
     294                        all_vecs.append(sp_vec)
     295                    else:
     296                        all_vecs.append(sp_vecs)
     297                except Exception:
     298                    all_vecs.append(None)
     299            if sorted_list == "Eigenvector":
     300                if (ts is None):
     301                    raise Exception("ts is required for the Eigenvector sorting method.")
     302                all_vecs = _sort_vectors(all_vecs, ts)
     303                all_vecs = [a[state] for a in all_vecs]
     304        else:
     305            raise Exception("Unkown value for 'sorted_list'.")
     306
     307        return all_vecs
     308
     309    def Eigenvalue(self, t0, ts=None, state=0, sorted_list=None):
     310        """Determines the eigenvalue of the GEVP by solving and projecting the correlator
     311
     312        Parameters
     313        ----------
     314        t0 : int
     315            The time t0 for G(t)v= lambda G(t_0)v
     316        ts : int
     317            fixed time G(t_s)v= lambda G(t_0)v  if return_list=False
     318            If return_list=True and sorting=Eigenvector it gives a reference point for the sorting method.
     319        state : int
     320            The state one is interested in ordered by energy. The lowest state is zero.
     321        sorted_list : string
     322            if this argument is set, a list of vectors (len=self.T) is returned. If it is left as None, only one vector is returned.
     323             "Eigenvalue"  -  The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice.
     324             "Eigenvector" -  Use the method described in arXiv:2004.10472 [hep-lat] to find the set of v(t) belonging to the state.
     325                              The reference state is identified by its eigenvalue at t=ts
     326        """
     327        vec = self.GEVP(t0, ts=ts, state=state, sorted_list=sorted_list)
     328        return self.projected(vec)
     329
     330    def Hankel(self, N, periodic=False):
     331        """Constructs an NxN Hankel matrix
     332
     333        C(t) c(t+1) ... c(t+n-1)
     334        C(t+1) c(t+2) ... c(t+n)
     335        .................
     336        C(t+(n-1)) c(t+n) ... c(t+2(n-1))
     337
     338        Parameters
     339        ----------
     340        N : int
     341            Dimension of the Hankel matrix
     342        periodic : bool, optional
     343            determines whether the matrix is extended periodically
     344        """
     345
     346        if self.N != 1:
     347            raise Exception("Multi-operator Prony not implemented!")
     348
     349        array = np.empty([N, N], dtype="object")
     350        new_content = []
     351        for t in range(self.T):
     352            new_content.append(array.copy())
     353
     354        def wrap(i):
     355            while i >= self.T:
     356                i -= self.T
     357            return i
     358
     359        for t in range(self.T):
     360            for i in range(N):
     361                for j in range(N):
     362                    if periodic:
     363                        new_content[t][i, j] = self.content[wrap(t + i + j)][0]
     364                    elif (t + i + j) >= self.T:
     365                        new_content[t] = None
     366                    else:
     367                        new_content[t][i, j] = self.content[t + i + j][0]
     368
     369        return Corr(new_content)
     370
     371    def roll(self, dt):
     372        """Periodically shift the correlator by dt timeslices
     373
     374        Parameters
     375        ----------
     376        dt : int
     377            number of timeslices
     378        """
     379        return Corr(list(np.roll(np.array(self.content, dtype=object), dt)))
     380
     381    def reverse(self):
     382        """Reverse the time ordering of the Corr"""
     383        return Corr(self.content[:: -1])
     384
     385    def thin(self, spacing=2, offset=0):
     386        """Thin out a correlator to suppress correlations
     387
     388        Parameters
     389        ----------
     390        spacing : int
     391            Keep only every 'spacing'th entry of the correlator
     392        offset : int
     393            Offset the equal spacing
     394        """
     395        new_content = []
     396        for t in range(self.T):
     397            if (offset + t) % spacing != 0:
     398                new_content.append(None)
     399            else:
     400                new_content.append(self.content[t])
     401        return Corr(new_content)
     402
     403    def correlate(self, partner):
     404        """Correlate the correlator with another correlator or Obs
     405
     406        Parameters
     407        ----------
     408        partner : Obs or Corr
     409            partner to correlate the correlator with.
     410            Can either be an Obs which is correlated with all entries of the
     411            correlator or a Corr of same length.
     412        """
     413        new_content = []
     414        for x0, t_slice in enumerate(self.content):
     415            if t_slice is None:
     416                new_content.append(None)
     417            else:
     418                if isinstance(partner, Corr):
     419                    if partner.content[x0] is None:
     420                        new_content.append(None)
     421                    else:
     422                        new_content.append(np.array([correlate(o, partner.content[x0][0]) for o in t_slice]))
     423                elif isinstance(partner, Obs):  # Should this include CObs?
     424                    new_content.append(np.array([correlate(o, partner) for o in t_slice]))
     425                else:
     426                    raise Exception("Can only correlate with an Obs or a Corr.")
     427
     428        return Corr(new_content)
     429
     430    def reweight(self, weight, **kwargs):
     431        """Reweight the correlator.
     432
     433        Parameters
     434        ----------
     435        weight : Obs
     436            Reweighting factor. An Observable that has to be defined on a superset of the
     437            configurations in obs[i].idl for all i.
     438        all_configs : bool
     439            if True, the reweighted observables are normalized by the average of
     440            the reweighting factor on all configurations in weight.idl and not
     441            on the configurations in obs[i].idl.
     442        """
     443        new_content = []
     444        for t_slice in self.content:
     445            if t_slice is None:
     446                new_content.append(None)
     447            else:
     448                new_content.append(np.array(reweight(weight, t_slice, **kwargs)))
     449        return Corr(new_content)
     450
     451    def T_symmetry(self, partner, parity=+1):
     452        """Return the time symmetry average of the correlator and its partner
     453
     454        Parameters
     455        ----------
     456        partner : Corr
     457            Time symmetry partner of the Corr
     458        partity : int
     459            Parity quantum number of the correlator, can be +1 or -1
     460        """
     461        if not isinstance(partner, Corr):
     462            raise Exception("T partner has to be a Corr object.")
     463        if parity not in [+1, -1]:
     464            raise Exception("Parity has to be +1 or -1.")
     465        T_partner = parity * partner.reverse()
     466
     467        t_slices = []
     468        test = (self - T_partner)
     469        test.gamma_method()
     470        for x0, t_slice in enumerate(test.content):
     471            if t_slice is not None:
     472                if not t_slice[0].is_zero_within_error(5):
     473                    t_slices.append(x0)
     474        if t_slices:
     475            warnings.warn("T symmetry partners do not agree within 5 sigma on time slices " + str(t_slices) + ".", RuntimeWarning)
     476
     477        return (self + T_partner) / 2
     478
     479    def deriv(self, variant="symmetric"):
     480        """Return the first derivative of the correlator with respect to x0.
     481
     482        Parameters
     483        ----------
     484        variant : str
     485            decides which definition of the finite differences derivative is used.
     486            Available choice: symmetric, forward, backward, improved, default: symmetric
     487        """
     488        if variant == "symmetric":
     489            newcontent = []
     490            for t in range(1, self.T - 1):
     491                if (self.content[t - 1] is None) or (self.content[t + 1] is None):
     492                    newcontent.append(None)
     493                else:
     494                    newcontent.append(0.5 * (self.content[t + 1] - self.content[t - 1]))
     495            if(all([x is None for x in newcontent])):
     496                raise Exception('Derivative is undefined at all timeslices')
     497            return Corr(newcontent, padding=[1, 1])
     498        elif variant == "forward":
     499            newcontent = []
     500            for t in range(self.T - 1):
     501                if (self.content[t] is None) or (self.content[t + 1] is None):
     502                    newcontent.append(None)
     503                else:
     504                    newcontent.append(self.content[t + 1] - self.content[t])
     505            if(all([x is None for x in newcontent])):
     506                raise Exception("Derivative is undefined at all timeslices")
     507            return Corr(newcontent, padding=[0, 1])
     508        elif variant == "backward":
     509            newcontent = []
     510            for t in range(1, self.T):
     511                if (self.content[t - 1] is None) or (self.content[t] is None):
     512                    newcontent.append(None)
     513                else:
     514                    newcontent.append(self.content[t] - self.content[t - 1])
     515            if(all([x is None for x in newcontent])):
     516                raise Exception("Derivative is undefined at all timeslices")
     517            return Corr(newcontent, padding=[1, 0])
     518        elif variant == "improved":
     519            newcontent = []
     520            for t in range(2, self.T - 2):
     521                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):
     522                    newcontent.append(None)
     523                else:
     524                    newcontent.append((1 / 12) * (self.content[t - 2] - 8 * self.content[t - 1] + 8 * self.content[t + 1] - self.content[t + 2]))
     525            if(all([x is None for x in newcontent])):
     526                raise Exception('Derivative is undefined at all timeslices')
     527            return Corr(newcontent, padding=[2, 2])
     528        else:
     529            raise Exception("Unknown variant.")
     530
     531    def second_deriv(self, variant="symmetric"):
     532        """Return the second derivative of the correlator with respect to x0.
     533
     534        Parameters
     535        ----------
     536        variant : str
     537            decides which definition of the finite differences derivative is used.
     538            Available choice: symmetric, improved, default: symmetric
     539        """
     540        if variant == "symmetric":
     541            newcontent = []
     542            for t in range(1, self.T - 1):
     543                if (self.content[t - 1] is None) or (self.content[t + 1] is None):
     544                    newcontent.append(None)
     545                else:
     546                    newcontent.append((self.content[t + 1] - 2 * self.content[t] + self.content[t - 1]))
     547            if(all([x is None for x in newcontent])):
     548                raise Exception("Derivative is undefined at all timeslices")
     549            return Corr(newcontent, padding=[1, 1])
     550        elif variant == "improved":
     551            newcontent = []
     552            for t in range(2, self.T - 2):
     553                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):
     554                    newcontent.append(None)
     555                else:
     556                    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]))
     557            if(all([x is None for x in newcontent])):
     558                raise Exception("Derivative is undefined at all timeslices")
     559            return Corr(newcontent, padding=[2, 2])
     560        else:
     561            raise Exception("Unknown variant.")
     562
     563    def m_eff(self, variant='log', guess=1.0):
     564        """Returns the effective mass of the correlator as correlator object
     565
     566        Parameters
     567        ----------
     568        variant : str
     569            log : uses the standard effective mass log(C(t) / C(t+1))
     570            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.
     571            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.
     572            See, e.g., arXiv:1205.5380
     573            arccosh : Uses the explicit form of the symmetrized correlator (not recommended)
     574        guess : float
     575            guess for the root finder, only relevant for the root variant
     576        """
     577        if self.N != 1:
     578            raise Exception('Correlator must be projected before getting m_eff')
     579        if variant == 'log':
     580            newcontent = []
     581            for t in range(self.T - 1):
     582                if (self.content[t] is None) or (self.content[t + 1] is None):
     583                    newcontent.append(None)
     584                else:
     585                    newcontent.append(self.content[t] / self.content[t + 1])
     586            if(all([x is None for x in newcontent])):
     587                raise Exception('m_eff is undefined at all timeslices')
     588
     589            return np.log(Corr(newcontent, padding=[0, 1]))
     590
     591        elif variant in ['periodic', 'cosh', 'sinh']:
     592            if variant in ['periodic', 'cosh']:
     593                func = anp.cosh
     594            else:
     595                func = anp.sinh
     596
     597            def root_function(x, d):
     598                return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d
     599
     600            newcontent = []
     601            for t in range(self.T - 1):
     602                if (self.content[t] is None) or (self.content[t + 1] is None):
     603                    newcontent.append(None)
     604                # Fill the two timeslices in the middle of the lattice with their predecessors
     605                elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]:
     606                    newcontent.append(newcontent[-1])
     607                else:
     608                    newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess)))
     609            if(all([x is None for x in newcontent])):
     610                raise Exception('m_eff is undefined at all timeslices')
     611
     612            return Corr(newcontent, padding=[0, 1])
     613
     614        elif variant == 'arccosh':
     615            newcontent = []
     616            for t in range(1, self.T - 1):
     617                if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t - 1] is None):
     618                    newcontent.append(None)
     619                else:
     620                    newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t]))
     621            if(all([x is None for x in newcontent])):
     622                raise Exception("m_eff is undefined at all timeslices")
     623            return np.arccosh(Corr(newcontent, padding=[1, 1]))
     624
     625        else:
     626            raise Exception('Unknown variant.')
     627
     628    def fit(self, function, fitrange=None, silent=False, **kwargs):
     629        r'''Fits function to the data
     630
     631        Parameters
     632        ----------
     633        function : obj
     634            function to fit to the data. See fits.least_squares for details.
     635        fitrange : list
     636            Two element list containing the timeslices on which the fit is supposed to start and stop.
     637            Caution: This range is inclusive as opposed to standard python indexing.
     638            `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6.
     639            If not specified, self.prange or all timeslices are used.
     640        silent : bool
     641            Decides whether output is printed to the standard output.
     642        '''
     643        if self.N != 1:
     644            raise Exception("Correlator must be projected before fitting")
     645
     646        if fitrange is None:
     647            if self.prange:
     648                fitrange = self.prange
     649            else:
     650                fitrange = [0, self.T - 1]
     651        else:
     652            if not isinstance(fitrange, list):
     653                raise Exception("fitrange has to be a list with two elements")
     654            if len(fitrange) != 2:
     655                raise Exception("fitrange has to have exactly two elements [fit_start, fit_stop]")
     656
     657        xs = [x for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None]
     658        ys = [self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None]
     659        result = least_squares(xs, ys, function, silent=silent, **kwargs)
     660        return result
     661
     662    def plateau(self, plateau_range=None, method="fit", auto_gamma=False):
     663        """ Extract a plateau value from a Corr object
     664
     665        Parameters
     666        ----------
     667        plateau_range : list
     668            list with two entries, indicating the first and the last timeslice
     669            of the plateau region.
     670        method : str
     671            method to extract the plateau.
     672                'fit' fits a constant to the plateau region
     673                'avg', 'average' or 'mean' just average over the given timeslices.
     674        auto_gamma : bool
     675            apply gamma_method with default parameters to the Corr. Defaults to None
     676        """
     677        if not plateau_range:
     678            if self.prange:
     679                plateau_range = self.prange
     680            else:
     681                raise Exception("no plateau range provided")
     682        if self.N != 1:
     683            raise Exception("Correlator must be projected before getting a plateau.")
     684        if(all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])):
     685            raise Exception("plateau is undefined at all timeslices in plateaurange.")
     686        if auto_gamma:
     687            self.gamma_method()
     688        if method == "fit":
     689            def const_func(a, t):
     690                return a[0]
     691            return self.fit(const_func, plateau_range)[0]
     692        elif method in ["avg", "average", "mean"]:
     693            returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None])
     694            return returnvalue
     695
     696        else:
     697            raise Exception("Unsupported plateau method: " + method)
     698
     699    def set_prange(self, prange):
     700        """Sets the attribute prange of the Corr object."""
     701        if not len(prange) == 2:
     702            raise Exception("prange must be a list or array with two values")
     703        if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))):
     704            raise Exception("Start and end point must be integers")
     705        if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] < prange[1]):
     706            raise Exception("Start and end point must define a range in the interval 0,T")
     707
     708        self.prange = prange
     709        return
     710
     711    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):
     712        """Plots the correlator using the tag of the correlator as label if available.
     713
     714        Parameters
     715        ----------
     716        x_range : list
     717            list of two values, determining the range of the x-axis e.g. [4, 8]
     718        comp : Corr or list of Corr
     719            Correlator or list of correlators which are plotted for comparison.
     720            The tags of these correlators are used as labels if available.
     721        logscale : bool
     722            Sets y-axis to logscale
     723        plateau : Obs
     724            Plateau value to be visualized in the figure
     725        fit_res : Fit_result
     726            Fit_result object to be visualized
     727        ylabel : str
     728            Label for the y-axis
     729        save : str
     730            path to file in which the figure should be saved
     731        auto_gamma : bool
     732            Apply the gamma method with standard parameters to all correlators and plateau values before plotting.
     733        hide_sigma : float
     734            Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors.
     735        references : list
     736            List of floating point values that are displayed as horizontal lines for reference.
     737        """
     738        if self.N != 1:
     739            raise Exception("Correlator must be projected before plotting")
     740
     741        if auto_gamma:
     742            self.gamma_method()
     743
     744        if x_range is None:
     745            x_range = [0, self.T - 1]
     746
     747        fig = plt.figure()
     748        ax1 = fig.add_subplot(111)
     749
     750        x, y, y_err = self.plottable()
     751        if hide_sigma:
     752            hide_from = np.argmax((hide_sigma * np.array(y_err)) > np.abs(y)) - 1
     753        else:
     754            hide_from = None
     755        ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag)
     756        if logscale:
     757            ax1.set_yscale('log')
     758        else:
     759            if y_range is None:
     760                try:
     761                    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)])
     762                    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)])
     763                    ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)])
     764                except Exception:
     765                    pass
     766            else:
     767                ax1.set_ylim(y_range)
     768        if comp:
     769            if isinstance(comp, (Corr, list)):
     770                for corr in comp if isinstance(comp, list) else [comp]:
     771                    if auto_gamma:
     772                        corr.gamma_method()
     773                    x, y, y_err = corr.plottable()
     774                    if hide_sigma:
     775                        hide_from = np.argmax((hide_sigma * np.array(y_err)) > np.abs(y)) - 1
     776                    else:
     777                        hide_from = None
     778                    plt.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor'])
     779            else:
     780                raise Exception("'comp' must be a correlator or a list of correlators.")
     781
     782        if plateau:
     783            if isinstance(plateau, Obs):
     784                if auto_gamma:
     785                    plateau.gamma_method()
     786                ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau))
     787                ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-')
     788            else:
     789                raise Exception("'plateau' must be an Obs")
     790
     791        if references:
     792            if isinstance(references, list):
     793                for ref in references:
     794                    ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--')
     795            else:
     796                raise Exception("'references' must be a list of floating pint values.")
     797
     798        if self.prange:
     799            ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',')
     800            ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',')
     801
     802        if fit_res:
     803            x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05)
     804            ax1.plot(x_samples,
     805                     fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples),
     806                     ls='-', marker=',', lw=2)
     807
     808        ax1.set_xlabel(r'$x_0 / a$')
     809        if ylabel:
     810            ax1.set_ylabel(ylabel)
     811        ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5])
     812
     813        handles, labels = ax1.get_legend_handles_labels()
     814        if labels:
     815            ax1.legend()
     816        plt.draw()
     817
     818        if save:
     819            if isinstance(save, str):
     820                fig.savefig(save)
     821            else:
     822                raise Exception("'save' has to be a string.")
     823
     824    def spaghetti_plot(self, logscale=True):
     825        """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations.
     826
     827        Parameters
     828        ----------
     829        logscale : bool
     830            Determines whether the scale of the y-axis is logarithmic or standard.
     831        """
     832        if self.N != 1:
     833            raise Exception("Correlator needs to be projected first.")
     834
     835        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]))
     836        x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None]
     837
     838        for name in mc_names:
     839            data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T
     840
     841            fig = plt.figure()
     842            ax = fig.add_subplot(111)
     843            for dat in data:
     844                ax.plot(x0_vals, dat, ls='-', marker='')
     845
     846            if logscale is True:
     847                ax.set_yscale('log')
     848
     849            ax.set_xlabel(r'$x_0 / a$')
     850            plt.title(name)
     851            plt.draw()
     852
     853    def dump(self, filename, datatype="json.gz", **kwargs):
     854        """Dumps the Corr into a file of chosen type
     855        Parameters
     856        ----------
     857        filename : str
     858            Name of the file to be saved.
     859        datatype : str
     860            Format of the exported file. Supported formats include
     861            "json.gz" and "pickle"
     862        path : str
     863            specifies a custom path for the file (default '.')
     864        """
     865        if datatype == "json.gz":
     866            from .input.json import dump_to_json
     867            if 'path' in kwargs:
     868                file_name = kwargs.get('path') + '/' + filename
     869            else:
     870                file_name = filename
     871            dump_to_json(self, file_name)
     872        elif datatype == "pickle":
     873            dump_object(self, filename, **kwargs)
     874        else:
     875            raise Exception("Unknown datatype " + str(datatype))
     876
     877    def print(self, range=[0, None]):
     878        print(self.__repr__(range))
     879
     880    def __repr__(self, range=[0, None]):
     881        content_string = ""
     882
     883        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
     884
     885        if self.tag is not None:
     886            content_string += "Description: " + self.tag + "\n"
     887        if self.N != 1:
     888            return content_string
     889
     890        if range[1]:
     891            range[1] += 1
     892        content_string += 'x0/a\tCorr(x0/a)\n------------------\n'
     893        for i, sub_corr in enumerate(self.content[range[0]:range[1]]):
     894            if sub_corr is None:
     895                content_string += str(i + range[0]) + '\n'
     896            else:
     897                content_string += str(i + range[0])
     898                for element in sub_corr:
     899                    content_string += '\t' + ' ' * int(element >= 0) + str(element)
     900                content_string += '\n'
     901        return content_string
     902
     903    def __str__(self):
     904        return self.__repr__()
     905
     906    # We define the basic operations, that can be performed with correlators.
     907    # While */+- get defined here, they only work for Corr*Obs and not Obs*Corr.
     908    # This is because Obs*Corr checks Obs.__mul__ first and does not catch an exception.
     909    # One could try and tell Obs to check if the y in __mul__ is a Corr and
     910
     911    def __add__(self, y):
     912        if isinstance(y, Corr):
     913            if ((self.N != y.N) or (self.T != y.T)):
     914                raise Exception("Addition of Corrs with different shape")
     915            newcontent = []
     916            for t in range(self.T):
     917                if (self.content[t] is None) or (y.content[t] is None):
     918                    newcontent.append(None)
     919                else:
     920                    newcontent.append(self.content[t] + y.content[t])
     921            return Corr(newcontent)
     922
     923        elif isinstance(y, (Obs, int, float, CObs)):
     924            newcontent = []
     925            for t in range(self.T):
     926                if (self.content[t] is None):
     927                    newcontent.append(None)
     928                else:
     929                    newcontent.append(self.content[t] + y)
     930            return Corr(newcontent, prange=self.prange)
     931        elif isinstance(y, np.ndarray):
     932            if y.shape == (self.T,):
     933                return Corr(list((np.array(self.content).T + y).T))
     934            else:
     935                raise ValueError("operands could not be broadcast together")
     936        else:
     937            raise TypeError("Corr + wrong type")
     938
     939    def __mul__(self, y):
     940        if isinstance(y, Corr):
     941            if not((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T):
     942                raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T")
     943            newcontent = []
     944            for t in range(self.T):
     945                if (self.content[t] is None) or (y.content[t] is None):
     946                    newcontent.append(None)
     947                else:
     948                    newcontent.append(self.content[t] * y.content[t])
     949            return Corr(newcontent)
     950
     951        elif isinstance(y, (Obs, int, float, CObs)):
     952            newcontent = []
     953            for t in range(self.T):
     954                if (self.content[t] is None):
     955                    newcontent.append(None)
     956                else:
     957                    newcontent.append(self.content[t] * y)
     958            return Corr(newcontent, prange=self.prange)
     959        elif isinstance(y, np.ndarray):
     960            if y.shape == (self.T,):
     961                return Corr(list((np.array(self.content).T * y).T))
     962            else:
     963                raise ValueError("operands could not be broadcast together")
     964        else:
     965            raise TypeError("Corr * wrong type")
     966
     967    def __truediv__(self, y):
     968        if isinstance(y, Corr):
     969            if not((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T):
     970                raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T")
     971            newcontent = []
     972            for t in range(self.T):
     973                if (self.content[t] is None) or (y.content[t] is None):
     974                    newcontent.append(None)
     975                else:
     976                    newcontent.append(self.content[t] / y.content[t])
     977            for t in range(self.T):
     978                if newcontent[t] is None:
     979                    continue
     980                if np.isnan(np.sum(newcontent[t]).value):
     981                    newcontent[t] = None
     982
     983            if all([item is None for item in newcontent]):
     984                raise Exception("Division returns completely undefined correlator")
     985            return Corr(newcontent)
     986
     987        elif isinstance(y, (Obs, CObs)):
     988            if isinstance(y, Obs):
     989                if y.value == 0:
     990                    raise Exception('Division by zero will return undefined correlator')
     991            if isinstance(y, CObs):
     992                if y.is_zero():
     993                    raise Exception('Division by zero will return undefined correlator')
     994
     995            newcontent = []
     996            for t in range(self.T):
     997                if (self.content[t] is None):
     998                    newcontent.append(None)
     999                else:
    1000                    newcontent.append(self.content[t] / y)
    1001            return Corr(newcontent, prange=self.prange)
    1002
    1003        elif isinstance(y, (int, float)):
    1004            if y == 0:
    1005                raise Exception('Division by zero will return undefined correlator')
    1006            newcontent = []
    1007            for t in range(self.T):
    1008                if (self.content[t] is None):
    1009                    newcontent.append(None)
    1010                else:
    1011                    newcontent.append(self.content[t] / y)
    1012            return Corr(newcontent, prange=self.prange)
    1013        elif isinstance(y, np.ndarray):
    1014            if y.shape == (self.T,):
    1015                return Corr(list((np.array(self.content).T / y).T))
    1016            else:
    1017                raise ValueError("operands could not be broadcast together")
    1018        else:
    1019            raise TypeError('Corr / wrong type')
    1020
    1021    def __neg__(self):
    1022        newcontent = [None if (item is None) else -1. * item for item in self.content]
    1023        return Corr(newcontent, prange=self.prange)
    1024
    1025    def __sub__(self, y):
    1026        return self + (-y)
    1027
    1028    def __pow__(self, y):
    1029        if isinstance(y, (Obs, int, float, CObs)):
    1030            newcontent = [None if (item is None) else item**y for item in self.content]
    1031            return Corr(newcontent, prange=self.prange)
    1032        else:
    1033            raise TypeError('Type of exponent not supported')
    1034
    1035    def __abs__(self):
    1036        newcontent = [None if (item is None) else np.abs(item) for item in self.content]
    1037        return Corr(newcontent, prange=self.prange)
    1038
    1039    # The numpy functions:
    1040    def sqrt(self):
    1041        return self**0.5
    1042
    1043    def log(self):
    1044        newcontent = [None if (item is None) else np.log(item) for item in self.content]
    1045        return Corr(newcontent, prange=self.prange)
    1046
    1047    def exp(self):
    1048        newcontent = [None if (item is None) else np.exp(item) for item in self.content]
    1049        return Corr(newcontent, prange=self.prange)
    1050
    1051    def _apply_func_to_corr(self, func):
    1052        newcontent = [None if (item is None) else func(item) for item in self.content]
    1053        for t in range(self.T):
    1054            if newcontent[t] is None:
    1055                continue
    1056            if np.isnan(np.sum(newcontent[t]).value):
    1057                newcontent[t] = None
    1058        if all([item is None for item in newcontent]):
    1059            raise Exception('Operation returns undefined correlator')
    1060        return Corr(newcontent)
    1061
    1062    def sin(self):
    1063        return self._apply_func_to_corr(np.sin)
    1064
    1065    def cos(self):
    1066        return self._apply_func_to_corr(np.cos)
    1067
    1068    def tan(self):
    1069        return self._apply_func_to_corr(np.tan)
    1070
    1071    def sinh(self):
    1072        return self._apply_func_to_corr(np.sinh)
    1073
    1074    def cosh(self):
    1075        return self._apply_func_to_corr(np.cosh)
    1076
    1077    def tanh(self):
    1078        return self._apply_func_to_corr(np.tanh)
    1079
    1080    def arcsin(self):
    1081        return self._apply_func_to_corr(np.arcsin)
    1082
    1083    def arccos(self):
    1084        return self._apply_func_to_corr(np.arccos)
    1085
    1086    def arctan(self):
    1087        return self._apply_func_to_corr(np.arctan)
    1088
    1089    def arcsinh(self):
    1090        return self._apply_func_to_corr(np.arcsinh)
    1091
    1092    def arccosh(self):
    1093        return self._apply_func_to_corr(np.arccosh)
    1094
    1095    def arctanh(self):
    1096        return self._apply_func_to_corr(np.arctanh)
    1097
    1098    # Right hand side operations (require tweak in main module to work)
    1099    def __radd__(self, y):
    1100        return self + y
    1101
    1102    def __rsub__(self, y):
    1103        return -self + y
    1104
    1105    def __rmul__(self, y):
    1106        return self * y
    1107
    1108    def __rtruediv__(self, y):
    1109        return (self / y) ** (-1)
    1110
    1111    @property
    1112    def real(self):
    1113        def return_real(obs_OR_cobs):
    1114            if isinstance(obs_OR_cobs, CObs):
    1115                return obs_OR_cobs.real
    1116            else:
    1117                return obs_OR_cobs
    1118
    1119        return self._apply_func_to_corr(return_real)
    1120
    1121    @property
    1122    def imag(self):
    1123        def return_imag(obs_OR_cobs):
    1124            if isinstance(obs_OR_cobs, CObs):
    1125                return obs_OR_cobs.imag
    1126            else:
    1127                return obs_OR_cobs * 0  # So it stays the right type
    1128
    1129        return self._apply_func_to_corr(return_imag)
    1130
    1131    def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
    1132        r''' Project large correlation matrix to lowest states
    1133
    1134        This method can be used to reduce the size of an (N x N) correlation matrix
    1135        to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise
    1136        is still small.
    1137
    1138        Parameters
    1139        ----------
    1140        Ntrunc: int
    1141            Rank of the target matrix.
    1142        tproj: int
    1143            Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method.
    1144            The default value is 3.
    1145        t0proj: int
    1146            Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly
    1147            discouraged for O(a) improved theories, since the correctness of the procedure
    1148            cannot be granted in this case. The default value is 2.
    1149        basematrix : Corr
    1150            Correlation matrix that is used to determine the eigenvectors of the
    1151            lowest states based on a GEVP. basematrix is taken to be the Corr itself if
    1152            is is not specified.
    1153
    1154        Notes
    1155        -----
    1156        We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving
    1157        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}$
    1158        and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the
    1159        resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via
    1160        $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large
    1161        correlation matrix and to remove some noise that is added by irrelevant operators.
    1162        This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated
    1163        bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$.
    1164        '''
    1165
    1166        if self.N == 1:
    1167            raise Exception('Method cannot be applied to one-dimensional correlators.')
    1168        if basematrix is None:
    1169            basematrix = self
    1170        if Ntrunc >= basematrix.N:
    1171            raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N))
    1172        if basematrix.N != self.N:
    1173            raise Exception('basematrix and targetmatrix have to be of the same size.')
    1174
    1175        evecs = []
    1176        for i in range(Ntrunc):
    1177            evecs.append(basematrix.GEVP(t0proj, tproj, state=i, sorted_list=None))
    1178
    1179        tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object)
    1180        rmat = []
    1181        for t in range(basematrix.T):
    1182            for i in range(Ntrunc):
    1183                for j in range(Ntrunc):
    1184                    tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j]
    1185            rmat.append(np.copy(tmpmat))
    1186
    1187        newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)]
    1188        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

    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

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.

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

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

Apply the gamma method to the content of the Corr.

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

    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 (item is None) 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 (self.content[t] is None 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)

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

    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)

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

    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

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

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

Symmetrize the correlator around x0=0.

def anti_symmetric(self):
View Source

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

Anti-symmetrize the correlator around x0=0.

def matrix_symmetric(self):
View Source

    236    def matrix_symmetric(self):
    237        """Symmetrizes the correlator matrices on every timeslice."""
    238        if self.N > 1:
    239            transposed = [None if (G is None) else G.T for G in self.content]
    240            return 0.5 * (Corr(transposed) + self)
    241        if self.N == 1:
    242            raise Exception("Trying to symmetrize a correlator matrix, that already has N=1.")

Symmetrizes the correlator matrices on every timeslice.

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

    244    def GEVP(self, t0, ts=None, state=0, sorted_list="Eigenvalue"):
    245        """Solve the general eigenvalue problem on the current correlator
    246
    247        Parameters
    248        ----------
    249        t0 : int
    250            The time t0 for G(t)v= lambda G(t_0)v
    251        ts : int
    252            fixed time G(t_s)v= lambda G(t_0)v  if return_list=False
    253            If return_list=True and sorting=Eigenvector it gives a reference point for the sorting method.
    254        state : int
    255            The state one is interested in, ordered by energy. The lowest state is zero.
    256        sorted_list : string
    257            if this argument is set, a list of vectors (len=self.T) is returned. If it is left as None, only one vector is returned.
    258             "Eigenvalue"  -  The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice.
    259             "Eigenvector" -  Use the method described in arXiv:2004.10472 [hep-lat] to find the set of v(t) belonging to the state.
    260                              The reference state is identified by its eigenvalue at t=ts
    261        """
    262
    263        if self.N == 1:
    264            raise Exception("GEVP methods only works on correlator matrices and not single correlators.")
    265
    266        symmetric_corr = self.matrix_symmetric()
    267        if sorted_list is None:
    268            if (ts is None):
    269                raise Exception("ts is required if sorted_list=None")
    270            if (self.content[t0] is None) or (self.content[ts] is None):
    271                raise Exception("Corr not defined at t0/ts")
    272            G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double")
    273            for i in range(self.N):
    274                for j in range(self.N):
    275                    G0[i, j] = symmetric_corr[t0][i, j].value
    276                    Gt[i, j] = symmetric_corr[ts][i, j].value
    277
    278            sp_vecs = _GEVP_solver(Gt, G0)
    279            sp_vec = sp_vecs[state]
    280            return sp_vec
    281        elif sorted_list in ["Eigenvalue", "Eigenvector"]:
    282            all_vecs = []
    283            for t in range(self.T):
    284                try:
    285                    G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double")
    286                    for i in range(self.N):
    287                        for j in range(self.N):
    288                            G0[i, j] = symmetric_corr[t0][i, j].value
    289                            Gt[i, j] = symmetric_corr[t][i, j].value
    290
    291                    sp_vecs = _GEVP_solver(Gt, G0)
    292                    if sorted_list == "Eigenvalue":
    293                        sp_vec = sp_vecs[state]
    294                        all_vecs.append(sp_vec)
    295                    else:
    296                        all_vecs.append(sp_vecs)
    297                except Exception:
    298                    all_vecs.append(None)
    299            if sorted_list == "Eigenvector":
    300                if (ts is None):
    301                    raise Exception("ts is required for the Eigenvector sorting method.")
    302                all_vecs = _sort_vectors(all_vecs, ts)
    303                all_vecs = [a[state] for a in all_vecs]
    304        else:
    305            raise Exception("Unkown value for 'sorted_list'.")
    306
    307        return all_vecs

Solve the general eigenvalue problem on the current correlator

###### Parameters

- t0 (int): The time t0 for G(t)v= lambda G(t_0)v
- ts (int): fixed time G(t_s)v= lambda G(t_0)v if return_list=False If return_list=True and sorting=Eigenvector it gives a reference point for the sorting method.
- state (int): The state one is interested in, ordered by energy. The lowest state is zero.
- sorted_list (string): if this argument is set, a list of vectors (len=self.T) is returned. If it is left as 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 [hep-lat] to find the set of v(t) belonging to the state. The reference state is identified by its eigenvalue at t=ts

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

    309    def Eigenvalue(self, t0, ts=None, state=0, sorted_list=None):
    310        """Determines the eigenvalue of the GEVP by solving and projecting the correlator
    311
    312        Parameters
    313        ----------
    314        t0 : int
    315            The time t0 for G(t)v= lambda G(t_0)v
    316        ts : int
    317            fixed time G(t_s)v= lambda G(t_0)v  if return_list=False
    318            If return_list=True and sorting=Eigenvector it gives a reference point for the sorting method.
    319        state : int
    320            The state one is interested in ordered by energy. The lowest state is zero.
    321        sorted_list : string
    322            if this argument is set, a list of vectors (len=self.T) is returned. If it is left as None, only one vector is returned.
    323             "Eigenvalue"  -  The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice.
    324             "Eigenvector" -  Use the method described in arXiv:2004.10472 [hep-lat] to find the set of v(t) belonging to the state.
    325                              The reference state is identified by its eigenvalue at t=ts
    326        """
    327        vec = self.GEVP(t0, ts=ts, state=state, sorted_list=sorted_list)
    328        return self.projected(vec)

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

###### Parameters

- t0 (int): The time t0 for G(t)v= lambda G(t_0)v
- ts (int): fixed time G(t_s)v= lambda G(t_0)v if return_list=False If return_list=True and sorting=Eigenvector it gives a reference point for the sorting method.
- state (int): The state one is interested in ordered by energy. The lowest state is zero.
- sorted_list (string): if this argument is set, a list of vectors (len=self.T) is returned. If it is left as 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 [hep-lat] to find the set of v(t) belonging to the state. The reference state is identified by its eigenvalue at t=ts

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

    330    def Hankel(self, N, periodic=False):
    331        """Constructs an NxN Hankel matrix
    332
    333        C(t) c(t+1) ... c(t+n-1)
    334        C(t+1) c(t+2) ... c(t+n)
    335        .................
    336        C(t+(n-1)) c(t+n) ... c(t+2(n-1))
    337
    338        Parameters
    339        ----------
    340        N : int
    341            Dimension of the Hankel matrix
    342        periodic : bool, optional
    343            determines whether the matrix is extended periodically
    344        """
    345
    346        if self.N != 1:
    347            raise Exception("Multi-operator Prony not implemented!")
    348
    349        array = np.empty([N, N], dtype="object")
    350        new_content = []
    351        for t in range(self.T):
    352            new_content.append(array.copy())
    353
    354        def wrap(i):
    355            while i >= self.T:
    356                i -= self.T
    357            return i
    358
    359        for t in range(self.T):
    360            for i in range(N):
    361                for j in range(N):
    362                    if periodic:
    363                        new_content[t][i, j] = self.content[wrap(t + i + j)][0]
    364                    elif (t + i + j) >= self.T:
    365                        new_content[t] = None
    366                    else:
    367                        new_content[t][i, j] = self.content[t + i + j][0]
    368
    369        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

    371    def roll(self, dt):
    372        """Periodically shift the correlator by dt timeslices
    373
    374        Parameters
    375        ----------
    376        dt : int
    377            number of timeslices
    378        """
    379        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

    381    def reverse(self):
    382        """Reverse the time ordering of the Corr"""
    383        return Corr(self.content[:: -1])

Reverse the time ordering of the Corr

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

    385    def thin(self, spacing=2, offset=0):
    386        """Thin out a correlator to suppress correlations
    387
    388        Parameters
    389        ----------
    390        spacing : int
    391            Keep only every 'spacing'th entry of the correlator
    392        offset : int
    393            Offset the equal spacing
    394        """
    395        new_content = []
    396        for t in range(self.T):
    397            if (offset + t) % spacing != 0:
    398                new_content.append(None)
    399            else:
    400                new_content.append(self.content[t])
    401        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

    403    def correlate(self, partner):
    404        """Correlate the correlator with another correlator or Obs
    405
    406        Parameters
    407        ----------
    408        partner : Obs or Corr
    409            partner to correlate the correlator with.
    410            Can either be an Obs which is correlated with all entries of the
    411            correlator or a Corr of same length.
    412        """
    413        new_content = []
    414        for x0, t_slice in enumerate(self.content):
    415            if t_slice is None:
    416                new_content.append(None)
    417            else:
    418                if isinstance(partner, Corr):
    419                    if partner.content[x0] is None:
    420                        new_content.append(None)
    421                    else:
    422                        new_content.append(np.array([correlate(o, partner.content[x0][0]) for o in t_slice]))
    423                elif isinstance(partner, Obs):  # Should this include CObs?
    424                    new_content.append(np.array([correlate(o, partner) for o in t_slice]))
    425                else:
    426                    raise Exception("Can only correlate with an Obs or a Corr.")
    427
    428        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

    430    def reweight(self, weight, **kwargs):
    431        """Reweight the correlator.
    432
    433        Parameters
    434        ----------
    435        weight : Obs
    436            Reweighting factor. An Observable that has to be defined on a superset of the
    437            configurations in obs[i].idl for all i.
    438        all_configs : bool
    439            if True, the reweighted observables are normalized by the average of
    440            the reweighting factor on all configurations in weight.idl and not
    441            on the configurations in obs[i].idl.
    442        """
    443        new_content = []
    444        for t_slice in self.content:
    445            if t_slice is None:
    446                new_content.append(None)
    447            else:
    448                new_content.append(np.array(reweight(weight, t_slice, **kwargs)))
    449        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

    451    def T_symmetry(self, partner, parity=+1):
    452        """Return the time symmetry average of the correlator and its partner
    453
    454        Parameters
    455        ----------
    456        partner : Corr
    457            Time symmetry partner of the Corr
    458        partity : int
    459            Parity quantum number of the correlator, can be +1 or -1
    460        """
    461        if not isinstance(partner, Corr):
    462            raise Exception("T partner has to be a Corr object.")
    463        if parity not in [+1, -1]:
    464            raise Exception("Parity has to be +1 or -1.")
    465        T_partner = parity * partner.reverse()
    466
    467        t_slices = []
    468        test = (self - T_partner)
    469        test.gamma_method()
    470        for x0, t_slice in enumerate(test.content):
    471            if t_slice is not None:
    472                if not t_slice[0].is_zero_within_error(5):
    473                    t_slices.append(x0)
    474        if t_slices:
    475            warnings.warn("T symmetry partners do not agree within 5 sigma on time slices " + str(t_slices) + ".", RuntimeWarning)
    476
    477        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

    479    def deriv(self, variant="symmetric"):
    480        """Return the first derivative of the correlator with respect to x0.
    481
    482        Parameters
    483        ----------
    484        variant : str
    485            decides which definition of the finite differences derivative is used.
    486            Available choice: symmetric, forward, backward, improved, default: symmetric
    487        """
    488        if variant == "symmetric":
    489            newcontent = []
    490            for t in range(1, self.T - 1):
    491                if (self.content[t - 1] is None) or (self.content[t + 1] is None):
    492                    newcontent.append(None)
    493                else:
    494                    newcontent.append(0.5 * (self.content[t + 1] - self.content[t - 1]))
    495            if(all([x is None for x in newcontent])):
    496                raise Exception('Derivative is undefined at all timeslices')
    497            return Corr(newcontent, padding=[1, 1])
    498        elif variant == "forward":
    499            newcontent = []
    500            for t in range(self.T - 1):
    501                if (self.content[t] is None) or (self.content[t + 1] is None):
    502                    newcontent.append(None)
    503                else:
    504                    newcontent.append(self.content[t + 1] - self.content[t])
    505            if(all([x is None for x in newcontent])):
    506                raise Exception("Derivative is undefined at all timeslices")
    507            return Corr(newcontent, padding=[0, 1])
    508        elif variant == "backward":
    509            newcontent = []
    510            for t in range(1, self.T):
    511                if (self.content[t - 1] is None) or (self.content[t] is None):
    512                    newcontent.append(None)
    513                else:
    514                    newcontent.append(self.content[t] - self.content[t - 1])
    515            if(all([x is None for x in newcontent])):
    516                raise Exception("Derivative is undefined at all timeslices")
    517            return Corr(newcontent, padding=[1, 0])
    518        elif variant == "improved":
    519            newcontent = []
    520            for t in range(2, self.T - 2):
    521                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):
    522                    newcontent.append(None)
    523                else:
    524                    newcontent.append((1 / 12) * (self.content[t - 2] - 8 * self.content[t - 1] + 8 * self.content[t + 1] - self.content[t + 2]))
    525            if(all([x is None for x in newcontent])):
    526                raise Exception('Derivative is undefined at all timeslices')
    527            return Corr(newcontent, padding=[2, 2])
    528        else:
    529            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, default: symmetric

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

    531    def second_deriv(self, variant="symmetric"):
    532        """Return the second derivative of the correlator with respect to x0.
    533
    534        Parameters
    535        ----------
    536        variant : str
    537            decides which definition of the finite differences derivative is used.
    538            Available choice: symmetric, improved, default: symmetric
    539        """
    540        if variant == "symmetric":
    541            newcontent = []
    542            for t in range(1, self.T - 1):
    543                if (self.content[t - 1] is None) or (self.content[t + 1] is None):
    544                    newcontent.append(None)
    545                else:
    546                    newcontent.append((self.content[t + 1] - 2 * self.content[t] + self.content[t - 1]))
    547            if(all([x is None for x in newcontent])):
    548                raise Exception("Derivative is undefined at all timeslices")
    549            return Corr(newcontent, padding=[1, 1])
    550        elif variant == "improved":
    551            newcontent = []
    552            for t in range(2, self.T - 2):
    553                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):
    554                    newcontent.append(None)
    555                else:
    556                    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]))
    557            if(all([x is None for x in newcontent])):
    558                raise Exception("Derivative is undefined at all timeslices")
    559            return Corr(newcontent, padding=[2, 2])
    560        else:
    561            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, default: symmetric

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

    563    def m_eff(self, variant='log', guess=1.0):
    564        """Returns the effective mass of the correlator as correlator object
    565
    566        Parameters
    567        ----------
    568        variant : str
    569            log : uses the standard effective mass log(C(t) / C(t+1))
    570            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.
    571            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.
    572            See, e.g., arXiv:1205.5380
    573            arccosh : Uses the explicit form of the symmetrized correlator (not recommended)
    574        guess : float
    575            guess for the root finder, only relevant for the root variant
    576        """
    577        if self.N != 1:
    578            raise Exception('Correlator must be projected before getting m_eff')
    579        if variant == 'log':
    580            newcontent = []
    581            for t in range(self.T - 1):
    582                if (self.content[t] is None) or (self.content[t + 1] is None):
    583                    newcontent.append(None)
    584                else:
    585                    newcontent.append(self.content[t] / self.content[t + 1])
    586            if(all([x is None for x in newcontent])):
    587                raise Exception('m_eff is undefined at all timeslices')
    588
    589            return np.log(Corr(newcontent, padding=[0, 1]))
    590
    591        elif variant in ['periodic', 'cosh', 'sinh']:
    592            if variant in ['periodic', 'cosh']:
    593                func = anp.cosh
    594            else:
    595                func = anp.sinh
    596
    597            def root_function(x, d):
    598                return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d
    599
    600            newcontent = []
    601            for t in range(self.T - 1):
    602                if (self.content[t] is None) or (self.content[t + 1] is None):
    603                    newcontent.append(None)
    604                # Fill the two timeslices in the middle of the lattice with their predecessors
    605                elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]:
    606                    newcontent.append(newcontent[-1])
    607                else:
    608                    newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess)))
    609            if(all([x is None for x in newcontent])):
    610                raise Exception('m_eff is undefined at all timeslices')
    611
    612            return Corr(newcontent, padding=[0, 1])
    613
    614        elif variant == 'arccosh':
    615            newcontent = []
    616            for t in range(1, self.T - 1):
    617                if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t - 1] is None):
    618                    newcontent.append(None)
    619                else:
    620                    newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t]))
    621            if(all([x is None for x in newcontent])):
    622                raise Exception("m_eff is undefined at all timeslices")
    623            return np.arccosh(Corr(newcontent, padding=[1, 1]))
    624
    625        else:
    626            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)
- guess (float): guess for the root finder, only relevant for the root variant

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

    628    def fit(self, function, fitrange=None, silent=False, **kwargs):
    629        r'''Fits function to the data
    630
    631        Parameters
    632        ----------
    633        function : obj
    634            function to fit to the data. See fits.least_squares for details.
    635        fitrange : list
    636            Two element list containing the timeslices on which the fit is supposed to start and stop.
    637            Caution: This range is inclusive as opposed to standard python indexing.
    638            `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6.
    639            If not specified, self.prange or all timeslices are used.
    640        silent : bool
    641            Decides whether output is printed to the standard output.
    642        '''
    643        if self.N != 1:
    644            raise Exception("Correlator must be projected before fitting")
    645
    646        if fitrange is None:
    647            if self.prange:
    648                fitrange = self.prange
    649            else:
    650                fitrange = [0, self.T - 1]
    651        else:
    652            if not isinstance(fitrange, list):
    653                raise Exception("fitrange has to be a list with two elements")
    654            if len(fitrange) != 2:
    655                raise Exception("fitrange has to have exactly two elements [fit_start, fit_stop]")
    656
    657        xs = [x for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None]
    658        ys = [self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None]
    659        result = least_squares(xs, ys, function, silent=silent, **kwargs)
    660        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

    662    def plateau(self, plateau_range=None, method="fit", auto_gamma=False):
    663        """ Extract a plateau value from a Corr object
    664
    665        Parameters
    666        ----------
    667        plateau_range : list
    668            list with two entries, indicating the first and the last timeslice
    669            of the plateau region.
    670        method : str
    671            method to extract the plateau.
    672                'fit' fits a constant to the plateau region
    673                'avg', 'average' or 'mean' just average over the given timeslices.
    674        auto_gamma : bool
    675            apply gamma_method with default parameters to the Corr. Defaults to None
    676        """
    677        if not plateau_range:
    678            if self.prange:
    679                plateau_range = self.prange
    680            else:
    681                raise Exception("no plateau range provided")
    682        if self.N != 1:
    683            raise Exception("Correlator must be projected before getting a plateau.")
    684        if(all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])):
    685            raise Exception("plateau is undefined at all timeslices in plateaurange.")
    686        if auto_gamma:
    687            self.gamma_method()
    688        if method == "fit":
    689            def const_func(a, t):
    690                return a[0]
    691            return self.fit(const_func, plateau_range)[0]
    692        elif method in ["avg", "average", "mean"]:
    693            returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None])
    694            return returnvalue
    695
    696        else:
    697            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

    699    def set_prange(self, prange):
    700        """Sets the attribute prange of the Corr object."""
    701        if not len(prange) == 2:
    702            raise Exception("prange must be a list or array with two values")
    703        if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))):
    704            raise Exception("Start and end point must be integers")
    705        if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] < prange[1]):
    706            raise Exception("Start and end point must define a range in the interval 0,T")
    707
    708        self.prange = prange
    709        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 ):
View Source

    711    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):
    712        """Plots the correlator using the tag of the correlator as label if available.
    713
    714        Parameters
    715        ----------
    716        x_range : list
    717            list of two values, determining the range of the x-axis e.g. [4, 8]
    718        comp : Corr or list of Corr
    719            Correlator or list of correlators which are plotted for comparison.
    720            The tags of these correlators are used as labels if available.
    721        logscale : bool
    722            Sets y-axis to logscale
    723        plateau : Obs
    724            Plateau value to be visualized in the figure
    725        fit_res : Fit_result
    726            Fit_result object to be visualized
    727        ylabel : str
    728            Label for the y-axis
    729        save : str
    730            path to file in which the figure should be saved
    731        auto_gamma : bool
    732            Apply the gamma method with standard parameters to all correlators and plateau values before plotting.
    733        hide_sigma : float
    734            Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors.
    735        references : list
    736            List of floating point values that are displayed as horizontal lines for reference.
    737        """
    738        if self.N != 1:
    739            raise Exception("Correlator must be projected before plotting")
    740
    741        if auto_gamma:
    742            self.gamma_method()
    743
    744        if x_range is None:
    745            x_range = [0, self.T - 1]
    746
    747        fig = plt.figure()
    748        ax1 = fig.add_subplot(111)
    749
    750        x, y, y_err = self.plottable()
    751        if hide_sigma:
    752            hide_from = np.argmax((hide_sigma * np.array(y_err)) > np.abs(y)) - 1
    753        else:
    754            hide_from = None
    755        ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag)
    756        if logscale:
    757            ax1.set_yscale('log')
    758        else:
    759            if y_range is None:
    760                try:
    761                    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)])
    762                    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)])
    763                    ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)])
    764                except Exception:
    765                    pass
    766            else:
    767                ax1.set_ylim(y_range)
    768        if comp:
    769            if isinstance(comp, (Corr, list)):
    770                for corr in comp if isinstance(comp, list) else [comp]:
    771                    if auto_gamma:
    772                        corr.gamma_method()
    773                    x, y, y_err = corr.plottable()
    774                    if hide_sigma:
    775                        hide_from = np.argmax((hide_sigma * np.array(y_err)) > np.abs(y)) - 1
    776                    else:
    777                        hide_from = None
    778                    plt.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor'])
    779            else:
    780                raise Exception("'comp' must be a correlator or a list of correlators.")
    781
    782        if plateau:
    783            if isinstance(plateau, Obs):
    784                if auto_gamma:
    785                    plateau.gamma_method()
    786                ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau))
    787                ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-')
    788            else:
    789                raise Exception("'plateau' must be an Obs")
    790
    791        if references:
    792            if isinstance(references, list):
    793                for ref in references:
    794                    ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--')
    795            else:
    796                raise Exception("'references' must be a list of floating pint values.")
    797
    798        if self.prange:
    799            ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',')
    800            ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',')
    801
    802        if fit_res:
    803            x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05)
    804            ax1.plot(x_samples,
    805                     fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples),
    806                     ls='-', marker=',', lw=2)
    807
    808        ax1.set_xlabel(r'$x_0 / a$')
    809        if ylabel:
    810            ax1.set_ylabel(ylabel)
    811        ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5])
    812
    813        handles, labels = ax1.get_legend_handles_labels()
    814        if labels:
    815            ax1.legend()
    816        plt.draw()
    817
    818        if save:
    819            if isinstance(save, str):
    820                fig.savefig(save)
    821            else:
    822                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.

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

    824    def spaghetti_plot(self, logscale=True):
    825        """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations.
    826
    827        Parameters
    828        ----------
    829        logscale : bool
    830            Determines whether the scale of the y-axis is logarithmic or standard.
    831        """
    832        if self.N != 1:
    833            raise Exception("Correlator needs to be projected first.")
    834
    835        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]))
    836        x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None]
    837
    838        for name in mc_names:
    839            data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T
    840
    841            fig = plt.figure()
    842            ax = fig.add_subplot(111)
    843            for dat in data:
    844                ax.plot(x0_vals, dat, ls='-', marker='')
    845
    846            if logscale is True:
    847                ax.set_yscale('log')
    848
    849            ax.set_xlabel(r'$x_0 / a$')
    850            plt.title(name)
    851            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

    853    def dump(self, filename, datatype="json.gz", **kwargs):
    854        """Dumps the Corr into a file of chosen type
    855        Parameters
    856        ----------
    857        filename : str
    858            Name of the file to be saved.
    859        datatype : str
    860            Format of the exported file. Supported formats include
    861            "json.gz" and "pickle"
    862        path : str
    863            specifies a custom path for the file (default '.')
    864        """
    865        if datatype == "json.gz":
    866            from .input.json import dump_to_json
    867            if 'path' in kwargs:
    868                file_name = kwargs.get('path') + '/' + filename
    869            else:
    870                file_name = filename
    871            dump_to_json(self, file_name)
    872        elif datatype == "pickle":
    873            dump_object(self, filename, **kwargs)
    874        else:
    875            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, range=[0, None]):
View Source

    877    def print(self, range=[0, None]):
    878        print(self.__repr__(range))

def sqrt(self):
View Source

    1040    def sqrt(self):
    1041        return self**0.5

def log(self):
View Source

    1043    def log(self):
    1044        newcontent = [None if (item is None) else np.log(item) for item in self.content]
    1045        return Corr(newcontent, prange=self.prange)

def exp(self):
View Source

    1047    def exp(self):
    1048        newcontent = [None if (item is None) else np.exp(item) for item in self.content]
    1049        return Corr(newcontent, prange=self.prange)

def sin(self):
View Source

    1062    def sin(self):
    1063        return self._apply_func_to_corr(np.sin)

def cos(self):
View Source

    1065    def cos(self):
    1066        return self._apply_func_to_corr(np.cos)

def tan(self):
View Source

    1068    def tan(self):
    1069        return self._apply_func_to_corr(np.tan)

def sinh(self):
View Source

    1071    def sinh(self):
    1072        return self._apply_func_to_corr(np.sinh)

def cosh(self):
View Source

    1074    def cosh(self):
    1075        return self._apply_func_to_corr(np.cosh)

def tanh(self):
View Source

    1077    def tanh(self):
    1078        return self._apply_func_to_corr(np.tanh)

def arcsin(self):
View Source

    1080    def arcsin(self):
    1081        return self._apply_func_to_corr(np.arcsin)

def arccos(self):
View Source

    1083    def arccos(self):
    1084        return self._apply_func_to_corr(np.arccos)

def arctan(self):
View Source

    1086    def arctan(self):
    1087        return self._apply_func_to_corr(np.arctan)

def arcsinh(self):
View Source

    1089    def arcsinh(self):
    1090        return self._apply_func_to_corr(np.arcsinh)

def arccosh(self):
View Source

    1092    def arccosh(self):
    1093        return self._apply_func_to_corr(np.arccosh)

def arctanh(self):
View Source

    1095    def arctanh(self):
    1096        return self._apply_func_to_corr(np.arctanh)

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

    1131    def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
    1132        r''' Project large correlation matrix to lowest states
    1133
    1134        This method can be used to reduce the size of an (N x N) correlation matrix
    1135        to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise
    1136        is still small.
    1137
    1138        Parameters
    1139        ----------
    1140        Ntrunc: int
    1141            Rank of the target matrix.
    1142        tproj: int
    1143            Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method.
    1144            The default value is 3.
    1145        t0proj: int
    1146            Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly
    1147            discouraged for O(a) improved theories, since the correctness of the procedure
    1148            cannot be granted in this case. The default value is 2.
    1149        basematrix : Corr
    1150            Correlation matrix that is used to determine the eigenvectors of the
    1151            lowest states based on a GEVP. basematrix is taken to be the Corr itself if
    1152            is is not specified.
    1153
    1154        Notes
    1155        -----
    1156        We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving
    1157        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}$
    1158        and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the
    1159        resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via
    1160        $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large
    1161        correlation matrix and to remove some noise that is added by irrelevant operators.
    1162        This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated
    1163        bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$.
    1164        '''
    1165
    1166        if self.N == 1:
    1167            raise Exception('Method cannot be applied to one-dimensional correlators.')
    1168        if basematrix is None:
    1169            basematrix = self
    1170        if Ntrunc >= basematrix.N:
    1171            raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N))
    1172        if basematrix.N != self.N:
    1173            raise Exception('basematrix and targetmatrix have to be of the same size.')
    1174
    1175        evecs = []
    1176        for i in range(Ntrunc):
    1177            evecs.append(basematrix.GEVP(t0proj, tproj, state=i, sorted_list=None))
    1178
    1179        tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object)
    1180        rmat = []
    1181        for t in range(basematrix.T):
    1182            for i in range(Ntrunc):
    1183                for j in range(Ntrunc):
    1184                    tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j]
    1185            rmat.append(np.copy(tmpmat))
    1186
    1187        newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)]
    1188        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)$.