pyerrors.correlators

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

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

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

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

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

Initialization

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

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

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

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

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

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

Initialize a Corr object.

Parameters
  • data_input (list or array): list of Obs or list of arrays of Obs or array of Corrs (see class docstring for details).
  • padding (list, optional): List with two entries where the first labels the padding at the front of the correlator and the second the padding at the back.
  • prange (list, optional): List containing the first and last timeslice of the plateau region identified for this correlator.
tag
content
T
prange
reweighted
def gamma_method(self, **kwargs):
141    def gamma_method(self, **kwargs):
142        """Apply the gamma method to the content of the Corr."""
143        for item in self.content:
144            if not (item is None):
145                if self.N == 1:
146                    item[0].gamma_method(**kwargs)
147                else:
148                    for i in range(self.N):
149                        for j in range(self.N):
150                            item[i, j].gamma_method(**kwargs)

Apply the gamma method to the content of the Corr.

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

Apply the gamma method to the content of the Corr.

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

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

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

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

Symmetrize the correlator around x0=0.

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

Anti-symmetrize the correlator around x0=0.

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

Checks whether a correlator matrices is symmetric on every timeslice.

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

Calculates the per-timeslice trace of a correlator matrix.

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

Symmetrizes the correlator matrices on every timeslice.

def GEVP(self, t0, ts=None, sort='Eigenvalue', **kwargs):
302    def GEVP(self, t0, ts=None, sort="Eigenvalue", **kwargs):
303        r'''Solve the generalized eigenvalue problem on the correlator matrix and returns the corresponding eigenvectors.
304
305        The eigenvectors are sorted according to the descending eigenvalues, the zeroth eigenvector(s) correspond to the
306        largest eigenvalue(s). The eigenvector(s) for the individual states can be accessed via slicing
307        ```python
308        C.GEVP(t0=2)[0]  # Ground state vector(s)
309        C.GEVP(t0=2)[:3]  # Vectors for the lowest three states
310        ```
311
312        Parameters
313        ----------
314        t0 : int
315            The time t0 for the right hand side of the GEVP according to $G(t)v_i=\lambda_i G(t_0)v_i$
316        ts : int
317            fixed time $G(t_s)v_i=\lambda_i G(t_0)v_i$ if sort=None.
318            If sort="Eigenvector" it gives a reference point for the sorting method.
319        sort : string
320            If this argument is set, a list of self.T vectors per state is returned. If it is set to None, only one vector is returned.
321            - "Eigenvalue": The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice.
322            - "Eigenvector": Use the method described in arXiv:2004.10472 to find the set of v(t) belonging to the state.
323              The reference state is identified by its eigenvalue at $t=t_s$.
324
325        Other Parameters
326        ----------------
327        state : int
328           Returns only the vector(s) for a specified state. The lowest state is zero.
329        '''
330
331        if self.N == 1:
332            raise Exception("GEVP methods only works on correlator matrices and not single correlators.")
333        if ts is not None:
334            if (ts <= t0):
335                raise Exception("ts has to be larger than t0.")
336
337        if "sorted_list" in kwargs:
338            warnings.warn("Argument 'sorted_list' is deprecated, use 'sort' instead.", DeprecationWarning)
339            sort = kwargs.get("sorted_list")
340
341        if self.is_matrix_symmetric():
342            symmetric_corr = self
343        else:
344            symmetric_corr = self.matrix_symmetric()
345
346        G0 = np.vectorize(lambda x: x.value)(symmetric_corr[t0])
347        np.linalg.cholesky(G0)  # Check if matrix G0 is positive-semidefinite.
348
349        if sort is None:
350            if (ts is None):
351                raise Exception("ts is required if sort=None.")
352            if (self.content[t0] is None) or (self.content[ts] is None):
353                raise Exception("Corr not defined at t0/ts.")
354            Gt = np.vectorize(lambda x: x.value)(symmetric_corr[ts])
355            reordered_vecs = _GEVP_solver(Gt, G0)
356
357        elif sort in ["Eigenvalue", "Eigenvector"]:
358            if sort == "Eigenvalue" and ts is not None:
359                warnings.warn("ts has no effect when sorting by eigenvalue is chosen.", RuntimeWarning)
360            all_vecs = [None] * (t0 + 1)
361            for t in range(t0 + 1, self.T):
362                try:
363                    Gt = np.vectorize(lambda x: x.value)(symmetric_corr[t])
364                    all_vecs.append(_GEVP_solver(Gt, G0))
365                except Exception:
366                    all_vecs.append(None)
367            if sort == "Eigenvector":
368                if ts is None:
369                    raise Exception("ts is required for the Eigenvector sorting method.")
370                all_vecs = _sort_vectors(all_vecs, ts)
371
372            reordered_vecs = [[v[s] if v is not None else None for v in all_vecs] for s in range(self.N)]
373        else:
374            raise Exception("Unkown value for 'sort'.")
375
376        if "state" in kwargs:
377            return reordered_vecs[kwargs.get("state")]
378        else:
379            return reordered_vecs

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

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

