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

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

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

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

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

Initialize a Corr object.

Parameters
  • data_input (list or array): list of Obs or list of arrays of Obs or array of Corrs
  • padding (list, optional): List with two entries where the first labels the padding at the front of the correlator and the second the padding at the back.
  • prange (list, optional): List containing the first and last timeslice of the plateau region indentified for this correlator.
reweighted
def gamma_method(self, **kwargs)
120    def gamma_method(self, **kwargs):
121        """Apply the gamma method to the content of the Corr."""
122        for item in self.content:
123            if not(item is None):
124                if self.N == 1:
125                    item[0].gamma_method(**kwargs)
126                else:
127                    for i in range(self.N):
128                        for j in range(self.N):
129                            item[i, j].gamma_method(**kwargs)

Apply the gamma method to the content of the Corr.

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

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

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

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

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

Parameters
  • i (int): First index to be picked.
  • j (int): Second index to be picked.
def plottable(self)
185    def plottable(self):
186        """Outputs the correlator in a plotable format.
187
188        Outputs three lists containing the timeslice index, the value on each
189        timeslice and the error on each timeslice.
190        """
191        if self.N != 1:
192            raise Exception("Can only make Corr[N=1] plottable")
193        x_list = [x for x in range(self.T) if not self.content[x] is None]
194        y_list = [y[0].value for y in self.content if y is not None]
195        y_err_list = [y[0].dvalue for y in self.content if y is not None]
196
197        return x_list, y_list, y_err_list

Outputs the correlator in a plotable format.

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

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

Symmetrize the correlator around x0=0.

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

Anti-symmetrize the correlator around x0=0.

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

Symmetrizes the correlator matrices on every timeslice.

