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

Checks whether a correlator matrices is symmetric on every timeslice.

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

Symmetrizes the correlator matrices on every timeslice.

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

Reverse the time ordering of the Corr

def thin(self, spacing=2, offset=0):
412    def thin(self, spacing=2, offset=0):
413        """Thin out a correlator to suppress correlations
414
415        Parameters
416        ----------
417        spacing : int
418            Keep only every 'spacing'th entry of the correlator
419        offset : int
420            Offset the equal spacing
421        """
422        new_content = []
423        for t in range(self.T):
424            if (offset + t) % spacing != 0:
425                new_content.append(None)
426            else:
427                new_content.append(self.content[t])
428        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):
430    def correlate(self, partner):
431        """Correlate the correlator with another correlator or Obs
432
433        Parameters
434        ----------
435        partner : Obs or Corr
436            partner to correlate the correlator with.
437            Can either be an Obs which is correlated with all entries of the
438            correlator or a Corr of same length.
439        """
440        if self.N != 1:
441            raise Exception("Only one-dimensional correlators can be safely correlated.")
442        new_content = []
443        for x0, t_slice in enumerate(self.content):
444            if _check_for_none(self, t_slice):
445                new_content.append(None)
446            else:
447                if isinstance(partner, Corr):
448                    if _check_for_none(partner, partner.content[x0]):
449                        new_content.append(None)
450                    else:
451                        new_content.append(np.array([correlate(o, partner.content[x0][0]) for o in t_slice]))
452                elif isinstance(partner, Obs):  # Should this include CObs?
453                    new_content.append(np.array([correlate(o, partner) for o in t_slice]))
454                else:
455                    raise Exception("Can only correlate with an Obs or a Corr.")
456
457        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):
459    def reweight(self, weight, **kwargs):
460        """Reweight the correlator.
461
462        Parameters
463        ----------
464        weight : Obs
465            Reweighting factor. An Observable that has to be defined on a superset of the
466            configurations in obs[i].idl for all i.
467        all_configs : bool
468            if True, the reweighted observables are normalized by the average of
469            the reweighting factor on all configurations in weight.idl and not
470            on the configurations in obs[i].idl.
471        """
472        if self.N != 1:
473            raise Exception("Reweighting only implemented for one-dimensional correlators.")
474        new_content = []
475        for t_slice in self.content:
476            if _check_for_none(self, t_slice):
477                new_content.append(None)
478            else:
479                new_content.append(np.array(reweight(weight, t_slice, **kwargs)))
480        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):
482    def T_symmetry(self, partner, parity=+1):
483        """Return the time symmetry average of the correlator and its partner
484
485        Parameters
486        ----------
487        partner : Corr
488            Time symmetry partner of the Corr
489        partity : int
490            Parity quantum number of the correlator, can be +1 or -1
491        """
492        if self.N != 1:
493            raise Exception("T_symmetry only implemented for one-dimensional correlators.")
494        if not isinstance(partner, Corr):
495            raise Exception("T partner has to be a Corr object.")
496        if parity not in [+1, -1]:
497            raise Exception("Parity has to be +1 or -1.")
498        T_partner = parity * partner.reverse()
499
500        t_slices = []
501        test = (self - T_partner)
502        test.gamma_method()
503        for x0, t_slice in enumerate(test.content):
504            if t_slice is not None:
505                if not t_slice[0].is_zero_within_error(5):
506                    t_slices.append(x0)
507        if t_slices:
508            warnings.warn("T symmetry partners do not agree within 5 sigma on time slices " + str(t_slices) + ".", RuntimeWarning)
509
510        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'):
512    def deriv(self, variant="symmetric"):
513        """Return the first derivative of the correlator with respect to x0.
514
515        Parameters
516        ----------
517        variant : str
518            decides which definition of the finite differences derivative is used.
519            Available choice: symmetric, forward, backward, improved, log, default: symmetric
520        """
521        if self.N != 1:
522            raise Exception("deriv only implemented for one-dimensional correlators.")
523        if variant == "symmetric":
524            newcontent = []
525            for t in range(1, self.T - 1):
526                if (self.content[t - 1] is None) or (self.content[t + 1] is None):
527                    newcontent.append(None)
528                else:
529                    newcontent.append(0.5 * (self.content[t + 1] - self.content[t - 1]))
530            if (all([x is None for x in newcontent])):
531                raise Exception('Derivative is undefined at all timeslices')
532            return Corr(newcontent, padding=[1, 1])
533        elif variant == "forward":
534            newcontent = []
535            for t in range(self.T - 1):
536                if (self.content[t] is None) or (self.content[t + 1] is None):
537                    newcontent.append(None)
538                else:
539                    newcontent.append(self.content[t + 1] - self.content[t])
540            if (all([x is None for x in newcontent])):
541                raise Exception("Derivative is undefined at all timeslices")
542            return Corr(newcontent, padding=[0, 1])
543        elif variant == "backward":
544            newcontent = []
545            for t in range(1, self.T):
546                if (self.content[t - 1] is None) or (self.content[t] is None):
547                    newcontent.append(None)
548                else:
549                    newcontent.append(self.content[t] - self.content[t - 1])
550            if (all([x is None for x in newcontent])):
551                raise Exception("Derivative is undefined at all timeslices")
552            return Corr(newcontent, padding=[1, 0])
553        elif variant == "improved":
554            newcontent = []
555            for t in range(2, self.T - 2):
556                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):
557                    newcontent.append(None)
558                else:
559                    newcontent.append((1 / 12) * (self.content[t - 2] - 8 * self.content[t - 1] + 8 * self.content[t + 1] - self.content[t + 2]))
560            if (all([x is None for x in newcontent])):
561                raise Exception('Derivative is undefined at all timeslices')
562            return Corr(newcontent, padding=[2, 2])
563        elif variant == 'log':
564            newcontent = []
565            for t in range(self.T):
566                if (self.content[t] is None) or (self.content[t] <= 0):
567                    newcontent.append(None)
568                else:
569                    newcontent.append(np.log(self.content[t]))
570            if (all([x is None for x in newcontent])):
571                raise Exception("Log is undefined at all timeslices")
572            logcorr = Corr(newcontent)
573            return self * logcorr.deriv('symmetric')
574        else:
575            raise Exception("Unknown variant.")

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

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

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

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

