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

Apply the gamma method to the content of the Corr.

def gm(self, **kwargs):
117    def gamma_method(self, **kwargs):
118        """Apply the gamma method to the content of the Corr."""
119        for item in self.content:
120            if not (item is None):
121                if self.N == 1:
122                    item[0].gamma_method(**kwargs)
123                else:
124                    for i in range(self.N):
125                        for j in range(self.N):
126                            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):
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)

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

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

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

Outputs the correlator in a plotable format.

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

def symmetric(self):
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 self.content[0] is not None:
206            if np.argmax(np.abs([o[0].value if o is not None else 0 for o in 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        G0 = np.vectorize(lambda x: x.value)(symmetric_corr[t0])
311        np.linalg.cholesky(G0)  # Check if matrix G0 is positive-semidefinite.
312
313        if sort is None:
314            if (ts is None):
315                raise Exception("ts is required if sort=None.")
316            if (self.content[t0] is None) or (self.content[ts] is None):
317                raise Exception("Corr not defined at t0/ts.")
318            Gt = np.vectorize(lambda x: x.value)(symmetric_corr[ts])
319            reordered_vecs = _GEVP_solver(Gt, G0)
320
321        elif sort in ["Eigenvalue", "Eigenvector"]:
322            if sort == "Eigenvalue" and ts is not None:
323                warnings.warn("ts has no effect when sorting by eigenvalue is chosen.", RuntimeWarning)
324            all_vecs = [None] * (t0 + 1)
325            for t in range(t0 + 1, self.T):
326                try:
327                    Gt = np.vectorize(lambda x: x.value)(symmetric_corr[t])
328                    all_vecs.append(_GEVP_solver(Gt, G0))
329                except Exception:
330                    all_vecs.append(None)
331            if sort == "Eigenvector":
332                if ts is None:
333                    raise Exception("ts is required for the Eigenvector sorting method.")
334                all_vecs = _sort_vectors(all_vecs, ts)
335
336            reordered_vecs = [[v[s] if v is not None else None for v in all_vecs] for s in range(self.N)]
337        else:
338            raise Exception("Unkown value for 'sort'.")
339
340        if "state" in kwargs:
341            return reordered_vecs[kwargs.get("state")]
342        else:
343            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'):
345    def Eigenvalue(self, t0, ts=None, state=0, sort="Eigenvalue"):
346        """Determines the eigenvalue of the GEVP by solving and projecting the correlator
347
348        Parameters
349        ----------
350        state : int
351            The state one is interested in ordered by energy. The lowest state is zero.
352
353        All other parameters are identical to the ones of Corr.GEVP.
354        """
355        vec = self.GEVP(t0, ts=ts, sort=sort)[state]
356        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):
358    def Hankel(self, N, periodic=False):
359        """Constructs an NxN Hankel matrix
360
361        C(t) c(t+1) ... c(t+n-1)
362        C(t+1) c(t+2) ... c(t+n)
363        .................
364        C(t+(n-1)) c(t+n) ... c(t+2(n-1))
365
366        Parameters
367        ----------
368        N : int
369            Dimension of the Hankel matrix
370        periodic : bool, optional
371            determines whether the matrix is extended periodically
372        """
373
374        if self.N != 1:
375            raise Exception("Multi-operator Prony not implemented!")
376
377        array = np.empty([N, N], dtype="object")
378        new_content = []
379        for t in range(self.T):
380            new_content.append(array.copy())
381
382        def wrap(i):
383            while i >= self.T:
384                i -= self.T
385            return i
386
387        for t in range(self.T):
388            for i in range(N):
389                for j in range(N):
390                    if periodic:
391                        new_content[t][i, j] = self.content[wrap(t + i + j)][0]
392                    elif (t + i + j) >= self.T:
393                        new_content[t] = None
394                    else:
395                        new_content[t][i, j] = self.content[t + i + j][0]
396
397        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):
399    def roll(self, dt):
400        """Periodically shift the correlator by dt timeslices
401
402        Parameters
403        ----------
404        dt : int
405            number of timeslices
406        """
407        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):
409    def reverse(self):
410        """Reverse the time ordering of the Corr"""
411        return Corr(self.content[:: -1])

Reverse the time ordering of the Corr