C.GEVP(t0=2)[0]  # Ground state vector(s)
C.GEVP(t0=2)[:3]  # Vectors for the lowest three states
Parameters
  • t0 (int): The time t0 for the right hand side of the GEVP according to $G(t)v_i=\lambda_i G(t_0)v_i$
  • ts (int): fixed time $G(t_s)v_i=\lambda_i G(t_0)v_i$ if sort=None. If sort="Eigenvector" it gives a reference point for the sorting method.
  • sort (string): If this argument is set, a list of self.T vectors per state is returned. If it is set to None, only one vector is returned.
    • "Eigenvalue": The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice.
    • "Eigenvector": Use the method described in arXiv:2004.10472 to find the set of v(t) belonging to the state. The reference state is identified by its eigenvalue at $t=t_s$.
Other Parameters
  • state (int): Returns only the vector(s) for a specified state. The lowest state is zero.
def Eigenvalue(self, t0, ts=None, state=0, sort='Eigenvalue'):
381    def Eigenvalue(self, t0, ts=None, state=0, sort="Eigenvalue"):
382        """Determines the eigenvalue of the GEVP by solving and projecting the correlator
383
384        Parameters
385        ----------
386        state : int
387            The state one is interested in ordered by energy. The lowest state is zero.
388
389        All other parameters are identical to the ones of Corr.GEVP.
390        """
391        vec = self.GEVP(t0, ts=ts, sort=sort)[state]
392        return self.projected(vec)

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

Parameters
  • state (int): The state one is interested in ordered by energy. The lowest state is zero.
  • All other parameters are identical to the ones of Corr.GEVP.
def Hankel(self, N, periodic=False):
394    def Hankel(self, N, periodic=False):
395        """Constructs an NxN Hankel matrix
396
397        C(t) c(t+1) ... c(t+n-1)
398        C(t+1) c(t+2) ... c(t+n)
399        .................
400        C(t+(n-1)) c(t+n) ... c(t+2(n-1))
401
402        Parameters
403        ----------
404        N : int
405            Dimension of the Hankel matrix
406        periodic : bool, optional
407            determines whether the matrix is extended periodically
408        """
409
410        if self.N != 1:
411            raise Exception("Multi-operator Prony not implemented!")
412
413        array = np.empty([N, N], dtype="object")
414        new_content = []
415        for t in range(self.T):
416            new_content.append(array.copy())
417
418        def wrap(i):
419            while i >= self.T:
420                i -= self.T
421            return i
422
423        for t in range(self.T):
424            for i in range(N):
425                for j in range(N):
426                    if periodic:
427                        new_content[t][i, j] = self.content[wrap(t + i + j)][0]
428                    elif (t + i + j) >= self.T:
429                        new_content[t] = None
430                    else:
431                        new_content[t][i, j] = self.content[t + i + j][0]
432
433        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):
435    def roll(self, dt):
436        """Periodically shift the correlator by dt timeslices
437
438        Parameters
439        ----------
440        dt : int
441            number of timeslices
442        """
443        return Corr(list(np.roll(np.array(self.content, dtype=object), dt, axis=0)))

Periodically shift the correlator by dt timeslices

Parameters
  • dt (int): number of timeslices
def reverse(self):
445    def reverse(self):
446        """Reverse the time ordering of the Corr"""
447        return Corr(self.content[:: -1])

Reverse the time ordering of the Corr