Returns the effective mass of the correlator as correlator object

Parameters
  • variant (str): log : uses the standard effective mass log(C(t) / C(t+1)) cosh, periodic : Use periodicitiy of the correlator by solving C(t) / C(t+1) = cosh(m * (t - T/2)) / cosh(m * (t + 1 - T/2)) for m. sinh : Use anti-periodicitiy of the correlator by solving C(t) / C(t+1) = sinh(m * (t - T/2)) / sinh(m * (t + 1 - T/2)) for m. See, e.g., arXiv:1205.5380 arccosh : Uses the explicit form of the symmetrized correlator (not recommended) logsym: uses the symmetric effective mass log(C(t-1) / C(t+1))/2
  • guess (float): guess for the root finder, only relevant for the root variant
def fit(self, function, fitrange=None, silent=False, **kwargs):
706    def fit(self, function, fitrange=None, silent=False, **kwargs):
707        r'''Fits function to the data
708
709        Parameters
710        ----------
711        function : obj
712            function to fit to the data. See fits.least_squares for details.
713        fitrange : list
714            Two element list containing the timeslices on which the fit is supposed to start and stop.
715            Caution: This range is inclusive as opposed to standard python indexing.
716            `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6.
717            If not specified, self.prange or all timeslices are used.
718        silent : bool
719            Decides whether output is printed to the standard output.
720        '''
721        if self.N != 1:
722            raise Exception("Correlator must be projected before fitting")
723
724        if fitrange is None:
725            if self.prange:
726                fitrange = self.prange
727            else:
728                fitrange = [0, self.T - 1]
729        else:
730            if not isinstance(fitrange, list):
731                raise Exception("fitrange has to be a list with two elements")
732            if len(fitrange) != 2:
733                raise Exception("fitrange has to have exactly two elements [fit_start, fit_stop]")
734
735        xs = [x for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None]
736        ys = [self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None]
737        result = least_squares(xs, ys, function, silent=silent, **kwargs)
738        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):
740    def plateau(self, plateau_range=None, method="fit", auto_gamma=False):
741        """ Extract a plateau value from a Corr object
742
743        Parameters
744        ----------
745        plateau_range : list
746            list with two entries, indicating the first and the last timeslice
747            of the plateau region.
748        method : str
749            method to extract the plateau.
750                'fit' fits a constant to the plateau region
751                'avg', 'average' or 'mean' just average over the given timeslices.
752        auto_gamma : bool
753            apply gamma_method with default parameters to the Corr. Defaults to None
754        """
755        if not plateau_range:
756            if self.prange:
757                plateau_range = self.prange
758            else:
759                raise Exception("no plateau range provided")
760        if self.N != 1:
761            raise Exception("Correlator must be projected before getting a plateau.")
762        if (all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])):
763            raise Exception("plateau is undefined at all timeslices in plateaurange.")
764        if auto_gamma:
765            self.gamma_method()
766        if method == "fit":
767            def const_func(a, t):
768                return a[0]
769            return self.fit(const_func, plateau_range)[0]
770        elif method in ["avg", "average", "mean"]:
771            returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None])
772            return returnvalue
773
774        else:
775            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):
777    def set_prange(self, prange):
778        """Sets the attribute prange of the Corr object."""
779        if not len(prange) == 2:
780            raise Exception("prange must be a list or array with two values")
781        if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))):
782            raise Exception("Start and end point must be integers")
783        if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] < prange[1]):
784            raise Exception("Start and end point must define a range in the interval 0,T")
785
786        self.prange = prange
787        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):
789    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):
790        """Plots the correlator using the tag of the correlator as label if available.
791
792        Parameters
793        ----------
794        x_range : list
795            list of two values, determining the range of the x-axis e.g. [4, 8].
796        comp : Corr or list of Corr
797            Correlator or list of correlators which are plotted for comparison.
798            The tags of these correlators are used as labels if available.
799        logscale : bool
800            Sets y-axis to logscale.
801        plateau : Obs
802            Plateau value to be visualized in the figure.
803        fit_res : Fit_result
804            Fit_result object to be visualized.
805        ylabel : str
806            Label for the y-axis.
807        save : str
808            path to file in which the figure should be saved.
809        auto_gamma : bool
810            Apply the gamma method with standard parameters to all correlators and plateau values before plotting.
811        hide_sigma : float
812            Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors.
813        references : list
814            List of floating point values that are displayed as horizontal lines for reference.
815        title : string
816            Optional title of the figure.
817        """
818        if self.N != 1:
819            raise Exception("Correlator must be projected before plotting")
820
821        if auto_gamma:
822            self.gamma_method()
823
824        if x_range is None:
825            x_range = [0, self.T - 1]
826
827        fig = plt.figure()
828        ax1 = fig.add_subplot(111)
829
830        x, y, y_err = self.plottable()
831        if hide_sigma:
832            hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
833        else:
834            hide_from = None
835        ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag)
836        if logscale:
837            ax1.set_yscale('log')
838        else:
839            if y_range is None:
840                try:
841                    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)])
842                    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)])
843                    ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)])
844                except Exception:
845                    pass
846            else:
847                ax1.set_ylim(y_range)
848        if comp:
849            if isinstance(comp, (Corr, list)):
850                for corr in comp if isinstance(comp, list) else [comp]:
851                    if auto_gamma:
852                        corr.gamma_method()
853                    x, y, y_err = corr.plottable()
854                    if hide_sigma:
855                        hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
856                    else:
857                        hide_from = None
858                    ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor'])
859            else:
860                raise Exception("'comp' must be a correlator or a list of correlators.")
861
862        if plateau:
863            if isinstance(plateau, Obs):
864                if auto_gamma:
865                    plateau.gamma_method()
866                ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau))
867                ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-')
868            else:
869                raise Exception("'plateau' must be an Obs")
870
871        if references:
872            if isinstance(references, list):
873                for ref in references:
874                    ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--')
875            else:
876                raise Exception("'references' must be a list of floating pint values.")
877
878        if self.prange:
879            ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',')
880            ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',')
881
882        if fit_res:
883            x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05)
884            ax1.plot(x_samples,
885                     fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples),
886                     ls='-', marker=',', lw=2)
887
888        ax1.set_xlabel(r'$x_0 / a$')
889        if ylabel:
890            ax1.set_ylabel(ylabel)
891        ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5])
892
893        handles, labels = ax1.get_legend_handles_labels()
894        if labels:
895            ax1.legend()
896
897        if title:
898            plt.title(title)
899
900        plt.draw()
901
902        if save:
903            if isinstance(save, str):
904                fig.savefig(save, bbox_inches='tight')
905            else:
906                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):
908    def spaghetti_plot(self, logscale=True):
909        """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations.
910
911        Parameters
912        ----------
913        logscale : bool
914            Determines whether the scale of the y-axis is logarithmic or standard.
915        """
916        if self.N != 1:
917            raise Exception("Correlator needs to be projected first.")
918
919        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]))
920        x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None]
921
922        for name in mc_names:
923            data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T
924
925            fig = plt.figure()
926            ax = fig.add_subplot(111)
927            for dat in data:
928                ax.plot(x0_vals, dat, ls='-', marker='')
929
930            if logscale is True:
931                ax.set_yscale('log')
932
933            ax.set_xlabel(r'$x_0 / a$')
934            plt.title(name)
935            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):
937    def dump(self, filename, datatype="json.gz", **kwargs):
938        """Dumps the Corr into a file of chosen type
939        Parameters
940        ----------
941        filename : str
942            Name of the file to be saved.
943        datatype : str
944            Format of the exported file. Supported formats include
945            "json.gz" and "pickle"
946        path : str
947            specifies a custom path for the file (default '.')
948        """
949        if datatype == "json.gz":
950            from .input.json import dump_to_json
951            if 'path' in kwargs:
952                file_name = kwargs.get('path') + '/' + filename
953            else:
954                file_name = filename
955            dump_to_json(self, file_name)
956        elif datatype == "pickle":
957            dump_object(self, filename, **kwargs)
958        else:
959            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):
961    def print(self, print_range=None):
962        print(self.__repr__(print_range))
def sqrt(self):
1126    def sqrt(self):
1127        return self ** 0.5
def log(self):
1129    def log(self):
1130        newcontent = [None if _check_for_none(self, item) else np.log(item) for item in self.content]
1131        return Corr(newcontent, prange=self.prange)
def exp(self):
1133    def exp(self):
1134        newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content]
1135        return Corr(newcontent, prange=self.prange)
def sin(self):
1148    def sin(self):
1149        return self._apply_func_to_corr(np.sin)
def cos(self):
1151    def cos(self):
1152        return self._apply_func_to_corr(np.cos)
def tan(self):
1154    def tan(self):
1155        return self._apply_func_to_corr(np.tan)
def sinh(self):
1157    def sinh(self):
1158        return self._apply_func_to_corr(np.sinh)
def cosh(self):
1160    def cosh(self):
1161        return self._apply_func_to_corr(np.cosh)
def tanh(self):
1163    def tanh(self):
1164        return self._apply_func_to_corr(np.tanh)
def arcsin(self):
1166    def arcsin(self):
1167        return self._apply_func_to_corr(np.arcsin)
def arccos(self):
1169    def arccos(self):
1170        return self._apply_func_to_corr(np.arccos)
def arctan(self):
1172    def arctan(self):
1173        return self._apply_func_to_corr(np.arctan)
def arcsinh(self):
1175    def arcsinh(self):
1176        return self._apply_func_to_corr(np.arcsinh)
def arccosh(self):
1178    def arccosh(self):
1179        return self._apply_func_to_corr(np.arccosh)
def arctanh(self):
1181    def arctanh(self):
1182        return self._apply_func_to_corr(np.arctanh)
def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
1217    def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
1218        r''' Project large correlation matrix to lowest states
1219
1220        This method can be used to reduce the size of an (N x N) correlation matrix
1221        to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise
1222        is still small.
1223
1224        Parameters
1225        ----------
1226        Ntrunc: int
1227            Rank of the target matrix.
1228        tproj: int
1229            Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method.
1230            The default value is 3.
1231        t0proj: int
1232            Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly
1233            discouraged for O(a) improved theories, since the correctness of the procedure
1234            cannot be granted in this case. The default value is 2.
1235        basematrix : Corr
1236            Correlation matrix that is used to determine the eigenvectors of the
1237            lowest states based on a GEVP. basematrix is taken to be the Corr itself if
1238            is is not specified.
1239
1240        Notes
1241        -----
1242        We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving
1243        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}$
1244        and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the
1245        resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via
1246        $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large
1247        correlation matrix and to remove some noise that is added by irrelevant operators.
1248        This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated
1249        bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$.
1250        '''
1251
1252        if self.N == 1:
1253            raise Exception('Method cannot be applied to one-dimensional correlators.')
1254        if basematrix is None:
1255            basematrix = self
1256        if Ntrunc >= basematrix.N:
1257            raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N))
1258        if basematrix.N != self.N:
1259            raise Exception('basematrix and targetmatrix have to be of the same size.')
1260
1261        evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc]
1262
1263        tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object)
1264        rmat = []
1265        for t in range(basematrix.T):
1266            for i in range(Ntrunc):
1267                for j in range(Ntrunc):
1268                    tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j]
1269            rmat.append(np.copy(tmpmat))
1270
1271        newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)]
1272        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)$.