def thin(self, spacing=2, offset=0):
413    def thin(self, spacing=2, offset=0):
414        """Thin out a correlator to suppress correlations
415
416        Parameters
417        ----------
418        spacing : int
419            Keep only every 'spacing'th entry of the correlator
420        offset : int
421            Offset the equal spacing
422        """
423        new_content = []
424        for t in range(self.T):
425            if (offset + t) % spacing != 0:
426                new_content.append(None)
427            else:
428                new_content.append(self.content[t])
429        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):
431    def correlate(self, partner):
432        """Correlate the correlator with another correlator or Obs
433
434        Parameters
435        ----------
436        partner : Obs or Corr
437            partner to correlate the correlator with.
438            Can either be an Obs which is correlated with all entries of the
439            correlator or a Corr of same length.
440        """
441        if self.N != 1:
442            raise Exception("Only one-dimensional correlators can be safely correlated.")
443        new_content = []
444        for x0, t_slice in enumerate(self.content):
445            if _check_for_none(self, t_slice):
446                new_content.append(None)
447            else:
448                if isinstance(partner, Corr):
449                    if _check_for_none(partner, partner.content[x0]):
450                        new_content.append(None)
451                    else:
452                        new_content.append(np.array([correlate(o, partner.content[x0][0]) for o in t_slice]))
453                elif isinstance(partner, Obs):  # Should this include CObs?
454                    new_content.append(np.array([correlate(o, partner) for o in t_slice]))
455                else:
456                    raise Exception("Can only correlate with an Obs or a Corr.")
457
458        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):
460    def reweight(self, weight, **kwargs):
461        """Reweight the correlator.
462
463        Parameters
464        ----------
465        weight : Obs
466            Reweighting factor. An Observable that has to be defined on a superset of the
467            configurations in obs[i].idl for all i.
468        all_configs : bool
469            if True, the reweighted observables are normalized by the average of
470            the reweighting factor on all configurations in weight.idl and not
471            on the configurations in obs[i].idl.
472        """
473        if self.N != 1:
474            raise Exception("Reweighting only implemented for one-dimensional correlators.")
475        new_content = []
476        for t_slice in self.content:
477            if _check_for_none(self, t_slice):
478                new_content.append(None)
479            else:
480                new_content.append(np.array(reweight(weight, t_slice, **kwargs)))
481        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):
483    def T_symmetry(self, partner, parity=+1):
484        """Return the time symmetry average of the correlator and its partner
485
486        Parameters
487        ----------
488        partner : Corr
489            Time symmetry partner of the Corr
490        partity : int
491            Parity quantum number of the correlator, can be +1 or -1
492        """
493        if self.N != 1:
494            raise Exception("T_symmetry only implemented for one-dimensional correlators.")
495        if not isinstance(partner, Corr):
496            raise Exception("T partner has to be a Corr object.")
497        if parity not in [+1, -1]:
498            raise Exception("Parity has to be +1 or -1.")
499        T_partner = parity * partner.reverse()
500
501        t_slices = []
502        test = (self - T_partner)
503        test.gamma_method()
504        for x0, t_slice in enumerate(test.content):
505            if t_slice is not None:
506                if not t_slice[0].is_zero_within_error(5):
507                    t_slices.append(x0)
508        if t_slices:
509            warnings.warn("T symmetry partners do not agree within 5 sigma on time slices " + str(t_slices) + ".", RuntimeWarning)
510
511        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'):
513    def deriv(self, variant="symmetric"):
514        """Return the first derivative of the correlator with respect to x0.
515
516        Parameters
517        ----------
518        variant : str
519            decides which definition of the finite differences derivative is used.
520            Available choice: symmetric, forward, backward, improved, log, default: symmetric
521        """
522        if self.N != 1:
523            raise Exception("deriv only implemented for one-dimensional correlators.")
524        if variant == "symmetric":
525            newcontent = []
526            for t in range(1, self.T - 1):
527                if (self.content[t - 1] is None) or (self.content[t + 1] is None):
528                    newcontent.append(None)
529                else:
530                    newcontent.append(0.5 * (self.content[t + 1] - self.content[t - 1]))
531            if (all([x is None for x in newcontent])):
532                raise Exception('Derivative is undefined at all timeslices')
533            return Corr(newcontent, padding=[1, 1])
534        elif variant == "forward":
535            newcontent = []
536            for t in range(self.T - 1):
537                if (self.content[t] is None) or (self.content[t + 1] is None):
538                    newcontent.append(None)
539                else:
540                    newcontent.append(self.content[t + 1] - self.content[t])
541            if (all([x is None for x in newcontent])):
542                raise Exception("Derivative is undefined at all timeslices")
543            return Corr(newcontent, padding=[0, 1])
544        elif variant == "backward":
545            newcontent = []
546            for t in range(1, self.T):
547                if (self.content[t - 1] is None) or (self.content[t] is None):
548                    newcontent.append(None)
549                else:
550                    newcontent.append(self.content[t] - self.content[t - 1])
551            if (all([x is None for x in newcontent])):
552                raise Exception("Derivative is undefined at all timeslices")
553            return Corr(newcontent, padding=[1, 0])
554        elif variant == "improved":
555            newcontent = []
556            for t in range(2, self.T - 2):
557                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):
558                    newcontent.append(None)
559                else:
560                    newcontent.append((1 / 12) * (self.content[t - 2] - 8 * self.content[t - 1] + 8 * self.content[t + 1] - self.content[t + 2]))
561            if (all([x is None for x in newcontent])):
562                raise Exception('Derivative is undefined at all timeslices')
563            return Corr(newcontent, padding=[2, 2])
564        elif variant == 'log':
565            newcontent = []
566            for t in range(self.T):
567                if (self.content[t] is None) or (self.content[t] <= 0):
568                    newcontent.append(None)
569                else:
570                    newcontent.append(np.log(self.content[t]))
571            if (all([x is None for x in newcontent])):
572                raise Exception("Log is undefined at all timeslices")
573            logcorr = Corr(newcontent)
574            return self * logcorr.deriv('symmetric')
575        else:
576            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'):
578    def second_deriv(self, variant="symmetric"):
579        """Return the second derivative of the correlator with respect to x0.
580
581        Parameters
582        ----------
583        variant : str
584            decides which definition of the finite differences derivative is used.
585            Available choice: symmetric, improved, log, default: symmetric
586        """
587        if self.N != 1:
588            raise Exception("second_deriv only implemented for one-dimensional correlators.")
589        if variant == "symmetric":
590            newcontent = []
591            for t in range(1, self.T - 1):
592                if (self.content[t - 1] is None) or (self.content[t + 1] is None):
593                    newcontent.append(None)
594                else:
595                    newcontent.append((self.content[t + 1] - 2 * self.content[t] + self.content[t - 1]))
596            if (all([x is None for x in newcontent])):
597                raise Exception("Derivative is undefined at all timeslices")
598            return Corr(newcontent, padding=[1, 1])
599        elif variant == "improved":
600            newcontent = []
601            for t in range(2, self.T - 2):
602                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):
603                    newcontent.append(None)
604                else:
605                    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]))
606            if (all([x is None for x in newcontent])):
607                raise Exception("Derivative is undefined at all timeslices")
608            return Corr(newcontent, padding=[2, 2])
609        elif variant == 'log':
610            newcontent = []
611            for t in range(self.T):
612                if (self.content[t] is None) or (self.content[t] <= 0):
613                    newcontent.append(None)
614                else:
615                    newcontent.append(np.log(self.content[t]))
616            if (all([x is None for x in newcontent])):
617                raise Exception("Log is undefined at all timeslices")
618            logcorr = Corr(newcontent)
619            return self * (logcorr.second_deriv('symmetric') + (logcorr.deriv('symmetric'))**2)
620        else:
621            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):
623    def m_eff(self, variant='log', guess=1.0):
624        """Returns the effective mass of the correlator as correlator object
625
626        Parameters
627        ----------
628        variant : str
629            log : uses the standard effective mass log(C(t) / C(t+1))
630            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.
631            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.
632            See, e.g., arXiv:1205.5380
633            arccosh : Uses the explicit form of the symmetrized correlator (not recommended)
634            logsym: uses the symmetric effective mass log(C(t-1) / C(t+1))/2
635        guess : float
636            guess for the root finder, only relevant for the root variant
637        """
638        if self.N != 1:
639            raise Exception('Correlator must be projected before getting m_eff')
640        if variant == 'log':
641            newcontent = []
642            for t in range(self.T - 1):
643                if ((self.content[t] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0):
644                    newcontent.append(None)
645                elif self.content[t][0].value / self.content[t + 1][0].value < 0:
646                    newcontent.append(None)
647                else:
648                    newcontent.append(self.content[t] / self.content[t + 1])
649            if (all([x is None for x in newcontent])):
650                raise Exception('m_eff is undefined at all timeslices')
651
652            return np.log(Corr(newcontent, padding=[0, 1]))
653
654        elif variant == 'logsym':
655            newcontent = []
656            for t in range(1, self.T - 1):
657                if ((self.content[t - 1] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0):
658                    newcontent.append(None)
659                elif self.content[t - 1][0].value / self.content[t + 1][0].value < 0:
660                    newcontent.append(None)
661                else:
662                    newcontent.append(self.content[t - 1] / self.content[t + 1])
663            if (all([x is None for x in newcontent])):
664                raise Exception('m_eff is undefined at all timeslices')
665
666            return np.log(Corr(newcontent, padding=[1, 1])) / 2
667
668        elif variant in ['periodic', 'cosh', 'sinh']:
669            if variant in ['periodic', 'cosh']:
670                func = anp.cosh
671            else:
672                func = anp.sinh
673
674            def root_function(x, d):
675                return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d
676
677            newcontent = []
678            for t in range(self.T - 1):
679                if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t + 1][0].value == 0):
680                    newcontent.append(None)
681                # Fill the two timeslices in the middle of the lattice with their predecessors
682                elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]:
683                    newcontent.append(newcontent[-1])
684                elif self.content[t][0].value / self.content[t + 1][0].value < 0:
685                    newcontent.append(None)
686                else:
687                    newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess)))
688            if (all([x is None for x in newcontent])):
689                raise Exception('m_eff is undefined at all timeslices')
690
691            return Corr(newcontent, padding=[0, 1])
692
693        elif variant == 'arccosh':
694            newcontent = []
695            for t in range(1, self.T - 1):
696                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):
697                    newcontent.append(None)
698                else:
699                    newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t]))
700            if (all([x is None for x in newcontent])):
701                raise Exception("m_eff is undefined at all timeslices")
702            return np.arccosh(Corr(newcontent, padding=[1, 1]))
703
704        else:
705            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):
707    def fit(self, function, fitrange=None, silent=False, **kwargs):
708        r'''Fits function to the data
709
710        Parameters
711        ----------
712        function : obj
713            function to fit to the data. See fits.least_squares for details.
714        fitrange : list
715            Two element list containing the timeslices on which the fit is supposed to start and stop.
716            Caution: This range is inclusive as opposed to standard python indexing.
717            `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6.
718            If not specified, self.prange or all timeslices are used.
719        silent : bool
720            Decides whether output is printed to the standard output.
721        '''
722        if self.N != 1:
723            raise Exception("Correlator must be projected before fitting")
724
725        if fitrange is None:
726            if self.prange:
727                fitrange = self.prange
728            else:
729                fitrange = [0, self.T - 1]
730        else:
731            if not isinstance(fitrange, list):
732                raise Exception("fitrange has to be a list with two elements")
733            if len(fitrange) != 2:
734                raise Exception("fitrange has to have exactly two elements [fit_start, fit_stop]")
735
736        xs = [x for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None]
737        ys = [self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None]
738        result = least_squares(xs, ys, function, silent=silent, **kwargs)
739        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):
741    def plateau(self, plateau_range=None, method="fit", auto_gamma=False):
742        """ Extract a plateau value from a Corr object
743
744        Parameters
745        ----------
746        plateau_range : list
747            list with two entries, indicating the first and the last timeslice
748            of the plateau region.
749        method : str
750            method to extract the plateau.
751                'fit' fits a constant to the plateau region
752                'avg', 'average' or 'mean' just average over the given timeslices.
753        auto_gamma : bool
754            apply gamma_method with default parameters to the Corr. Defaults to None
755        """
756        if not plateau_range:
757            if self.prange:
758                plateau_range = self.prange
759            else:
760                raise Exception("no plateau range provided")
761        if self.N != 1:
762            raise Exception("Correlator must be projected before getting a plateau.")
763        if (all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])):
764            raise Exception("plateau is undefined at all timeslices in plateaurange.")
765        if auto_gamma:
766            self.gamma_method()
767        if method == "fit":
768            def const_func(a, t):
769                return a[0]
770            return self.fit(const_func, plateau_range)[0]
771        elif method in ["avg", "average", "mean"]:
772            returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None])
773            return returnvalue
774
775        else:
776            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):
778    def set_prange(self, prange):
779        """Sets the attribute prange of the Corr object."""
780        if not len(prange) == 2:
781            raise Exception("prange must be a list or array with two values")
782        if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))):
783            raise Exception("Start and end point must be integers")
784        if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] < prange[1]):
785            raise Exception("Start and end point must define a range in the interval 0,T")
786
787        self.prange = prange
788        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):
790    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):
791        """Plots the correlator using the tag of the correlator as label if available.
792
793        Parameters
794        ----------
795        x_range : list
796            list of two values, determining the range of the x-axis e.g. [4, 8].
797        comp : Corr or list of Corr
798            Correlator or list of correlators which are plotted for comparison.
799            The tags of these correlators are used as labels if available.
800        logscale : bool
801            Sets y-axis to logscale.
802        plateau : Obs
803            Plateau value to be visualized in the figure.
804        fit_res : Fit_result
805            Fit_result object to be visualized.
806        ylabel : str
807            Label for the y-axis.
808        save : str
809            path to file in which the figure should be saved.
810        auto_gamma : bool
811            Apply the gamma method with standard parameters to all correlators and plateau values before plotting.
812        hide_sigma : float
813            Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors.
814        references : list
815            List of floating point values that are displayed as horizontal lines for reference.
816        title : string
817            Optional title of the figure.
818        """
819        if self.N != 1:
820            raise Exception("Correlator must be projected before plotting")
821
822        if auto_gamma:
823            self.gamma_method()
824
825        if x_range is None:
826            x_range = [0, self.T - 1]
827
828        fig = plt.figure()
829        ax1 = fig.add_subplot(111)
830
831        x, y, y_err = self.plottable()
832        if hide_sigma:
833            hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
834        else:
835            hide_from = None
836        ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag)
837        if logscale:
838            ax1.set_yscale('log')
839        else:
840            if y_range is None:
841                try:
842                    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)])
843                    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)])
844                    ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)])
845                except Exception:
846                    pass
847            else:
848                ax1.set_ylim(y_range)
849        if comp:
850            if isinstance(comp, (Corr, list)):
851                for corr in comp if isinstance(comp, list) else [comp]:
852                    if auto_gamma:
853                        corr.gamma_method()
854                    x, y, y_err = corr.plottable()
855                    if hide_sigma:
856                        hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
857                    else:
858                        hide_from = None
859                    ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor'])
860            else:
861                raise Exception("'comp' must be a correlator or a list of correlators.")
862
863        if plateau:
864            if isinstance(plateau, Obs):
865                if auto_gamma:
866                    plateau.gamma_method()
867                ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau))
868                ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-')
869            else:
870                raise Exception("'plateau' must be an Obs")
871
872        if references:
873            if isinstance(references, list):
874                for ref in references:
875                    ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--')
876            else:
877                raise Exception("'references' must be a list of floating pint values.")
878
879        if self.prange:
880            ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',')
881            ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',')
882
883        if fit_res:
884            x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05)
885            ax1.plot(x_samples,
886                     fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples),
887                     ls='-', marker=',', lw=2)
888
889        ax1.set_xlabel(r'$x_0 / a$')
890        if ylabel:
891            ax1.set_ylabel(ylabel)
892        ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5])
893
894        handles, labels = ax1.get_legend_handles_labels()
895        if labels:
896            ax1.legend()
897
898        if title:
899            plt.title(title)
900
901        plt.draw()
902
903        if save:
904            if isinstance(save, str):
905                fig.savefig(save, bbox_inches='tight')
906            else:
907                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):
909    def spaghetti_plot(self, logscale=True):
910        """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations.
911
912        Parameters
913        ----------
914        logscale : bool
915            Determines whether the scale of the y-axis is logarithmic or standard.
916        """
917        if self.N != 1:
918            raise Exception("Correlator needs to be projected first.")
919
920        mc_names = list(set([item for sublist in [sum(map(o[0].e_content.get, o[0].mc_names), []) for o in self.content if o is not None] for item in sublist]))
921        x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None]
922
923        for name in mc_names:
924            data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T
925
926            fig = plt.figure()
927            ax = fig.add_subplot(111)
928            for dat in data:
929                ax.plot(x0_vals, dat, ls='-', marker='')
930
931            if logscale is True:
932                ax.set_yscale('log')
933
934            ax.set_xlabel(r'$x_0 / a$')
935            plt.title(name)
936            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):
938    def dump(self, filename, datatype="json.gz", **kwargs):
939        """Dumps the Corr into a file of chosen type
940        Parameters
941        ----------
942        filename : str
943            Name of the file to be saved.
944        datatype : str
945            Format of the exported file. Supported formats include
946            "json.gz" and "pickle"
947        path : str
948            specifies a custom path for the file (default '.')
949        """
950        if datatype == "json.gz":
951            from .input.json import dump_to_json
952            if 'path' in kwargs:
953                file_name = kwargs.get('path') + '/' + filename
954            else:
955                file_name = filename
956            dump_to_json(self, file_name)
957        elif datatype == "pickle":
958            dump_object(self, filename, **kwargs)
959        else:
960            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):
962    def print(self, print_range=None):
963        print(self.__repr__(print_range))
def sqrt(self):
1129    def sqrt(self):
1130        return self ** 0.5
def log(self):
1132    def log(self):
1133        newcontent = [None if _check_for_none(self, item) else np.log(item) for item in self.content]
1134        return Corr(newcontent, prange=self.prange)
def exp(self):
1136    def exp(self):
1137        newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content]
1138        return Corr(newcontent, prange=self.prange)
def sin(self):
1153    def sin(self):
1154        return self._apply_func_to_corr(np.sin)
def cos(self):
1156    def cos(self):
1157        return self._apply_func_to_corr(np.cos)
def tan(self):
1159    def tan(self):
1160        return self._apply_func_to_corr(np.tan)
def sinh(self):
1162    def sinh(self):
1163        return self._apply_func_to_corr(np.sinh)
def cosh(self):
1165    def cosh(self):
1166        return self._apply_func_to_corr(np.cosh)
def tanh(self):
1168    def tanh(self):
1169        return self._apply_func_to_corr(np.tanh)
def arcsin(self):
1171    def arcsin(self):
1172        return self._apply_func_to_corr(np.arcsin)
def arccos(self):
1174    def arccos(self):
1175        return self._apply_func_to_corr(np.arccos)
def arctan(self):
1177    def arctan(self):
1178        return self._apply_func_to_corr(np.arctan)
def arcsinh(self):
1180    def arcsinh(self):
1181        return self._apply_func_to_corr(np.arcsinh)
def arccosh(self):
1183    def arccosh(self):
1184        return self._apply_func_to_corr(np.arccosh)
def arctanh(self):
1186    def arctanh(self):
1187        return self._apply_func_to_corr(np.arctanh)
def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
1222    def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
1223        r''' Project large correlation matrix to lowest states
1224
1225        This method can be used to reduce the size of an (N x N) correlation matrix
1226        to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise
1227        is still small.
1228
1229        Parameters
1230        ----------
1231        Ntrunc: int
1232            Rank of the target matrix.
1233        tproj: int
1234            Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method.
1235            The default value is 3.
1236        t0proj: int
1237            Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly
1238            discouraged for O(a) improved theories, since the correctness of the procedure
1239            cannot be granted in this case. The default value is 2.
1240        basematrix : Corr
1241            Correlation matrix that is used to determine the eigenvectors of the
1242            lowest states based on a GEVP. basematrix is taken to be the Corr itself if
1243            is is not specified.
1244
1245        Notes
1246        -----
1247        We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving
1248        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}$
1249        and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the
1250        resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via
1251        $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large
1252        correlation matrix and to remove some noise that is added by irrelevant operators.
1253        This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated
1254        bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$.
1255        '''
1256
1257        if self.N == 1:
1258            raise Exception('Method cannot be applied to one-dimensional correlators.')
1259        if basematrix is None:
1260            basematrix = self
1261        if Ntrunc >= basematrix.N:
1262            raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N))
1263        if basematrix.N != self.N:
1264            raise Exception('basematrix and targetmatrix have to be of the same size.')
1265
1266        evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc]
1267
1268        tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object)
1269        rmat = []
1270        for t in range(basematrix.T):
1271            for i in range(Ntrunc):
1272                for j in range(Ntrunc):
1273                    tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j]
1274            rmat.append(np.copy(tmpmat))
1275
1276        newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)]
1277        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)$.