def thin(self, spacing=2, offset=0):
449    def thin(self, spacing=2, offset=0):
450        """Thin out a correlator to suppress correlations
451
452        Parameters
453        ----------
454        spacing : int
455            Keep only every 'spacing'th entry of the correlator
456        offset : int
457            Offset the equal spacing
458        """
459        new_content = []
460        for t in range(self.T):
461            if (offset + t) % spacing != 0:
462                new_content.append(None)
463            else:
464                new_content.append(self.content[t])
465        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):
467    def correlate(self, partner):
468        """Correlate the correlator with another correlator or Obs
469
470        Parameters
471        ----------
472        partner : Obs or Corr
473            partner to correlate the correlator with.
474            Can either be an Obs which is correlated with all entries of the
475            correlator or a Corr of same length.
476        """
477        if self.N != 1:
478            raise Exception("Only one-dimensional correlators can be safely correlated.")
479        new_content = []
480        for x0, t_slice in enumerate(self.content):
481            if _check_for_none(self, t_slice):
482                new_content.append(None)
483            else:
484                if isinstance(partner, Corr):
485                    if _check_for_none(partner, partner.content[x0]):
486                        new_content.append(None)
487                    else:
488                        new_content.append(np.array([correlate(o, partner.content[x0][0]) for o in t_slice]))
489                elif isinstance(partner, Obs):  # Should this include CObs?
490                    new_content.append(np.array([correlate(o, partner) for o in t_slice]))
491                else:
492                    raise Exception("Can only correlate with an Obs or a Corr.")
493
494        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):
496    def reweight(self, weight, **kwargs):
497        """Reweight the correlator.
498
499        Parameters
500        ----------
501        weight : Obs
502            Reweighting factor. An Observable that has to be defined on a superset of the
503            configurations in obs[i].idl for all i.
504        all_configs : bool
505            if True, the reweighted observables are normalized by the average of
506            the reweighting factor on all configurations in weight.idl and not
507            on the configurations in obs[i].idl.
508        """
509        if self.N != 1:
510            raise Exception("Reweighting only implemented for one-dimensional correlators.")
511        new_content = []
512        for t_slice in self.content:
513            if _check_for_none(self, t_slice):
514                new_content.append(None)
515            else:
516                new_content.append(np.array(reweight(weight, t_slice, **kwargs)))
517        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):
519    def T_symmetry(self, partner, parity=+1):
520        """Return the time symmetry average of the correlator and its partner
521
522        Parameters
523        ----------
524        partner : Corr
525            Time symmetry partner of the Corr
526        parity : int
527            Parity quantum number of the correlator, can be +1 or -1
528        """
529        if self.N != 1:
530            raise Exception("T_symmetry only implemented for one-dimensional correlators.")
531        if not isinstance(partner, Corr):
532            raise Exception("T partner has to be a Corr object.")
533        if parity not in [+1, -1]:
534            raise Exception("Parity has to be +1 or -1.")
535        T_partner = parity * partner.reverse()
536
537        t_slices = []
538        test = (self - T_partner)
539        test.gamma_method()
540        for x0, t_slice in enumerate(test.content):
541            if t_slice is not None:
542                if not t_slice[0].is_zero_within_error(5):
543                    t_slices.append(x0)
544        if t_slices:
545            warnings.warn("T symmetry partners do not agree within 5 sigma on time slices " + str(t_slices) + ".", RuntimeWarning)
546
547        return (self + T_partner) / 2

Return the time symmetry average of the correlator and its partner

Parameters
  • partner (Corr): Time symmetry partner of the Corr
  • parity (int): Parity quantum number of the correlator, can be +1 or -1
def deriv(self, variant='symmetric'):
549    def deriv(self, variant="symmetric"):
550        """Return the first derivative of the correlator with respect to x0.
551
552        Parameters
553        ----------
554        variant : str
555            decides which definition of the finite differences derivative is used.
556            Available choice: symmetric, forward, backward, improved, log, default: symmetric
557        """
558        if self.N != 1:
559            raise Exception("deriv only implemented for one-dimensional correlators.")
560        if variant == "symmetric":
561            newcontent = []
562            for t in range(1, self.T - 1):
563                if (self.content[t - 1] is None) or (self.content[t + 1] is None):
564                    newcontent.append(None)
565                else:
566                    newcontent.append(0.5 * (self.content[t + 1] - self.content[t - 1]))
567            if (all([x is None for x in newcontent])):
568                raise Exception('Derivative is undefined at all timeslices')
569            return Corr(newcontent, padding=[1, 1])
570        elif variant == "forward":
571            newcontent = []
572            for t in range(self.T - 1):
573                if (self.content[t] is None) or (self.content[t + 1] is None):
574                    newcontent.append(None)
575                else:
576                    newcontent.append(self.content[t + 1] - self.content[t])
577            if (all([x is None for x in newcontent])):
578                raise Exception("Derivative is undefined at all timeslices")
579            return Corr(newcontent, padding=[0, 1])
580        elif variant == "backward":
581            newcontent = []
582            for t in range(1, self.T):
583                if (self.content[t - 1] is None) or (self.content[t] is None):
584                    newcontent.append(None)
585                else:
586                    newcontent.append(self.content[t] - self.content[t - 1])
587            if (all([x is None for x in newcontent])):
588                raise Exception("Derivative is undefined at all timeslices")
589            return Corr(newcontent, padding=[1, 0])
590        elif variant == "improved":
591            newcontent = []
592            for t in range(2, self.T - 2):
593                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):
594                    newcontent.append(None)
595                else:
596                    newcontent.append((1 / 12) * (self.content[t - 2] - 8 * self.content[t - 1] + 8 * self.content[t + 1] - self.content[t + 2]))
597            if (all([x is None for x in newcontent])):
598                raise Exception('Derivative is undefined at all timeslices')
599            return Corr(newcontent, padding=[2, 2])
600        elif variant == 'log':
601            newcontent = []
602            for t in range(self.T):
603                if (self.content[t] is None) or (self.content[t] <= 0):
604                    newcontent.append(None)
605                else:
606                    newcontent.append(np.log(self.content[t]))
607            if (all([x is None for x in newcontent])):
608                raise Exception("Log is undefined at all timeslices")
609            logcorr = Corr(newcontent)
610            return self * logcorr.deriv('symmetric')
611        else:
612            raise Exception("Unknown variant.")

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