def GEVP(self, t0, ts=None, sort='Eigenvalue', **kwargs)
249    def GEVP(self, t0, ts=None, sort="Eigenvalue", **kwargs):
250        r'''Solve the generalized eigenvalue problem on the correlator matrix and returns the corresponding eigenvectors.
251
252        The eigenvectors are sorted according to the descending eigenvalues, the zeroth eigenvector(s) correspond to the
253        largest eigenvalue(s). The eigenvector(s) for the individual states can be accessed via slicing
254        ```python
255        C.GEVP(t0=2)[0]  # Ground state vector(s)
256        C.GEVP(t0=2)[:3]  # Vectors for the lowest three states
257        ```
258
259        Parameters
260        ----------
261        t0 : int
262            The time t0 for the right hand side of the GEVP according to $G(t)v_i=\lambda_i G(t_0)v_i$
263        ts : int
264            fixed time $G(t_s)v_i=\lambda_i G(t_0)v_i$ if sort=None.
265            If sort="Eigenvector" it gives a reference point for the sorting method.
266        sort : string
267            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.
268            - "Eigenvalue": The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice.
269            - "Eigenvector": Use the method described in arXiv:2004.10472 to find the set of v(t) belonging to the state.
270              The reference state is identified by its eigenvalue at $t=t_s$.
271
272        Other Parameters
273        ----------------
274        state : int
275           Returns only the vector(s) for a specified state. The lowest state is zero.
276        '''
277
278        if self.N == 1:
279            raise Exception("GEVP methods only works on correlator matrices and not single correlators.")
280        if ts is not None:
281            if (ts <= t0):
282                raise Exception("ts has to be larger than t0.")
283
284        if "sorted_list" in kwargs:
285            warnings.warn("Argument 'sorted_list' is deprecated, use 'sort' instead.", DeprecationWarning)
286            sort = kwargs.get("sorted_list")
287
288        symmetric_corr = self.matrix_symmetric()
289        if sort is None:
290            if (ts is None):
291                raise Exception("ts is required if sort=None.")
292            if (self.content[t0] is None) or (self.content[ts] is None):
293                raise Exception("Corr not defined at t0/ts.")
294            G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double")
295            for i in range(self.N):
296                for j in range(self.N):
297                    G0[i, j] = symmetric_corr[t0][i, j].value
298                    Gt[i, j] = symmetric_corr[ts][i, j].value
299
300            reordered_vecs = _GEVP_solver(Gt, G0)
301
302        elif sort in ["Eigenvalue", "Eigenvector"]:
303            if sort == "Eigenvalue" and ts is not None:
304                warnings.warn("ts has no effect when sorting by eigenvalue is chosen.", RuntimeWarning)
305            all_vecs = [None] * (t0 + 1)
306            for t in range(t0 + 1, self.T):
307                try:
308                    G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double")
309                    for i in range(self.N):
310                        for j in range(self.N):
311                            G0[i, j] = symmetric_corr[t0][i, j].value
312                            Gt[i, j] = symmetric_corr[t][i, j].value
313
314                    all_vecs.append(_GEVP_solver(Gt, G0))
315                except Exception:
316                    all_vecs.append(None)
317            if sort == "Eigenvector":
318                if (ts is None):
319                    raise Exception("ts is required for the Eigenvector sorting method.")
320                all_vecs = _sort_vectors(all_vecs, ts)
321
322            reordered_vecs = [[v[s] if v is not None else None for v in all_vecs] for s in range(self.N)]
323        else:
324            raise Exception("Unkown value for 'sort'.")
325
326        if "state" in kwargs:
327            return reordered_vecs[kwargs.get("state")]
328        else:
329            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')
331    def Eigenvalue(self, t0, ts=None, state=0, sort="Eigenvalue"):
332        """Determines the eigenvalue of the GEVP by solving and projecting the correlator
333
334        Parameters
335        ----------
336        state : int
337            The state one is interested in ordered by energy. The lowest state is zero.
338
339        All other parameters are identical to the ones of Corr.GEVP.
340        """
341        vec = self.GEVP(t0, ts=ts, sort=sort)[state]
342        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)
344    def Hankel(self, N, periodic=False):
345        """Constructs an NxN Hankel matrix
346
347        C(t) c(t+1) ... c(t+n-1)
348        C(t+1) c(t+2) ... c(t+n)
349        .................
350        C(t+(n-1)) c(t+n) ... c(t+2(n-1))
351
352        Parameters
353        ----------
354        N : int
355            Dimension of the Hankel matrix
356        periodic : bool, optional
357            determines whether the matrix is extended periodically
358        """
359
360        if self.N != 1:
361            raise Exception("Multi-operator Prony not implemented!")
362
363        array = np.empty([N, N], dtype="object")
364        new_content = []
365        for t in range(self.T):
366            new_content.append(array.copy())
367
368        def wrap(i):
369            while i >= self.T:
370                i -= self.T
371            return i
372
373        for t in range(self.T):
374            for i in range(N):
375                for j in range(N):
376                    if periodic:
377                        new_content[t][i, j] = self.content[wrap(t + i + j)][0]
378                    elif (t + i + j) >= self.T:
379                        new_content[t] = None
380                    else:
381                        new_content[t][i, j] = self.content[t + i + j][0]
382
383        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)
385    def roll(self, dt):
386        """Periodically shift the correlator by dt timeslices
387
388        Parameters
389        ----------
390        dt : int
391            number of timeslices
392        """
393        return Corr(list(np.roll(np.array(self.content, dtype=object), dt)))

Periodically shift the correlator by dt timeslices

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

Reverse the time ordering of the Corr

