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, title=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        title : string
 761            Optional title of the figure.
 762        """
 763        if self.N != 1:
 764            raise Exception("Correlator must be projected before plotting")
 765
 766        if auto_gamma:
 767            self.gamma_method()
 768
 769        if x_range is None:
 770            x_range = [0, self.T - 1]
 771
 772        fig = plt.figure()
 773        ax1 = fig.add_subplot(111)
 774
 775        x, y, y_err = self.plottable()
 776        if hide_sigma:
 777            hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
 778        else:
 779            hide_from = None
 780        ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag)
 781        if logscale:
 782            ax1.set_yscale('log')
 783        else:
 784            if y_range is None:
 785                try:
 786                    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)])
 787                    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)])
 788                    ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)])
 789                except Exception:
 790                    pass
 791            else:
 792                ax1.set_ylim(y_range)
 793        if comp:
 794            if isinstance(comp, (Corr, list)):
 795                for corr in comp if isinstance(comp, list) else [comp]:
 796                    if auto_gamma:
 797                        corr.gamma_method()
 798                    x, y, y_err = corr.plottable()
 799                    if hide_sigma:
 800                        hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
 801                    else:
 802                        hide_from = None
 803                    plt.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor'])
 804            else:
 805                raise Exception("'comp' must be a correlator or a list of correlators.")
 806
 807        if plateau:
 808            if isinstance(plateau, Obs):
 809                if auto_gamma:
 810                    plateau.gamma_method()
 811                ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau))
 812                ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-')
 813            else:
 814                raise Exception("'plateau' must be an Obs")
 815
 816        if references:
 817            if isinstance(references, list):
 818                for ref in references:
 819                    ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--')
 820            else:
 821                raise Exception("'references' must be a list of floating pint values.")
 822
 823        if self.prange:
 824            ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',')
 825            ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',')
 826
 827        if fit_res:
 828            x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05)
 829            ax1.plot(x_samples,
 830                     fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples),
 831                     ls='-', marker=',', lw=2)
 832
 833        ax1.set_xlabel(r'$x_0 / a$')
 834        if ylabel:
 835            ax1.set_ylabel(ylabel)
 836        ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5])
 837
 838        handles, labels = ax1.get_legend_handles_labels()
 839        if labels:
 840            ax1.legend()
 841
 842        if title:
 843            plt.title(title)
 844
 845        plt.draw()
 846
 847        if save:
 848            if isinstance(save, str):
 849                fig.savefig(save, bbox_inches='tight')
 850            else:
 851                raise Exception("'save' has to be a string.")
 852
 853    def spaghetti_plot(self, logscale=True):
 854        """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations.
 855
 856        Parameters
 857        ----------
 858        logscale : bool
 859            Determines whether the scale of the y-axis is logarithmic or standard.
 860        """
 861        if self.N != 1:
 862            raise Exception("Correlator needs to be projected first.")
 863
 864        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]))
 865        x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None]
 866
 867        for name in mc_names:
 868            data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T
 869
 870            fig = plt.figure()
 871            ax = fig.add_subplot(111)
 872            for dat in data:
 873                ax.plot(x0_vals, dat, ls='-', marker='')
 874
 875            if logscale is True:
 876                ax.set_yscale('log')
 877
 878            ax.set_xlabel(r'$x_0 / a$')
 879            plt.title(name)
 880            plt.draw()
 881
 882    def dump(self, filename, datatype="json.gz", **kwargs):
 883        """Dumps the Corr into a file of chosen type
 884        Parameters
 885        ----------
 886        filename : str
 887            Name of the file to be saved.
 888        datatype : str
 889            Format of the exported file. Supported formats include
 890            "json.gz" and "pickle"
 891        path : str
 892            specifies a custom path for the file (default '.')
 893        """
 894        if datatype == "json.gz":
 895            from .input.json import dump_to_json
 896            if 'path' in kwargs:
 897                file_name = kwargs.get('path') + '/' + filename
 898            else:
 899                file_name = filename
 900            dump_to_json(self, file_name)
 901        elif datatype == "pickle":
 902            dump_object(self, filename, **kwargs)
 903        else:
 904            raise Exception("Unknown datatype " + str(datatype))
 905
 906    def print(self, print_range=None):
 907        print(self.__repr__(print_range))
 908
 909    def __repr__(self, print_range=None):
 910        if print_range is None:
 911            print_range = [0, None]
 912
 913        content_string = ""
 914        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
 915
 916        if self.tag is not None:
 917            content_string += "Description: " + self.tag + "\n"
 918        if self.N != 1:
 919            return content_string
 920
 921        if print_range[1]:
 922            print_range[1] += 1
 923        content_string += 'x0/a\tCorr(x0/a)\n------------------\n'
 924        for i, sub_corr in enumerate(self.content[print_range[0]:print_range[1]]):
 925            if sub_corr is None:
 926                content_string += str(i + print_range[0]) + '\n'
 927            else:
 928                content_string += str(i + print_range[0])
 929                for element in sub_corr:
 930                    content_string += '\t' + ' ' * int(element >= 0) + str(element)
 931                content_string += '\n'
 932        return content_string
 933
 934    def __str__(self):
 935        return self.__repr__()
 936
 937    # We define the basic operations, that can be performed with correlators.
 938    # While */+- get defined here, they only work for Corr*Obs and not Obs*Corr.
 939    # This is because Obs*Corr checks Obs.__mul__ first and does not catch an exception.
 940    # One could try and tell Obs to check if the y in __mul__ is a Corr and
 941
 942    def __add__(self, y):
 943        if isinstance(y, Corr):
 944            if ((self.N != y.N) or (self.T != y.T)):
 945                raise Exception("Addition of Corrs with different shape")
 946            newcontent = []
 947            for t in range(self.T):
 948                if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]):
 949                    newcontent.append(None)
 950                else:
 951                    newcontent.append(self.content[t] + y.content[t])
 952            return Corr(newcontent)
 953
 954        elif isinstance(y, (Obs, int, float, CObs)):
 955            newcontent = []
 956            for t in range(self.T):
 957                if _check_for_none(self, self.content[t]):
 958                    newcontent.append(None)
 959                else:
 960                    newcontent.append(self.content[t] + y)
 961            return Corr(newcontent, prange=self.prange)
 962        elif isinstance(y, np.ndarray):
 963            if y.shape == (self.T,):
 964                return Corr(list((np.array(self.content).T + y).T))
 965            else:
 966                raise ValueError("operands could not be broadcast together")
 967        else:
 968            raise TypeError("Corr + wrong type")
 969
 970    def __mul__(self, y):
 971        if isinstance(y, Corr):
 972            if not((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T):
 973                raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T")
 974            newcontent = []
 975            for t in range(self.T):
 976                if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]):
 977                    newcontent.append(None)
 978                else:
 979                    newcontent.append(self.content[t] * y.content[t])
 980            return Corr(newcontent)
 981
 982        elif isinstance(y, (Obs, int, float, CObs)):
 983            newcontent = []
 984            for t in range(self.T):
 985                if _check_for_none(self, self.content[t]):
 986                    newcontent.append(None)
 987                else:
 988                    newcontent.append(self.content[t] * y)
 989            return Corr(newcontent, prange=self.prange)
 990        elif isinstance(y, np.ndarray):
 991            if y.shape == (self.T,):
 992                return Corr(list((np.array(self.content).T * y).T))
 993            else:
 994                raise ValueError("operands could not be broadcast together")
 995        else:
 996            raise TypeError("Corr * wrong type")
 997
 998    def __truediv__(self, y):
 999        if isinstance(y, Corr):
1000            if not((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T):
1001                raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T")
1002            newcontent = []
1003            for t in range(self.T):
1004                if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]):
1005                    newcontent.append(None)
1006                else:
1007                    newcontent.append(self.content[t] / y.content[t])
1008            for t in range(self.T):
1009                if _check_for_none(self, newcontent[t]):
1010                    continue
1011                if np.isnan(np.sum(newcontent[t]).value):
1012                    newcontent[t] = None
1013
1014            if all([item is None for item in newcontent]):
1015                raise Exception("Division returns completely undefined correlator")
1016            return Corr(newcontent)
1017
1018        elif isinstance(y, (Obs, CObs)):
1019            if isinstance(y, Obs):
1020                if y.value == 0:
1021                    raise Exception('Division by zero will return undefined correlator')
1022            if isinstance(y, CObs):
1023                if y.is_zero():
1024                    raise Exception('Division by zero will return undefined correlator')
1025
1026            newcontent = []
1027            for t in range(self.T):
1028                if _check_for_none(self, self.content[t]):
1029                    newcontent.append(None)
1030                else:
1031                    newcontent.append(self.content[t] / y)
1032            return Corr(newcontent, prange=self.prange)
1033
1034        elif isinstance(y, (int, float)):
1035            if y == 0:
1036                raise Exception('Division by zero will return undefined correlator')
1037            newcontent = []
1038            for t in range(self.T):
1039                if _check_for_none(self, self.content[t]):
1040                    newcontent.append(None)
1041                else:
1042                    newcontent.append(self.content[t] / y)
1043            return Corr(newcontent, prange=self.prange)
1044        elif isinstance(y, np.ndarray):
1045            if y.shape == (self.T,):
1046                return Corr(list((np.array(self.content).T / y).T))
1047            else:
1048                raise ValueError("operands could not be broadcast together")
1049        else:
1050            raise TypeError('Corr / wrong type')
1051
1052    def __neg__(self):
1053        newcontent = [None if _check_for_none(self, item) else -1. * item for item in self.content]
1054        return Corr(newcontent, prange=self.prange)
1055
1056    def __sub__(self, y):
1057        return self + (-y)
1058
1059    def __pow__(self, y):
1060        if isinstance(y, (Obs, int, float, CObs)):
1061            newcontent = [None if _check_for_none(self, item) else item**y for item in self.content]
1062            return Corr(newcontent, prange=self.prange)
1063        else:
1064            raise TypeError('Type of exponent not supported')
1065
1066    def __abs__(self):
1067        newcontent = [None if _check_for_none(self, item) else np.abs(item) for item in self.content]
1068        return Corr(newcontent, prange=self.prange)
1069
1070    # The numpy functions:
1071    def sqrt(self):
1072        return self ** 0.5
1073
1074    def log(self):
1075        newcontent = [None if _check_for_none(self, item) else np.log(item) for item in self.content]
1076        return Corr(newcontent, prange=self.prange)
1077
1078    def exp(self):
1079        newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content]
1080        return Corr(newcontent, prange=self.prange)
1081
1082    def _apply_func_to_corr(self, func):
1083        newcontent = [None if _check_for_none(self, item) else func(item) for item in self.content]
1084        for t in range(self.T):
1085            if _check_for_none(self, newcontent[t]):
1086                continue
1087            if np.isnan(np.sum(newcontent[t]).value):
1088                newcontent[t] = None
1089        if all([item is None for item in newcontent]):
1090            raise Exception('Operation returns undefined correlator')
1091        return Corr(newcontent)
1092
1093    def sin(self):
1094        return self._apply_func_to_corr(np.sin)
1095
1096    def cos(self):
1097        return self._apply_func_to_corr(np.cos)
1098
1099    def tan(self):
1100        return self._apply_func_to_corr(np.tan)
1101
1102    def sinh(self):
1103        return self._apply_func_to_corr(np.sinh)
1104
1105    def cosh(self):
1106        return self._apply_func_to_corr(np.cosh)
1107
1108    def tanh(self):
1109        return self._apply_func_to_corr(np.tanh)
1110
1111    def arcsin(self):
1112        return self._apply_func_to_corr(np.arcsin)
1113
1114    def arccos(self):
1115        return self._apply_func_to_corr(np.arccos)
1116
1117    def arctan(self):
1118        return self._apply_func_to_corr(np.arctan)
1119
1120    def arcsinh(self):
1121        return self._apply_func_to_corr(np.arcsinh)
1122
1123    def arccosh(self):
1124        return self._apply_func_to_corr(np.arccosh)
1125
1126    def arctanh(self):
1127        return self._apply_func_to_corr(np.arctanh)
1128
1129    # Right hand side operations (require tweak in main module to work)
1130    def __radd__(self, y):
1131        return self + y
1132
1133    def __rsub__(self, y):
1134        return -self + y
1135
1136    def __rmul__(self, y):
1137        return self * y
1138
1139    def __rtruediv__(self, y):
1140        return (self / y) ** (-1)
1141
1142    @property
1143    def real(self):
1144        def return_real(obs_OR_cobs):
1145            if isinstance(obs_OR_cobs, CObs):
1146                return obs_OR_cobs.real
1147            else:
1148                return obs_OR_cobs
1149
1150        return self._apply_func_to_corr(return_real)
1151
1152    @property
1153    def imag(self):
1154        def return_imag(obs_OR_cobs):
1155            if isinstance(obs_OR_cobs, CObs):
1156                return obs_OR_cobs.imag
1157            else:
1158                return obs_OR_cobs * 0  # So it stays the right type
1159
1160        return self._apply_func_to_corr(return_imag)
1161
1162    def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
1163        r''' Project large correlation matrix to lowest states
1164
1165        This method can be used to reduce the size of an (N x N) correlation matrix
1166        to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise
1167        is still small.
1168
1169        Parameters
1170        ----------
1171        Ntrunc: int
1172            Rank of the target matrix.
1173        tproj: int
1174            Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method.
1175            The default value is 3.
1176        t0proj: int
1177            Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly
1178            discouraged for O(a) improved theories, since the correctness of the procedure
1179            cannot be granted in this case. The default value is 2.
1180        basematrix : Corr
1181            Correlation matrix that is used to determine the eigenvectors of the
1182            lowest states based on a GEVP. basematrix is taken to be the Corr itself if
1183            is is not specified.
1184
1185        Notes
1186        -----
1187        We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving
1188        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}$
1189        and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the
1190        resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via
1191        $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large
1192        correlation matrix and to remove some noise that is added by irrelevant operators.
1193        This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated
1194        bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$.
1195        '''
1196
1197        if self.N == 1:
1198            raise Exception('Method cannot be applied to one-dimensional correlators.')
1199        if basematrix is None:
1200            basematrix = self
1201        if Ntrunc >= basematrix.N:
1202            raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N))
1203        if basematrix.N != self.N:
1204            raise Exception('basematrix and targetmatrix have to be of the same size.')
1205
1206        evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc]
1207
1208        tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object)
1209        rmat = []
1210        for t in range(basematrix.T):
1211            for i in range(Ntrunc):
1212                for j in range(Ntrunc):
1213                    tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j]
1214            rmat.append(np.copy(tmpmat))
1215
1216        newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)]
1217        return Corr(newcontent)
1218
1219
1220def _sort_vectors(vec_set, ts):
1221    """Helper function used to find a set of Eigenvectors consistent over all timeslices"""
1222    reference_sorting = np.array(vec_set[ts])
1223    N = reference_sorting.shape[0]
1224    sorted_vec_set = []
1225    for t in range(len(vec_set)):
1226        if vec_set[t] is None:
1227            sorted_vec_set.append(None)
1228        elif not t == ts:
1229            perms = [list(o) for o in permutations([i for i in range(N)], N)]
1230            best_score = 0
1231            for perm in perms:
1232                current_score = 1
1233                for k in range(N):
1234                    new_sorting = reference_sorting.copy()
1235                    new_sorting[perm[k], :] = vec_set[t][k]
1236                    current_score *= abs(np.linalg.det(new_sorting))
1237                if current_score > best_score:
1238                    best_score = current_score
1239                    best_perm = perm
1240            sorted_vec_set.append([vec_set[t][k] for k in best_perm])
1241        else:
1242            sorted_vec_set.append(vec_set[t])
1243
1244    return sorted_vec_set
1245
1246
1247def _check_for_none(corr, entry):
1248    """Checks if entry for correlator corr is None"""
1249    return len(list(filter(None, np.asarray(entry).flatten()))) < corr.N ** 2
1250
1251
1252def _GEVP_solver(Gt, G0):
1253    """Helper function for solving the GEVP and sorting the eigenvectors.
1254
1255    The helper function assumes that both provided matrices are symmetric and
1256    only processes the lower triangular part of both matrices. In case the matrices
1257    are not symmetric the upper triangular parts are effectively discarded."""
1258    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, title=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        title : string
 762            Optional title of the figure.
 763        """
 764        if self.N != 1:
 765            raise Exception("Correlator must be projected before plotting")
 766
 767        if auto_gamma:
 768            self.gamma_method()
 769
 770        if x_range is None:
 771            x_range = [0, self.T - 1]
 772
 773        fig = plt.figure()
 774        ax1 = fig.add_subplot(111)
 775
 776        x, y, y_err = self.plottable()
 777        if hide_sigma:
 778            hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
 779        else:
 780            hide_from = None
 781        ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag)
 782        if logscale:
 783            ax1.set_yscale('log')
 784        else:
 785            if y_range is None:
 786                try:
 787                    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)])
 788                    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)])
 789                    ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)])
 790                except Exception:
 791                    pass
 792            else:
 793                ax1.set_ylim(y_range)
 794        if comp:
 795            if isinstance(comp, (Corr, list)):
 796                for corr in comp if isinstance(comp, list) else [comp]:
 797                    if auto_gamma:
 798                        corr.gamma_method()
 799                    x, y, y_err = corr.plottable()
 800                    if hide_sigma:
 801                        hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
 802                    else:
 803                        hide_from = None
 804                    plt.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor'])
 805            else:
 806                raise Exception("'comp' must be a correlator or a list of correlators.")
 807
 808        if plateau:
 809            if isinstance(plateau, Obs):
 810                if auto_gamma:
 811                    plateau.gamma_method()
 812                ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau))
 813                ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-')
 814            else:
 815                raise Exception("'plateau' must be an Obs")
 816
 817        if references:
 818            if isinstance(references, list):
 819                for ref in references:
 820                    ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--')
 821            else:
 822                raise Exception("'references' must be a list of floating pint values.")
 823
 824        if self.prange:
 825            ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',')
 826            ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',')
 827
 828        if fit_res:
 829            x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05)
 830            ax1.plot(x_samples,
 831                     fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples),
 832                     ls='-', marker=',', lw=2)
 833
 834        ax1.set_xlabel(r'$x_0 / a$')
 835        if ylabel:
 836            ax1.set_ylabel(ylabel)
 837        ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5])
 838
 839        handles, labels = ax1.get_legend_handles_labels()
 840        if labels:
 841            ax1.legend()
 842
 843        if title:
 844            plt.title(title)
 845
 846        plt.draw()
 847
 848        if save:
 849            if isinstance(save, str):
 850                fig.savefig(save, bbox_inches='tight')
 851            else:
 852                raise Exception("'save' has to be a string.")
 853
 854    def spaghetti_plot(self, logscale=True):
 855        """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations.
 856
 857        Parameters
 858        ----------
 859        logscale : bool
 860            Determines whether the scale of the y-axis is logarithmic or standard.
 861        """
 862        if self.N != 1:
 863            raise Exception("Correlator needs to be projected first.")
 864
 865        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]))
 866        x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None]
 867
 868        for name in mc_names:
 869            data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T
 870
 871            fig = plt.figure()
 872            ax = fig.add_subplot(111)
 873            for dat in data:
 874                ax.plot(x0_vals, dat, ls='-', marker='')
 875
 876            if logscale is True:
 877                ax.set_yscale('log')
 878
 879            ax.set_xlabel(r'$x_0 / a$')
 880            plt.title(name)
 881            plt.draw()
 882
 883    def dump(self, filename, datatype="json.gz", **kwargs):
 884        """Dumps the Corr into a file of chosen type
 885        Parameters
 886        ----------
 887        filename : str
 888            Name of the file to be saved.
 889        datatype : str
 890            Format of the exported file. Supported formats include
 891            "json.gz" and "pickle"
 892        path : str
 893            specifies a custom path for the file (default '.')
 894        """
 895        if datatype == "json.gz":
 896            from .input.json import dump_to_json
 897            if 'path' in kwargs:
 898                file_name = kwargs.get('path') + '/' + filename
 899            else:
 900                file_name = filename
 901            dump_to_json(self, file_name)
 902        elif datatype == "pickle":
 903            dump_object(self, filename, **kwargs)
 904        else:
 905            raise Exception("Unknown datatype " + str(datatype))
 906
 907    def print(self, print_range=None):
 908        print(self.__repr__(print_range))
 909
 910    def __repr__(self, print_range=None):
 911        if print_range is None:
 912            print_range = [0, None]
 913
 914        content_string = ""
 915        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
 916
 917        if self.tag is not None:
 918            content_string += "Description: " + self.tag + "\n"
 919        if self.N != 1:
 920            return content_string
 921
 922        if print_range[1]:
 923            print_range[1] += 1
 924        content_string += 'x0/a\tCorr(x0/a)\n------------------\n'
 925        for i, sub_corr in enumerate(self.content[print_range[0]:print_range[1]]):
 926            if sub_corr is None:
 927                content_string += str(i + print_range[0]) + '\n'
 928            else:
 929                content_string += str(i + print_range[0])
 930                for element in sub_corr:
 931                    content_string += '\t' + ' ' * int(element >= 0) + str(element)
 932                content_string += '\n'
 933        return content_string
 934
 935    def __str__(self):
 936        return self.__repr__()
 937
 938    # We define the basic operations, that can be performed with correlators.
 939    # While */+- get defined here, they only work for Corr*Obs and not Obs*Corr.
 940    # This is because Obs*Corr checks Obs.__mul__ first and does not catch an exception.
 941    # One could try and tell Obs to check if the y in __mul__ is a Corr and
 942
 943    def __add__(self, y):
 944        if isinstance(y, Corr):
 945            if ((self.N != y.N) or (self.T != y.T)):
 946                raise Exception("Addition of Corrs with different shape")
 947            newcontent = []
 948            for t in range(self.T):
 949                if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]):
 950                    newcontent.append(None)
 951                else:
 952                    newcontent.append(self.content[t] + y.content[t])
 953            return Corr(newcontent)
 954
 955        elif isinstance(y, (Obs, int, float, CObs)):
 956            newcontent = []
 957            for t in range(self.T):
 958                if _check_for_none(self, self.content[t]):
 959                    newcontent.append(None)
 960                else:
 961                    newcontent.append(self.content[t] + y)
 962            return Corr(newcontent, prange=self.prange)
 963        elif isinstance(y, np.ndarray):
 964            if y.shape == (self.T,):
 965                return Corr(list((np.array(self.content).T + y).T))
 966            else:
 967                raise ValueError("operands could not be broadcast together")
 968        else:
 969            raise TypeError("Corr + wrong type")
 970
 971    def __mul__(self, y):
 972        if isinstance(y, Corr):
 973            if not((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T):
 974                raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T")
 975            newcontent = []
 976            for t in range(self.T):
 977                if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]):
 978                    newcontent.append(None)
 979                else:
 980                    newcontent.append(self.content[t] * y.content[t])
 981            return Corr(newcontent)
 982
 983        elif isinstance(y, (Obs, int, float, CObs)):
 984            newcontent = []
 985            for t in range(self.T):
 986                if _check_for_none(self, self.content[t]):
 987                    newcontent.append(None)
 988                else:
 989                    newcontent.append(self.content[t] * y)
 990            return Corr(newcontent, prange=self.prange)
 991        elif isinstance(y, np.ndarray):
 992            if y.shape == (self.T,):
 993                return Corr(list((np.array(self.content).T * y).T))
 994            else:
 995                raise ValueError("operands could not be broadcast together")
 996        else:
 997            raise TypeError("Corr * wrong type")
 998
 999    def __truediv__(self, y):