Parameters
  • variant (str): decides which definition of the finite differences derivative is used. Available choice: symmetric, forward, backward, improved, log, default: symmetric
def second_deriv(self, variant='symmetric'):
614    def second_deriv(self, variant="symmetric"):
615        r"""Return the second derivative of the correlator with respect to x0.
616
617        Parameters
618        ----------
619        variant : str
620            decides which definition of the finite differences derivative is used.
621            Available choice:
622                - symmetric (default)
623                    $$\tilde{\partial}^2_0 f(x_0) = f(x_0+1)-2f(x_0)+f(x_0-1)$$
624                - big_symmetric
625                    $$\partial^2_0 f(x_0) = \frac{f(x_0+2)-2f(x_0)+f(x_0-2)}{4}$$
626                - improved
627                    $$\partial^2_0 f(x_0) = \frac{-f(x_0+2) + 16 * f(x_0+1) - 30 * f(x_0) + 16 * f(x_0-1) - f(x_0-2)}{12}$$
628                - log
629                    $$f(x) = \tilde{\partial}^2_0 log(f(x_0))+(\tilde{\partial}_0 log(f(x_0)))^2$$
630        """
631        if self.N != 1:
632            raise Exception("second_deriv only implemented for one-dimensional correlators.")
633        if variant == "symmetric":
634            newcontent = []
635            for t in range(1, self.T - 1):
636                if (self.content[t - 1] is None) or (self.content[t + 1] is None):
637                    newcontent.append(None)
638                else:
639                    newcontent.append((self.content[t + 1] - 2 * self.content[t] + self.content[t - 1]))
640            if (all([x is None for x in newcontent])):
641                raise Exception("Derivative is undefined at all timeslices")
642            return Corr(newcontent, padding=[1, 1])
643        elif variant == "big_symmetric":
644            newcontent = []
645            for t in range(2, self.T - 2):
646                if (self.content[t - 2] is None) or (self.content[t + 2] is None):
647                    newcontent.append(None)
648                else:
649                    newcontent.append((self.content[t + 2] - 2 * self.content[t] + self.content[t - 2]) / 4)
650            if (all([x is None for x in newcontent])):
651                raise Exception("Derivative is undefined at all timeslices")
652            return Corr(newcontent, padding=[2, 2])
653        elif variant == "improved":
654            newcontent = []
655            for t in range(2, self.T - 2):
656                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):
657                    newcontent.append(None)
658                else:
659                    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]))
660            if (all([x is None for x in newcontent])):
661                raise Exception("Derivative is undefined at all timeslices")
662            return Corr(newcontent, padding=[2, 2])
663        elif variant == 'log':
664            newcontent = []
665            for t in range(self.T):
666                if (self.content[t] is None) or (self.content[t] <= 0):
667                    newcontent.append(None)
668                else:
669                    newcontent.append(np.log(self.content[t]))
670            if (all([x is None for x in newcontent])):
671                raise Exception("Log is undefined at all timeslices")
672            logcorr = Corr(newcontent)
673            return self * (logcorr.second_deriv('symmetric') + (logcorr.deriv('symmetric'))**2)
674        else:
675            raise Exception("Unknown variant.")

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