def thin(self, spacing=2, offset=0)
399    def thin(self, spacing=2, offset=0):
400        """Thin out a correlator to suppress correlations
401
402        Parameters
403        ----------
404        spacing : int
405            Keep only every 'spacing'th entry of the correlator
406        offset : int
407            Offset the equal spacing
408        """
409        new_content = []
410        for t in range(self.T):
411            if (offset + t) % spacing != 0:
412                new_content.append(None)
413            else:
414                new_content.append(self.content[t])
415        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)
417    def correlate(self, partner):
418        """Correlate the correlator with another correlator or Obs
419
420        Parameters
421        ----------
422        partner : Obs or Corr
423            partner to correlate the correlator with.
424            Can either be an Obs which is correlated with all entries of the
425            correlator or a Corr of same length.
426        """
427        if self.N != 1:
428            raise Exception("Only one-dimensional correlators can be safely correlated.")
429        new_content = []
430        for x0, t_slice in enumerate(self.content):
431            if _check_for_none(self, t_slice):
432                new_content.append(None)
433            else:
434                if isinstance(partner, Corr):
435                    if _check_for_none(partner, partner.content[x0]):
436                        new_content.append(None)
437                    else:
438                        new_content.append(np.array([correlate(o, partner.content[x0][0]) for o in t_slice]))
439                elif isinstance(partner, Obs):  # Should this include CObs?
440                    new_content.append(np.array([correlate(o, partner) for o in t_slice]))
441                else:
442                    raise Exception("Can only correlate with an Obs or a Corr.")
443
444        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)
446    def reweight(self, weight, **kwargs):
447        """Reweight the correlator.
448
449        Parameters
450        ----------
451        weight : Obs
452            Reweighting factor. An Observable that has to be defined on a superset of the
453            configurations in obs[i].idl for all i.
454        all_configs : bool
455            if True, the reweighted observables are normalized by the average of
456            the reweighting factor on all configurations in weight.idl and not
457            on the configurations in obs[i].idl.
458        """
459        if self.N != 1:
460            raise Exception("Reweighting only implemented for one-dimensional correlators.")
461        new_content = []
462        for t_slice in self.content:
463            if _check_for_none(self, t_slice):
464                new_content.append(None)
465            else:
466                new_content.append(np.array(reweight(weight, t_slice, **kwargs)))
467        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)
469    def T_symmetry(self, partner, parity=+1):
470        """Return the time symmetry average of the correlator and its partner
471
472        Parameters
473        ----------
474        partner : Corr
475            Time symmetry partner of the Corr
476        partity : int
477            Parity quantum number of the correlator, can be +1 or -1
478        """
479        if self.N != 1:
480            raise Exception("T_symmetry only implemented for one-dimensional correlators.")
481        if not isinstance(partner, Corr):
482            raise Exception("T partner has to be a Corr object.")
483        if parity not in [+1, -1]:
484            raise Exception("Parity has to be +1 or -1.")
485        T_partner = parity * partner.reverse()
486
487        t_slices = []
488        test = (self - T_partner)
489        test.gamma_method()
490        for x0, t_slice in enumerate(test.content):
491            if t_slice is not None:
492                if not t_slice[0].is_zero_within_error(5):
493                    t_slices.append(x0)
494        if t_slices:
495            warnings.warn("T symmetry partners do not agree within 5 sigma on time slices " + str(t_slices) + ".", RuntimeWarning)
496
497        return (self + T_partner) / 2

Return the time symmetry average of the correlator and its partner

Parameters
  • partner (Corr): Time symmetry partner of the Corr
  • partity (int): Parity quantum number of the correlator, can be +1 or -1
def deriv(self, variant='symmetric')
499    def deriv(self, variant="symmetric"):
500        """Return the first derivative of the correlator with respect to x0.
501
502        Parameters
503        ----------
504        variant : str
505            decides which definition of the finite differences derivative is used.
506            Available choice: symmetric, forward, backward, improved, default: symmetric
507        """
508        if self.N != 1:
509            raise Exception("deriv only implemented for one-dimensional correlators.")
510        if variant == "symmetric":
511            newcontent = []
512            for t in range(1, self.T - 1):
513                if (self.content[t - 1] is None) or (self.content[t + 1] is None):
514                    newcontent.append(None)
515                else:
516                    newcontent.append(0.5 * (self.content[t + 1] - self.content[t - 1]))
517            if(all([x is None for x in newcontent])):
518                raise Exception('Derivative is undefined at all timeslices')
519            return Corr(newcontent, padding=[1, 1])
520        elif variant == "forward":
521            newcontent = []
522            for t in range(self.T - 1):
523                if (self.content[t] is None) or (self.content[t + 1] is None):
524                    newcontent.append(None)
525                else:
526                    newcontent.append(self.content[t + 1] - self.content[t])
527            if(all([x is None for x in newcontent])):
528                raise Exception("Derivative is undefined at all timeslices")
529            return Corr(newcontent, padding=[0, 1])
530        elif variant == "backward":
531            newcontent = []
532            for t in range(1, self.T):
533                if (self.content[t - 1] is None) or (self.content[t] is None):
534                    newcontent.append(None)
535                else:
536                    newcontent.append(self.content[t] - self.content[t - 1])
537            if(all([x is None for x in newcontent])):
538                raise Exception("Derivative is undefined at all timeslices")
539            return Corr(newcontent, padding=[1, 0])
540        elif variant == "improved":
541            newcontent = []
542            for t in range(2, self.T - 2):
543                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):
544                    newcontent.append(None)
545                else:
546                    newcontent.append((1 / 12) * (self.content[t - 2] - 8 * self.content[t - 1] + 8 * self.content[t + 1] - self.content[t + 2]))
547            if(all([x is None for x in newcontent])):
548                raise Exception('Derivative is undefined at all timeslices')
549            return Corr(newcontent, padding=[2, 2])
550        else:
551            raise Exception("Unknown variant.")

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