1000        if isinstance(y, Corr):
1001            if not((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T):
1002                raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T")
1003            newcontent = []
1004            for t in range(self.T):
1005                if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]):
1006                    newcontent.append(None)
1007                else:
1008                    newcontent.append(self.content[t] / y.content[t])
1009            for t in range(self.T):
1010                if _check_for_none(self, newcontent[t]):
1011                    continue
1012                if np.isnan(np.sum(newcontent[t]).value):
1013                    newcontent[t] = None
1014
1015            if all([item is None for item in newcontent]):
1016                raise Exception("Division returns completely undefined correlator")
1017            return Corr(newcontent)
1018
1019        elif isinstance(y, (Obs, CObs)):
1020            if isinstance(y, Obs):
1021                if y.value == 0:
1022                    raise Exception('Division by zero will return undefined correlator')
1023            if isinstance(y, CObs):
1024                if y.is_zero():
1025                    raise Exception('Division by zero will return undefined correlator')
1026
1027            newcontent = []
1028            for t in range(self.T):
1029                if _check_for_none(self, self.content[t]):
1030                    newcontent.append(None)
1031                else:
1032                    newcontent.append(self.content[t] / y)
1033            return Corr(newcontent, prange=self.prange)
1034
1035        elif isinstance(y, (int, float)):
1036            if y == 0:
1037                raise Exception('Division by zero will return undefined correlator')
1038            newcontent = []
1039            for t in range(self.T):
1040                if _check_for_none(self, self.content[t]):
1041                    newcontent.append(None)
1042                else:
1043                    newcontent.append(self.content[t] / y)
1044            return Corr(newcontent, prange=self.prange)
1045        elif isinstance(y, np.ndarray):
1046            if y.shape == (self.T,):
1047                return Corr(list((np.array(self.content).T / y).T))
1048            else:
1049                raise ValueError("operands could not be broadcast together")
1050        else:
1051            raise TypeError('Corr / wrong type')
1052
1053    def __neg__(self):
1054        newcontent = [None if _check_for_none(self, item) else -1. * item for item in self.content]
1055        return Corr(newcontent, prange=self.prange)
1056
1057    def __sub__(self, y):
1058        return self + (-y)
1059
1060    def __pow__(self, y):
1061        if isinstance(y, (Obs, int, float, CObs)):
1062            newcontent = [None if _check_for_none(self, item) else item**y for item in self.content]
1063            return Corr(newcontent, prange=self.prange)
1064        else:
1065            raise TypeError('Type of exponent not supported')
1066
1067    def __abs__(self):
1068        newcontent = [None if _check_for_none(self, item) else np.abs(item) for item in self.content]
1069        return Corr(newcontent, prange=self.prange)
1070
1071    # The numpy functions:
1072    def sqrt(self):
1073        return self ** 0.5
1074
1075    def log(self):
1076        newcontent = [None if _check_for_none(self, item) else np.log(item) for item in self.content]
1077        return Corr(newcontent, prange=self.prange)
1078
1079    def exp(self):
1080        newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content]
1081        return Corr(newcontent, prange=self.prange)
1082
1083    def _apply_func_to_corr(self, func):
1084        newcontent = [None if _check_for_none(self, item) else func(item) for item in self.content]
1085        for t in range(self.T):
1086            if _check_for_none(self, newcontent[t]):
1087                continue
1088            if np.isnan(np.sum(newcontent[t]).value):
1089                newcontent[t] = None
1090        if all([item is None for item in newcontent]):
1091            raise Exception('Operation returns undefined correlator')
1092        return Corr(newcontent)
1093
1094    def sin(self):
1095        return self._apply_func_to_corr(np.sin)
1096
1097    def cos(self):
1098        return self._apply_func_to_corr(np.cos)
1099
1100    def tan(self):
1101        return self._apply_func_to_corr(np.tan)
1102
1103    def sinh(self):
1104        return self._apply_func_to_corr(np.sinh)
1105
1106    def cosh(self):
1107        return self._apply_func_to_corr(np.cosh)
1108
1109    def tanh(self):
1110        return self._apply_func_to_corr(np.tanh)
1111
1112    def arcsin(self):
1113        return self._apply_func_to_corr(np.arcsin)
1114
1115    def arccos(self):
1116        return self._apply_func_to_corr(np.arccos)
1117
1118    def arctan(self):
1119        return self._apply_func_to_corr(np.arctan)
1120
1121    def arcsinh(self):
1122        return self._apply_func_to_corr(np.arcsinh)
1123
1124    def arccosh(self):
1125        return self._apply_func_to_corr(np.arccosh)
1126
1127    def arctanh(self):
1128        return self._apply_func_to_corr(np.arctanh)
1129
1130    # Right hand side operations (require tweak in main module to work)
1131    def __radd__(self, y):
1132        return self + y
1133
1134    def __rsub__(self, y):
1135        return -self + y
1136
1137    def __rmul__(self, y):
1138        return self * y
1139
1140    def __rtruediv__(self, y):
1141        return (self / y) ** (-1)
1142
1143    @property
1144    def real(self):
1145        def return_real(obs_OR_cobs):
1146            if isinstance(obs_OR_cobs, CObs):
1147                return obs_OR_cobs.real
1148            else:
1149                return obs_OR_cobs
1150
1151        return self._apply_func_to_corr(return_real)
1152
1153    @property
1154    def imag(self):
1155        def return_imag(obs_OR_cobs):
1156            if isinstance(obs_OR_cobs, CObs):
1157                return obs_OR_cobs.imag
1158            else:
1159                return obs_OR_cobs * 0  # So it stays the right type
1160
1161        return self._apply_func_to_corr(return_imag)
1162
1163    def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
1164        r''' Project large correlation matrix to lowest states
1165
1166        This method can be used to reduce the size of an (N x N) correlation matrix
1167        to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise
1168        is still small.
1169
1170        Parameters
1171        ----------
1172        Ntrunc: int
1173            Rank of the target matrix.
1174        tproj: int
1175            Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method.
1176            The default value is 3.
1177        t0proj: int
1178            Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly
1179            discouraged for O(a) improved theories, since the correctness of the procedure
1180            cannot be granted in this case. The default value is 2.
1181        basematrix : Corr
1182            Correlation matrix that is used to determine the eigenvectors of the
1183            lowest states based on a GEVP. basematrix is taken to be the Corr itself if
1184            is is not specified.
1185
1186        Notes
1187        -----
1188        We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving
1189        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}$
1190        and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the
1191        resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via
1192        $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large
1193        correlation matrix and to remove some noise that is added by irrelevant operators.
1194        This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated
1195        bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$.
1196        '''
1197
1198        if self.N == 1:
1199            raise Exception('Method cannot be applied to one-dimensional correlators.')
1200        if basematrix is None:
1201            basematrix = self
1202        if Ntrunc >= basematrix.N:
1203            raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N))
1204        if basematrix.N != self.N:
1205            raise Exception('basematrix and targetmatrix have to be of the same size.')
1206
1207        evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc]
1208
1209        tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object)
1210        rmat = []
1211        for t in range(basematrix.T):
1212            for i in range(Ntrunc):
1213                for j in range(Ntrunc):
1214                    tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j]
1215            rmat.append(np.copy(tmpmat))
1216
1217        newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)]
1218        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, title=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, title=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        title : string
762            Optional title of the figure.
763        """
764        if self.N != 1:
765            raise Exception("Correlator must be projected before plotting")
766
767        if auto_gamma:
768            self.gamma_method()
769
770        if x_range is None:
771            x_range = [0, self.T - 1]
772
773        fig = plt.figure()
774        ax1 = fig.add_subplot(111)
775
776        x, y, y_err = self.plottable()
777        if hide_sigma:
778            hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
779        else:
780            hide_from = None
781        ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag)
782        if logscale:
783            ax1.set_yscale('log')
784        else:
785            if y_range is None:
786                try:
787                    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)])
788                    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)])
789                    ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)])
790                except Exception:
791                    pass
792            else:
793                ax1.set_ylim(y_range)
794        if comp:
795            if isinstance(comp, (Corr, list)):
796                for corr in comp if isinstance(comp, list) else [comp]:
797                    if auto_gamma:
798                        corr.gamma_method()
799                    x, y, y_err = corr.plottable()
800                    if hide_sigma:
801                        hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
802                    else:
803                        hide_from = None
804                    plt.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor'])
805            else:
806                raise Exception("'comp' must be a correlator or a list of correlators.")
807
808        if plateau:
809            if isinstance(plateau, Obs):
810                if auto_gamma:
811                    plateau.gamma_method()
812                ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau))
813                ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-')
814            else:
815                raise Exception("'plateau' must be an Obs")
816
817        if references:
818            if isinstance(references, list):
819                for ref in references:
820                    ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--')
821            else:
822                raise Exception("'references' must be a list of floating pint values.")
823
824        if self.prange:
825            ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',')
826            ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',')
827
828        if fit_res:
829            x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05)
830            ax1.plot(x_samples,
831                     fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples),
832                     ls='-', marker=',', lw=2)
833
834        ax1.set_xlabel(r'$x_0 / a$')
835        if ylabel:
836            ax1.set_ylabel(ylabel)
837        ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5])
838
839        handles, labels = ax1.get_legend_handles_labels()
840        if labels:
841            ax1.legend()
842
843        if title:
844            plt.title(title)
845
846        plt.draw()
847
848        if save:
849            if isinstance(save, str):
850                fig.savefig(save, bbox_inches='tight')
851            else:
852                raise Exception("'save' has to be a string.")

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

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