Parameters
  • variant (str): decides which definition of the finite differences derivative is used. Available choice: - symmetric (default) $$\tilde{\partial}^2_0 f(x_0) = f(x_0+1)-2f(x_0)+f(x_0-1)$$ - big_symmetric $$\partial^2_0 f(x_0) = \frac{f(x_0+2)-2f(x_0)+f(x_0-2)}{4}$$ - improved $$\partial^2_0 f(x_0) = \frac{-f(x_0+2) + 16 * f(x_0+1) - 30 * f(x_0) + 16 * f(x_0-1) - f(x_0-2)}{12}$$ - log $$f(x) = \tilde{\partial}^2_0 log(f(x_0))+(\tilde{\partial}_0 log(f(x_0)))^2$$
def m_eff(self, variant='log', guess=1.0):
677    def m_eff(self, variant='log', guess=1.0):
678        """Returns the effective mass of the correlator as correlator object
679
680        Parameters
681        ----------
682        variant : str
683            log : uses the standard effective mass log(C(t) / C(t+1))
684            cosh, periodic : Use periodicity of the correlator by solving C(t) / C(t+1) = cosh(m * (t - T/2)) / cosh(m * (t + 1 - T/2)) for m.
685            sinh : Use anti-periodicity of the correlator by solving C(t) / C(t+1) = sinh(m * (t - T/2)) / sinh(m * (t + 1 - T/2)) for m.
686            See, e.g., arXiv:1205.5380
687            arccosh : Uses the explicit form of the symmetrized correlator (not recommended)
688            logsym: uses the symmetric effective mass log(C(t-1) / C(t+1))/2
689        guess : float
690            guess for the root finder, only relevant for the root variant
691        """
692        if self.N != 1:
693            raise Exception('Correlator must be projected before getting m_eff')
694        if variant == 'log':
695            newcontent = []
696            for t in range(self.T - 1):
697                if ((self.content[t] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0):
698                    newcontent.append(None)
699                elif self.content[t][0].value / self.content[t + 1][0].value < 0:
700                    newcontent.append(None)
701                else:
702                    newcontent.append(self.content[t] / self.content[t + 1])
703            if (all([x is None for x in newcontent])):
704                raise Exception('m_eff is undefined at all timeslices')
705
706            return np.log(Corr(newcontent, padding=[0, 1]))
707
708        elif variant == 'logsym':
709            newcontent = []
710            for t in range(1, self.T - 1):
711                if ((self.content[t - 1] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0):
712                    newcontent.append(None)
713                elif self.content[t - 1][0].value / self.content[t + 1][0].value < 0:
714                    newcontent.append(None)
715                else:
716                    newcontent.append(self.content[t - 1] / self.content[t + 1])
717            if (all([x is None for x in newcontent])):
718                raise Exception('m_eff is undefined at all timeslices')
719
720            return np.log(Corr(newcontent, padding=[1, 1])) / 2
721
722        elif variant in ['periodic', 'cosh', 'sinh']:
723            if variant in ['periodic', 'cosh']:
724                func = anp.cosh
725            else:
726                func = anp.sinh
727
728            def root_function(x, d):
729                return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d
730
731            newcontent = []
732            for t in range(self.T - 1):
733                if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t + 1][0].value == 0):
734                    newcontent.append(None)
735                # Fill the two timeslices in the middle of the lattice with their predecessors
736                elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]:
737                    newcontent.append(newcontent[-1])
738                elif self.content[t][0].value / self.content[t + 1][0].value < 0:
739                    newcontent.append(None)
740                else:
741                    newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess)))
742            if (all([x is None for x in newcontent])):
743                raise Exception('m_eff is undefined at all timeslices')
744
745            return Corr(newcontent, padding=[0, 1])
746
747        elif variant == 'arccosh':
748            newcontent = []
749            for t in range(1, self.T - 1):
750                if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t - 1] is None) or (self.content[t][0].value == 0):
751                    newcontent.append(None)
752                else:
753                    newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t]))
754            if (all([x is None for x in newcontent])):
755                raise Exception("m_eff is undefined at all timeslices")
756            return np.arccosh(Corr(newcontent, padding=[1, 1]))
757
758        else:
759            raise Exception('Unknown variant.')

Returns the effective mass of the correlator as correlator object

Parameters
  • variant (str): log : uses the standard effective mass log(C(t) / C(t+1)) cosh, periodic : Use periodicity of the correlator by solving C(t) / C(t+1) = cosh(m * (t - T/2)) / cosh(m * (t + 1 - T/2)) for m. sinh : Use anti-periodicity of the correlator by solving C(t) / C(t+1) = sinh(m * (t - T/2)) / sinh(m * (t + 1 - T/2)) for m. See, e.g., arXiv:1205.5380 arccosh : Uses the explicit form of the symmetrized correlator (not recommended) logsym: uses the symmetric effective mass log(C(t-1) / C(t+1))/2
  • guess (float): guess for the root finder, only relevant for the root variant