Parameters
  • variant (str): decides which definition of the finite differences derivative is used. Available choice: symmetric, forward, backward, improved, default: symmetric
def second_deriv(self, variant='symmetric')
553    def second_deriv(self, variant="symmetric"):
554        """Return the second derivative of the correlator with respect to x0.
555
556        Parameters
557        ----------
558        variant : str
559            decides which definition of the finite differences derivative is used.
560            Available choice: symmetric, improved, default: symmetric
561        """
562        if self.N != 1:
563            raise Exception("second_deriv only implemented for one-dimensional correlators.")
564        if variant == "symmetric":
565            newcontent = []
566            for t in range(1, self.T - 1):
567                if (self.content[t - 1] is None) or (self.content[t + 1] is None):
568                    newcontent.append(None)
569                else:
570                    newcontent.append((self.content[t + 1] - 2 * self.content[t] + self.content[t - 1]))
571            if(all([x is None for x in newcontent])):
572                raise Exception("Derivative is undefined at all timeslices")
573            return Corr(newcontent, padding=[1, 1])
574        elif variant == "improved":
575            newcontent = []
576            for t in range(2, self.T - 2):
577                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):
578                    newcontent.append(None)
579                else:
580                    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]))
581            if(all([x is None for x in newcontent])):
582                raise Exception("Derivative is undefined at all timeslices")
583            return Corr(newcontent, padding=[2, 2])
584        else:
585            raise Exception("Unknown variant.")

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

Parameters
  • variant (str): decides which definition of the finite differences derivative is used. Available choice: symmetric, improved, default: symmetric
def m_eff(self, variant='log', guess=1.0)
587    def m_eff(self, variant='log', guess=1.0):
588        """Returns the effective mass of the correlator as correlator object
589
590        Parameters
591        ----------
592        variant : str
593            log : uses the standard effective mass log(C(t) / C(t+1))
594            cosh, periodic : Use periodicitiy of the correlator by solving C(t) / C(t+1) = cosh(m * (t - T/2)) / cosh(m * (t + 1 - T/2)) for m.
595            sinh : Use anti-periodicitiy of the correlator by solving C(t) / C(t+1) = sinh(m * (t - T/2)) / sinh(m * (t + 1 - T/2)) for m.
596            See, e.g., arXiv:1205.5380
597            arccosh : Uses the explicit form of the symmetrized correlator (not recommended)
598        guess : float
599            guess for the root finder, only relevant for the root variant
600        """
601        if self.N != 1:
602            raise Exception('Correlator must be projected before getting m_eff')
603        if variant == 'log':
604            newcontent = []
605            for t in range(self.T - 1):
606                if ((self.content[t] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0):
607                    newcontent.append(None)
608                else:
609                    newcontent.append(self.content[t] / self.content[t + 1])
610            if(all([x is None for x in newcontent])):
611                raise Exception('m_eff is undefined at all timeslices')
612
613            return np.log(Corr(newcontent, padding=[0, 1]))
614
615        elif variant in ['periodic', 'cosh', 'sinh']:
616            if variant in ['periodic', 'cosh']:
617                func = anp.cosh
618            else:
619                func = anp.sinh
620
621            def root_function(x, d):
622                return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d
623
624            newcontent = []
625            for t in range(self.T - 1):
626                if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t + 1][0].value == 0):
627                    newcontent.append(None)
628                # Fill the two timeslices in the middle of the lattice with their predecessors
629                elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]:
630                    newcontent.append(newcontent[-1])
631                else:
632                    newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess)))
633            if(all([x is None for x in newcontent])):
634                raise Exception('m_eff is undefined at all timeslices')
635
636            return Corr(newcontent, padding=[0, 1])
637
638        elif variant == 'arccosh':
639            newcontent = []
640            for t in range(1, self.T - 1):
641                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):
642                    newcontent.append(None)
643                else:
644                    newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t]))
645            if(all([x is None for x in newcontent])):
646                raise Exception("m_eff is undefined at all timeslices")
647            return np.arccosh(Corr(newcontent, padding=[1, 1]))
648
649        else:
650            raise Exception('Unknown variant.')

Returns the effective mass of the correlator as correlator object

Parameters
  • variant (str): log : uses the standard effective mass log(C(t) / C(t+1)) cosh, periodic : Use periodicitiy of the correlator by solving C(t) / C(t+1) = cosh(m * (t - T/2)) / cosh(m * (t + 1 - T/2)) for m. sinh : Use anti-periodicitiy of the correlator by solving C(t) / C(t+1) = sinh(m * (t - T/2)) / sinh(m * (t + 1 - T/2)) for m. See, e.g., arXiv:1205.5380 arccosh : Uses the explicit form of the symmetrized correlator (not recommended)
  • guess (float): guess for the root finder, only relevant for the root variant