def fit(self, function, fitrange=None, silent=False, **kwargs):
761    def fit(self, function, fitrange=None, silent=False, **kwargs):
762        r'''Fits function to the data
763
764        Parameters
765        ----------
766        function : obj
767            function to fit to the data. See fits.least_squares for details.
768        fitrange : list
769            Two element list containing the timeslices on which the fit is supposed to start and stop.
770            Caution: This range is inclusive as opposed to standard python indexing.
771            `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6.
772            If not specified, self.prange or all timeslices are used.
773        silent : bool
774            Decides whether output is printed to the standard output.
775        '''
776        if self.N != 1:
777            raise Exception("Correlator must be projected before fitting")
778
779        if fitrange is None:
780            if self.prange:
781                fitrange = self.prange
782            else:
783                fitrange = [0, self.T - 1]
784        else:
785            if not isinstance(fitrange, list):
786                raise Exception("fitrange has to be a list with two elements")
787            if len(fitrange) != 2:
788                raise Exception("fitrange has to have exactly two elements [fit_start, fit_stop]")
789
790        xs = np.array([x for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None])
791        ys = np.array([self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None])
792        result = least_squares(xs, ys, function, silent=silent, **kwargs)
793        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):
795    def plateau(self, plateau_range=None, method="fit", auto_gamma=False):
796        """ Extract a plateau value from a Corr object
797
798        Parameters
799        ----------
800        plateau_range : list
801            list with two entries, indicating the first and the last timeslice
802            of the plateau region.
803        method : str
804            method to extract the plateau.
805                'fit' fits a constant to the plateau region
806                'avg', 'average' or 'mean' just average over the given timeslices.
807        auto_gamma : bool
808            apply gamma_method with default parameters to the Corr. Defaults to None
809        """
810        if not plateau_range:
811            if self.prange:
812                plateau_range = self.prange
813            else:
814                raise Exception("no plateau range provided")
815        if self.N != 1:
816            raise Exception("Correlator must be projected before getting a plateau.")
817        if (all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])):
818            raise Exception("plateau is undefined at all timeslices in plateaurange.")
819        if auto_gamma:
820            self.gamma_method()
821        if method == "fit":
822            def const_func(a, t):
823                return a[0]
824            return self.fit(const_func, plateau_range)[0]
825        elif method in ["avg", "average", "mean"]:
826            returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None])
827            return returnvalue
828
829        else:
830            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):
832    def set_prange(self, prange):
833        """Sets the attribute prange of the Corr object."""
834        if not len(prange) == 2:
835            raise Exception("prange must be a list or array with two values")
836        if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))):
837            raise Exception("Start and end point must be integers")
838        if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] < prange[1]):
839            raise Exception("Start and end point must define a range in the interval 0,T")
840
841        self.prange = prange
842        return

Sets the attribute prange of the Corr object.

def show( self, x_range=None, comp=None, y_range=None, logscale=False, plateau=None, fit_res=None, fit_key=None, ylabel=None, save=None, auto_gamma=False, hide_sigma=None, references=None, title=None):
844    def show(self, x_range=None, comp=None, y_range=None, logscale=False, plateau=None, fit_res=None, fit_key=None, ylabel=None, save=None, auto_gamma=False, hide_sigma=None, references=None, title=None):
845        """Plots the correlator using the tag of the correlator as label if available.
846
847        Parameters
848        ----------
849        x_range : list
850            list of two values, determining the range of the x-axis e.g. [4, 8].
851        comp : Corr or list of Corr
852            Correlator or list of correlators which are plotted for comparison.
853            The tags of these correlators are used as labels if available.
854        logscale : bool
855            Sets y-axis to logscale.
856        plateau : Obs
857            Plateau value to be visualized in the figure.
858        fit_res : Fit_result
859            Fit_result object to be visualized.
860        fit_key : str
861            Key for the fit function in Fit_result.fit_function (for combined fits).
862        ylabel : str
863            Label for the y-axis.
864        save : str
865            path to file in which the figure should be saved.
866        auto_gamma : bool
867            Apply the gamma method with standard parameters to all correlators and plateau values before plotting.
868        hide_sigma : float
869            Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors.
870        references : list
871            List of floating point values that are displayed as horizontal lines for reference.
872        title : string
873            Optional title of the figure.
874        """
875        if self.N != 1:
876            raise Exception("Correlator must be projected before plotting")
877
878        if auto_gamma:
879            self.gamma_method()
880
881        if x_range is None:
882            x_range = [0, self.T - 1]
883
884        fig = plt.figure()
885        ax1 = fig.add_subplot(111)
886
887        x, y, y_err = self.plottable()
888        if hide_sigma:
889            hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
890        else:
891            hide_from = None
892        ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag)
893        if logscale:
894            ax1.set_yscale('log')
895        else:
896            if y_range is None:
897                try:
898                    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)])
899                    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)])
900                    ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)])
901                except Exception:
902                    pass
903            else:
904                ax1.set_ylim(y_range)
905        if comp:
906            if isinstance(comp, (Corr, list)):
907                for corr in comp if isinstance(comp, list) else [comp]:
908                    if auto_gamma:
909                        corr.gamma_method()
910                    x, y, y_err = corr.plottable()
911                    if hide_sigma:
912                        hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
913                    else:
914                        hide_from = None
915                    ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor'])
916            else:
917                raise Exception("'comp' must be a correlator or a list of correlators.")
918
919        if plateau:
920            if isinstance(plateau, Obs):
921                if auto_gamma:
922                    plateau.gamma_method()
923                ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau))
924                ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-')
925            else:
926                raise Exception("'plateau' must be an Obs")
927
928        if references:
929            if isinstance(references, list):
930                for ref in references:
931                    ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--')
932            else:
933                raise Exception("'references' must be a list of floating pint values.")
934
935        if self.prange:
936            ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',', color="black", zorder=0)
937            ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',', color="black", zorder=0)
938
939        if fit_res:
940            x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05)
941            if isinstance(fit_res.fit_function, dict):
942                if fit_key:
943                    ax1.plot(x_samples, fit_res.fit_function[fit_key]([o.value for o in fit_res.fit_parameters], x_samples), ls='-', marker=',', lw=2)
944                else:
945                    raise ValueError("Please provide a 'fit_key' for visualizing combined fits.")
946            else:
947                ax1.plot(x_samples, fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples), ls='-', marker=',', lw=2)
948
949        ax1.set_xlabel(r'$x_0 / a$')
950        if ylabel:
951            ax1.set_ylabel(ylabel)
952        ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5])
953
954        handles, labels = ax1.get_legend_handles_labels()
955        if labels:
956            ax1.legend()
957
958        if title:
959            plt.title(title)
960
961        plt.draw()
962
963        if save:
964            if isinstance(save, str):
965                fig.savefig(save, bbox_inches='tight')
966            else:
967                raise Exception("'save' has to be a string.")

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

Parameters
  • x_range (list): list of two values, determining the range of the x-axis e.g. [4, 8].
  • comp (Corr or list of Corr): Correlator or list of correlators which are plotted for comparison. The tags of these correlators are used as labels if available.
  • logscale (bool): Sets y-axis to logscale.
  • plateau (Obs): Plateau value to be visualized in the figure.
  • fit_res (Fit_result): Fit_result object to be visualized.
  • fit_key (str): Key for the fit function in Fit_result.fit_function (for combined fits).
  • ylabel (str): Label for the y-axis.
  • save (str): path to file in which the figure should be saved.
  • auto_gamma (bool): Apply the gamma method with standard parameters to all correlators and plateau values before plotting.
  • hide_sigma (float): Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors.
  • references (list): List of floating point values that are displayed as horizontal lines for reference.
  • title (string): Optional title of the figure.
def spaghetti_plot(self, logscale=True):
969    def spaghetti_plot(self, logscale=True):
970        """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations.
971
972        Parameters
973        ----------
974        logscale : bool
975            Determines whether the scale of the y-axis is logarithmic or standard.
976        """
977        if self.N != 1:
978            raise Exception("Correlator needs to be projected first.")
979
980        mc_names = list(set([item for sublist in [sum(map(o[0].e_content.get, o[0].mc_names), []) for o in self.content if o is not None] for item in sublist]))
981        x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None]
982
983        for name in mc_names:
984            data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T
985
986            fig = plt.figure()
987            ax = fig.add_subplot(111)
988            for dat in data:
989                ax.plot(x0_vals, dat, ls='-', marker='')
990
991            if logscale is True:
992                ax.set_yscale('log')
993
994            ax.set_xlabel(r'$x_0 / a$')
995            plt.title(name)
996            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):
 998    def dump(self, filename, datatype="json.gz", **kwargs):
 999        """Dumps the Corr into a file of chosen type
1000        Parameters
1001        ----------
1002        filename : str
1003            Name of the file to be saved.
1004        datatype : str
1005            Format of the exported file. Supported formats include
1006            "json.gz" and "pickle"
1007        path : str
1008            specifies a custom path for the file (default '.')
1009        """
1010        if datatype == "json.gz":
1011            from .input.json import dump_to_json
1012            if 'path' in kwargs:
1013                file_name = kwargs.get('path') + '/' + filename
1014            else:
1015                file_name = filename
1016            dump_to_json(self, file_name)
1017        elif datatype == "pickle":
1018            dump_object(self, filename, **kwargs)
1019        else:
1020            raise Exception("Unknown datatype " + str(datatype))