def fit(self, function, fitrange=None, silent=False, **kwargs)
652    def fit(self, function, fitrange=None, silent=False, **kwargs):
653        r'''Fits function to the data
654
655        Parameters
656        ----------
657        function : obj
658            function to fit to the data. See fits.least_squares for details.
659        fitrange : list
660            Two element list containing the timeslices on which the fit is supposed to start and stop.
661            Caution: This range is inclusive as opposed to standard python indexing.
662            `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6.
663            If not specified, self.prange or all timeslices are used.
664        silent : bool
665            Decides whether output is printed to the standard output.
666        '''
667        if self.N != 1:
668            raise Exception("Correlator must be projected before fitting")
669
670        if fitrange is None:
671            if self.prange:
672                fitrange = self.prange
673            else:
674                fitrange = [0, self.T - 1]
675        else:
676            if not isinstance(fitrange, list):
677                raise Exception("fitrange has to be a list with two elements")
678            if len(fitrange) != 2:
679                raise Exception("fitrange has to have exactly two elements [fit_start, fit_stop]")
680
681        xs = [x for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None]
682        ys = [self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None]
683        result = least_squares(xs, ys, function, silent=silent, **kwargs)
684        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)
686    def plateau(self, plateau_range=None, method="fit", auto_gamma=False):
687        """ Extract a plateau value from a Corr object
688
689        Parameters
690        ----------
691        plateau_range : list
692            list with two entries, indicating the first and the last timeslice
693            of the plateau region.
694        method : str
695            method to extract the plateau.
696                'fit' fits a constant to the plateau region
697                'avg', 'average' or 'mean' just average over the given timeslices.
698        auto_gamma : bool
699            apply gamma_method with default parameters to the Corr. Defaults to None
700        """
701        if not plateau_range:
702            if self.prange:
703                plateau_range = self.prange
704            else:
705                raise Exception("no plateau range provided")
706        if self.N != 1:
707            raise Exception("Correlator must be projected before getting a plateau.")
708        if(all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])):
709            raise Exception("plateau is undefined at all timeslices in plateaurange.")
710        if auto_gamma:
711            self.gamma_method()
712        if method == "fit":
713            def const_func(a, t):
714                return a[0]
715            return self.fit(const_func, plateau_range)[0]
716        elif method in ["avg", "average", "mean"]:
717            returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None])
718            return returnvalue
719
720        else:
721            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)
723    def set_prange(self, prange):
724        """Sets the attribute prange of the Corr object."""
725        if not len(prange) == 2:
726            raise Exception("prange must be a list or array with two values")
727        if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))):
728            raise Exception("Start and end point must be integers")
729        if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] < prange[1]):
730            raise Exception("Start and end point must define a range in the interval 0,T")
731
732        self.prange = prange
733        return

Sets the attribute prange of the Corr object.

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

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

Parameters
  • x_range (list): list of two values, determining the range of the x-axis e.g. [4, 8]
  • comp (Corr or list of Corr): Correlator or list of correlators which are plotted for comparison. The tags of these correlators are used as labels if available.
  • logscale (bool): Sets y-axis to logscale
  • plateau (Obs): Plateau value to be visualized in the figure
  • fit_res (Fit_result): Fit_result object to be visualized
  • ylabel (str): Label for the y-axis
  • save (str): path to file in which the figure should be saved
  • auto_gamma (bool): Apply the gamma method with standard parameters to all correlators and plateau values before plotting.
  • hide_sigma (float): Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors.
  • references (list): List of floating point values that are displayed as horizontal lines for reference.