Dumps the Corr into a file of chosen type

Parameters
  • filename (str): Name of the file to be saved.
  • datatype (str): Format of the exported file. Supported formats include "json.gz" and "pickle"
  • path (str): specifies a custom path for the file (default '.')
def print(self, print_range=None):
1022    def print(self, print_range=None):
1023        print(self.__repr__(print_range))
def sqrt(self):
1232    def sqrt(self):
1233        return self ** 0.5
def log(self):
1235    def log(self):
1236        newcontent = [None if _check_for_none(self, item) else np.log(item) for item in self.content]
1237        return Corr(newcontent, prange=self.prange)
def exp(self):
1239    def exp(self):
1240        newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content]
1241        return Corr(newcontent, prange=self.prange)
def sin(self):
1256    def sin(self):
1257        return self._apply_func_to_corr(np.sin)
def cos(self):
1259    def cos(self):
1260        return self._apply_func_to_corr(np.cos)
def tan(self):
1262    def tan(self):
1263        return self._apply_func_to_corr(np.tan)
def sinh(self):
1265    def sinh(self):
1266        return self._apply_func_to_corr(np.sinh)
def cosh(self):
1268    def cosh(self):
1269        return self._apply_func_to_corr(np.cosh)
def tanh(self):
1271    def tanh(self):
1272        return self._apply_func_to_corr(np.tanh)
def arcsin(self):
1274    def arcsin(self):
1275        return self._apply_func_to_corr(np.arcsin)
def arccos(self):
1277    def arccos(self):
1278        return self._apply_func_to_corr(np.arccos)
def arctan(self):
1280    def arctan(self):
1281        return self._apply_func_to_corr(np.arctan)
def arcsinh(self):
1283    def arcsinh(self):
1284        return self._apply_func_to_corr(np.arcsinh)
def arccosh(self):
1286    def arccosh(self):
1287        return self._apply_func_to_corr(np.arccosh)
def arctanh(self):
1289    def arctanh(self):
1290        return self._apply_func_to_corr(np.arctanh)
real
imag
def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
1325    def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
1326        r''' Project large correlation matrix to lowest states
1327
1328        This method can be used to reduce the size of an (N x N) correlation matrix
1329        to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise
1330        is still small.
1331
1332        Parameters
1333        ----------
1334        Ntrunc: int
1335            Rank of the target matrix.
1336        tproj: int
1337            Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method.
1338            The default value is 3.
1339        t0proj: int
1340            Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly
1341            discouraged for O(a) improved theories, since the correctness of the procedure
1342            cannot be granted in this case. The default value is 2.
1343        basematrix : Corr
1344            Correlation matrix that is used to determine the eigenvectors of the
1345            lowest states based on a GEVP. basematrix is taken to be the Corr itself if
1346            is is not specified.
1347
1348        Notes
1349        -----
1350        We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving
1351        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}$
1352        and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the
1353        resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via
1354        $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large
1355        correlation matrix and to remove some noise that is added by irrelevant operators.
1356        This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated
1357        bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$.
1358        '''
1359
1360        if self.N == 1:
1361            raise Exception('Method cannot be applied to one-dimensional correlators.')
1362        if basematrix is None:
1363            basematrix = self
1364        if Ntrunc >= basematrix.N:
1365            raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N))
1366        if basematrix.N != self.N:
1367            raise Exception('basematrix and targetmatrix have to be of the same size.')
1368
1369        evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc]
1370
1371        tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object)
1372        rmat = []
1373        for t in range(basematrix.T):
1374            for i in range(Ntrunc):
1375                for j in range(Ntrunc):
1376                    tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j]
1377            rmat.append(np.copy(tmpmat))
1378
1379        newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)]
1380        return Corr(newcontent)

Project large correlation matrix to lowest states

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

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

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

N