def spaghetti_plot(self, logscale=True)
848    def spaghetti_plot(self, logscale=True):
849        """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations.
850
851        Parameters
852        ----------
853        logscale : bool
854            Determines whether the scale of the y-axis is logarithmic or standard.
855        """
856        if self.N != 1:
857            raise Exception("Correlator needs to be projected first.")
858
859        mc_names = list(set([item for sublist in [o[0].mc_names for o in self.content if o is not None] for item in sublist]))
860        x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None]
861
862        for name in mc_names:
863            data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T
864
865            fig = plt.figure()
866            ax = fig.add_subplot(111)
867            for dat in data:
868                ax.plot(x0_vals, dat, ls='-', marker='')
869
870            if logscale is True:
871                ax.set_yscale('log')
872
873            ax.set_xlabel(r'$x_0 / a$')
874            plt.title(name)
875            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)
877    def dump(self, filename, datatype="json.gz", **kwargs):
878        """Dumps the Corr into a file of chosen type
879        Parameters
880        ----------
881        filename : str
882            Name of the file to be saved.
883        datatype : str
884            Format of the exported file. Supported formats include
885            "json.gz" and "pickle"
886        path : str
887            specifies a custom path for the file (default '.')
888        """
889        if datatype == "json.gz":
890            from .input.json import dump_to_json
891            if 'path' in kwargs:
892                file_name = kwargs.get('path') + '/' + filename
893            else:
894                file_name = filename
895            dump_to_json(self, file_name)
896        elif datatype == "pickle":
897            dump_object(self, filename, **kwargs)
898        else:
899            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)
901    def print(self, print_range=None):
902        print(self.__repr__(print_range))
def sqrt(self)
1066    def sqrt(self):
1067        return self ** 0.5
def log(self)
1069    def log(self):
1070        newcontent = [None if _check_for_none(self, item) else np.log(item) for item in self.content]
1071        return Corr(newcontent, prange=self.prange)
def exp(self)
1073    def exp(self):
1074        newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content]
1075        return Corr(newcontent, prange=self.prange)
def sin(self)
1088    def sin(self):
1089        return self._apply_func_to_corr(np.sin)
def cos(self)
1091    def cos(self):
1092        return self._apply_func_to_corr(np.cos)
def tan(self)
1094    def tan(self):
1095        return self._apply_func_to_corr(np.tan)
def sinh(self)
1097    def sinh(self):
1098        return self._apply_func_to_corr(np.sinh)
def cosh(self)
1100    def cosh(self):
1101        return self._apply_func_to_corr(np.cosh)
def tanh(self)
1103    def tanh(self):
1104        return self._apply_func_to_corr(np.tanh)
def arcsin(self)
1106    def arcsin(self):
1107        return self._apply_func_to_corr(np.arcsin)
def arccos(self)
1109    def arccos(self):
1110        return self._apply_func_to_corr(np.arccos)
def arctan(self)
1112    def arctan(self):
1113        return self._apply_func_to_corr(np.arctan)
def arcsinh(self)
1115    def arcsinh(self):
1116        return self._apply_func_to_corr(np.arcsinh)
def arccosh(self)
1118    def arccosh(self):
1119        return self._apply_func_to_corr(np.arccosh)
def arctanh(self)
1121    def arctanh(self):
1122        return self._apply_func_to_corr(np.arctanh)
real
imag
def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None)
1157    def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
1158        r''' Project large correlation matrix to lowest states
1159
1160        This method can be used to reduce the size of an (N x N) correlation matrix
1161        to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise
1162        is still small.
1163
1164        Parameters
1165        ----------
1166        Ntrunc: int
1167            Rank of the target matrix.
1168        tproj: int
1169            Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method.
1170            The default value is 3.
1171        t0proj: int
1172            Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly
1173            discouraged for O(a) improved theories, since the correctness of the procedure
1174            cannot be granted in this case. The default value is 2.
1175        basematrix : Corr
1176            Correlation matrix that is used to determine the eigenvectors of the
1177            lowest states based on a GEVP. basematrix is taken to be the Corr itself if
1178            is is not specified.
1179
1180        Notes
1181        -----
1182        We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving
1183        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}$
1184        and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the
1185        resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via
1186        $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large
1187        correlation matrix and to remove some noise that is added by irrelevant operators.
1188        This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated
1189        bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$.
1190        '''
1191
1192        if self.N == 1:
1193            raise Exception('Method cannot be applied to one-dimensional correlators.')
1194        if basematrix is None:
1195            basematrix = self
1196        if Ntrunc >= basematrix.N:
1197            raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N))
1198        if basematrix.N != self.N:
1199            raise Exception('basematrix and targetmatrix have to be of the same size.')
1200
1201        evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc]
1202
1203        tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object)
1204        rmat = []
1205        for t in range(basematrix.T):
1206            for i in range(Ntrunc):
1207                for j in range(Ntrunc):
1208                    tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j]
1209            rmat.append(np.copy(tmpmat))
1210
1211        newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)]
1212        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)$.