pyerrors.correlators

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

Initialize a Corr object.

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

Apply the gamma method to the content of the Corr.

#   def projected(self, vector_l=None, vector_r=None, normalize=False):
View Source
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 (item is None) 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 (self.content[t] is None 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):
View Source
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):
View Source
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):
View Source
198    def symmetric(self):
199        """ Symmetrize the correlator around x0=0."""
200        if self.T % 2 != 0:
201            raise Exception("Can not symmetrize odd T")
202
203        if np.argmax(np.abs(self.content)) != 0:
204            warnings.warn("Correlator does not seem to be symmetric around x0=0.", RuntimeWarning)
205
206        newcontent = [self.content[0]]
207        for t in range(1, self.T):
208            if (self.content[t] is None) or (self.content[self.T - t] is None):
209                newcontent.append(None)
210            else:
211                newcontent.append(0.5 * (self.content[t] + self.content[self.T - t]))
212        if(all([x is None for x in newcontent])):
213            raise Exception("Corr could not be symmetrized: No redundant values")
214        return Corr(newcontent, prange=self.prange)

Symmetrize the correlator around x0=0.

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

Anti-symmetrize the correlator around x0=0.

#   def matrix_symmetric(self):
View Source
236    def matrix_symmetric(self):
237        """Symmetrizes the correlator matrices on every timeslice."""
238        if self.N > 1:
239            transposed = [None if (G is None) else G.T for G in self.content]
240            return 0.5 * (Corr(transposed) + self)
241        if self.N == 1:
242            raise Exception("Trying to symmetrize a correlator matrix, that already has N=1.")

Symmetrizes the correlator matrices on every timeslice.

#   def GEVP(self, t0, ts=None, state=0, sorted_list='Eigenvalue'):
View Source
244    def GEVP(self, t0, ts=None, state=0, sorted_list="Eigenvalue"):
245        """Solve the general eigenvalue problem on the current correlator
246
247        Parameters
248        ----------
249        t0 : int
250            The time t0 for G(t)v= lambda G(t_0)v
251        ts : int
252            fixed time G(t_s)v= lambda G(t_0)v  if return_list=False
253            If return_list=True and sorting=Eigenvector it gives a reference point for the sorting method.
254        state : int
255            The state one is interested in, ordered by energy. The lowest state is zero.
256        sorted_list : string
257            if this argument is set, a list of vectors (len=self.T) is returned. If it is left as None, only one vector is returned.
258             "Eigenvalue"  -  The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice.
259             "Eigenvector" -  Use the method described in arXiv:2004.10472 [hep-lat] to find the set of v(t) belonging to the state.
260                              The reference state is identified by its eigenvalue at t=ts
261        """
262        symmetric_corr = self.matrix_symmetric()
263        if sorted_list is None:
264            if (ts is None):
265                raise Exception("ts is required if sorted_list=None")
266            if (self.content[t0] is None) or (self.content[ts] is None):
267                raise Exception("Corr not defined at t0/ts")
268            G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double")
269            for i in range(self.N):
270                for j in range(self.N):
271                    G0[i, j] = symmetric_corr[t0][i, j].value
272                    Gt[i, j] = symmetric_corr[ts][i, j].value
273
274            sp_vecs = _GEVP_solver(Gt, G0)
275            sp_vec = sp_vecs[state]
276            return sp_vec
277        elif sorted_list in ["Eigenvalue", "Eigenvector"]:
278            all_vecs = []
279            for t in range(self.T):
280                try:
281                    G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double")
282                    for i in range(self.N):
283                        for j in range(self.N):
284                            G0[i, j] = symmetric_corr[t0][i, j].value
285                            Gt[i, j] = symmetric_corr[t][i, j].value
286
287                    sp_vecs = _GEVP_solver(Gt, G0)
288                    if sorted_list == "Eigenvalue":
289                        sp_vec = sp_vecs[state]
290                        all_vecs.append(sp_vec)
291                    else:
292                        all_vecs.append(sp_vecs)
293                except Exception:
294                    all_vecs.append(None)
295            if sorted_list == "Eigenvector":
296                if (ts is None):
297                    raise Exception("ts is required for the Eigenvector sorting method.")
298                all_vecs = _sort_vectors(all_vecs, ts)
299                all_vecs = [a[state] for a in all_vecs]
300        else:
301            raise Exception("Unkown value for 'sorted_list'.")
302
303        return all_vecs

Solve the general eigenvalue problem on the current correlator

Parameters
  • t0 (int): The time t0 for G(t)v= lambda G(t_0)v
  • ts (int): fixed time G(t_s)v= lambda G(t_0)v if return_list=False If return_list=True and sorting=Eigenvector it gives a reference point for the sorting method.
  • state (int): The state one is interested in, ordered by energy. The lowest state is zero.
  • sorted_list (string): if this argument is set, a list of vectors (len=self.T) is returned. If it is left as 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 [hep-lat] to find the set of v(t) belonging to the state. The reference state is identified by its eigenvalue at t=ts
#   def Eigenvalue(self, t0, ts=None, state=0, sorted_list=None):
View Source
305    def Eigenvalue(self, t0, ts=None, state=0, sorted_list=None):
306        """Determines the eigenvalue of the GEVP by solving and projecting the correlator
307
308        Parameters
309        ----------
310        t0 : int
311            The time t0 for G(t)v= lambda G(t_0)v
312        ts : int
313            fixed time G(t_s)v= lambda G(t_0)v  if return_list=False
314            If return_list=True and sorting=Eigenvector it gives a reference point for the sorting method.
315        state : int
316            The state one is interested in ordered by energy. The lowest state is zero.
317        sorted_list : string
318            if this argument is set, a list of vectors (len=self.T) is returned. If it is left as None, only one vector is returned.
319             "Eigenvalue"  -  The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice.
320             "Eigenvector" -  Use the method described in arXiv:2004.10472 [hep-lat] to find the set of v(t) belonging to the state.
321                              The reference state is identified by its eigenvalue at t=ts
322        """
323        vec = self.GEVP(t0, ts=ts, state=state, sorted_list=sorted_list)
324        return self.projected(vec)

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

Parameters
  • t0 (int): The time t0 for G(t)v= lambda G(t_0)v
  • ts (int): fixed time G(t_s)v= lambda G(t_0)v if return_list=False If return_list=True and sorting=Eigenvector it gives a reference point for the sorting method.
  • state (int): The state one is interested in ordered by energy. The lowest state is zero.
  • sorted_list (string): if this argument is set, a list of vectors (len=self.T) is returned. If it is left as 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 [hep-lat] to find the set of v(t) belonging to the state. The reference state is identified by its eigenvalue at t=ts
#   def Hankel(self, N, periodic=False):
View Source
326    def Hankel(self, N, periodic=False):
327        """Constructs an NxN Hankel matrix
328
329        C(t) c(t+1) ... c(t+n-1)
330        C(t+1) c(t+2) ... c(t+n)
331        .................
332        C(t+(n-1)) c(t+n) ... c(t+2(n-1))
333
334        Parameters
335        ----------
336        N : int
337            Dimension of the Hankel matrix
338        periodic : bool, optional
339            determines whether the matrix is extended periodically
340        """
341
342        if self.N != 1:
343            raise Exception("Multi-operator Prony not implemented!")
344
345        array = np.empty([N, N], dtype="object")
346        new_content = []
347        for t in range(self.T):
348            new_content.append(array.copy())
349
350        def wrap(i):
351            while i >= self.T:
352                i -= self.T
353            return i
354
355        for t in range(self.T):
356            for i in range(N):
357                for j in range(N):
358                    if periodic:
359                        new_content[t][i, j] = self.content[wrap(t + i + j)][0]
360                    elif (t + i + j) >= self.T:
361                        new_content[t] = None
362                    else:
363                        new_content[t][i, j] = self.content[t + i + j][0]
364
365        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):
View Source
367    def roll(self, dt):
368        """Periodically shift the correlator by dt timeslices
369
370        Parameters
371        ----------
372        dt : int
373            number of timeslices
374        """
375        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):
View Source
377    def reverse(self):
378        """Reverse the time ordering of the Corr"""
379        return Corr(self.content[:: -1])

Reverse the time ordering of the Corr

#   def thin(self, spacing=2, offset=0):
View Source
381    def thin(self, spacing=2, offset=0):
382        """Thin out a correlator to suppress correlations
383
384        Parameters
385        ----------
386        spacing : int
387            Keep only every 'spacing'th entry of the correlator
388        offset : int
389            Offset the equal spacing
390        """
391        new_content = []
392        for t in range(self.T):
393            if (offset + t) % spacing != 0:
394                new_content.append(None)
395            else:
396                new_content.append(self.content[t])
397        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):
View Source
399    def correlate(self, partner):
400        """Correlate the correlator with another correlator or Obs
401
402        Parameters
403        ----------
404        partner : Obs or Corr
405            partner to correlate the correlator with.
406            Can either be an Obs which is correlated with all entries of the
407            correlator or a Corr of same length.
408        """
409        new_content = []
410        for x0, t_slice in enumerate(self.content):
411            if t_slice is None:
412                new_content.append(None)
413            else:
414                if isinstance(partner, Corr):
415                    if partner.content[x0] is None:
416                        new_content.append(None)
417                    else:
418                        new_content.append(np.array([correlate(o, partner.content[x0][0]) for o in t_slice]))
419                elif isinstance(partner, Obs):  # Should this include CObs?
420                    new_content.append(np.array([correlate(o, partner) for o in t_slice]))
421                else:
422                    raise Exception("Can only correlate with an Obs or a Corr.")
423
424        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):
View Source
426    def reweight(self, weight, **kwargs):
427        """Reweight the correlator.
428
429        Parameters
430        ----------
431        weight : Obs
432            Reweighting factor. An Observable that has to be defined on a superset of the
433            configurations in obs[i].idl for all i.
434        all_configs : bool
435            if True, the reweighted observables are normalized by the average of
436            the reweighting factor on all configurations in weight.idl and not
437            on the configurations in obs[i].idl.
438        """
439        new_content = []
440        for t_slice in self.content:
441            if t_slice is None:
442                new_content.append(None)
443            else:
444                new_content.append(np.array(reweight(weight, t_slice, **kwargs)))
445        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):
View Source
447    def T_symmetry(self, partner, parity=+1):
448        """Return the time symmetry average of the correlator and its partner
449
450        Parameters
451        ----------
452        partner : Corr
453            Time symmetry partner of the Corr
454        partity : int
455            Parity quantum number of the correlator, can be +1 or -1
456        """
457        if not isinstance(partner, Corr):
458            raise Exception("T partner has to be a Corr object.")
459        if parity not in [+1, -1]:
460            raise Exception("Parity has to be +1 or -1.")
461        T_partner = parity * partner.reverse()
462
463        t_slices = []
464        test = (self - T_partner)
465        test.gamma_method()
466        for x0, t_slice in enumerate(test.content):
467            if t_slice is not None:
468                if not t_slice[0].is_zero_within_error(5):
469                    t_slices.append(x0)
470        if t_slices:
471            warnings.warn("T symmetry partners do not agree within 5 sigma on time slices " + str(t_slices) + ".", RuntimeWarning)
472
473        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'):
View Source
475    def deriv(self, variant="symmetric"):
476        """Return the first derivative of the correlator with respect to x0.
477
478        Parameters
479        ----------
480        variant : str
481            decides which definition of the finite differences derivative is used.
482            Available choice: symmetric, forward, backward, improved, default: symmetric
483        """
484        if variant == "symmetric":
485            newcontent = []
486            for t in range(1, self.T - 1):
487                if (self.content[t - 1] is None) or (self.content[t + 1] is None):
488                    newcontent.append(None)
489                else:
490                    newcontent.append(0.5 * (self.content[t + 1] - self.content[t - 1]))
491            if(all([x is None for x in newcontent])):
492                raise Exception('Derivative is undefined at all timeslices')
493            return Corr(newcontent, padding=[1, 1])
494        elif variant == "forward":
495            newcontent = []
496            for t in range(self.T - 1):
497                if (self.content[t] is None) or (self.content[t + 1] is None):
498                    newcontent.append(None)
499                else:
500                    newcontent.append(self.content[t + 1] - self.content[t])
501            if(all([x is None for x in newcontent])):
502                raise Exception("Derivative is undefined at all timeslices")
503            return Corr(newcontent, padding=[0, 1])
504        elif variant == "backward":
505            newcontent = []
506            for t in range(1, self.T):
507                if (self.content[t - 1] is None) or (self.content[t] is None):
508                    newcontent.append(None)
509                else:
510                    newcontent.append(self.content[t] - self.content[t - 1])
511            if(all([x is None for x in newcontent])):
512                raise Exception("Derivative is undefined at all timeslices")
513            return Corr(newcontent, padding=[1, 0])
514        elif variant == "improved":
515            newcontent = []
516            for t in range(2, self.T - 2):
517                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):
518                    newcontent.append(None)
519                else:
520                    newcontent.append((1 / 12) * (self.content[t - 2] - 8 * self.content[t - 1] + 8 * self.content[t + 1] - self.content[t + 2]))
521            if(all([x is None for x in newcontent])):
522                raise Exception('Derivative is undefined at all timeslices')
523            return Corr(newcontent, padding=[2, 2])
524        else:
525            raise Exception("Unknown variant.")

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

Parameters
  • variant (str): decides which definition of the finite differences derivative is used. Available choice: symmetric, forward, backward, improved, default: symmetric
#   def second_deriv(self, variant='symmetric'):
View Source
527    def second_deriv(self, variant="symmetric"):
528        """Return the second derivative of the correlator with respect to x0.
529
530        Parameters
531        ----------
532        variant : str
533            decides which definition of the finite differences derivative is used.
534            Available choice: symmetric, improved, default: symmetric
535        """
536        if variant == "symmetric":
537            newcontent = []
538            for t in range(1, self.T - 1):
539                if (self.content[t - 1] is None) or (self.content[t + 1] is None):
540                    newcontent.append(None)
541                else:
542                    newcontent.append((self.content[t + 1] - 2 * self.content[t] + self.content[t - 1]))
543            if(all([x is None for x in newcontent])):
544                raise Exception("Derivative is undefined at all timeslices")
545            return Corr(newcontent, padding=[1, 1])
546        elif variant == "improved":
547            newcontent = []
548            for t in range(2, self.T - 2):
549                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):
550                    newcontent.append(None)
551                else:
552                    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]))
553            if(all([x is None for x in newcontent])):
554                raise Exception("Derivative is undefined at all timeslices")
555            return Corr(newcontent, padding=[2, 2])
556        else:
557            raise Exception("Unknown variant.")

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

Parameters
  • variant (str): decides which definition of the finite differences derivative is used. Available choice: symmetric, improved, default: symmetric
#   def m_eff(self, variant='log', guess=1.0):
View Source
559    def m_eff(self, variant='log', guess=1.0):
560        """Returns the effective mass of the correlator as correlator object
561
562        Parameters
563        ----------
564        variant : str
565            log : uses the standard effective mass log(C(t) / C(t+1))
566            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.
567            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.
568            See, e.g., arXiv:1205.5380
569            arccosh : Uses the explicit form of the symmetrized correlator (not recommended)
570        guess : float
571            guess for the root finder, only relevant for the root variant
572        """
573        if self.N != 1:
574            raise Exception('Correlator must be projected before getting m_eff')
575        if variant == 'log':
576            newcontent = []
577            for t in range(self.T - 1):
578                if (self.content[t] is None) or (self.content[t + 1] is None):
579                    newcontent.append(None)
580                else:
581                    newcontent.append(self.content[t] / self.content[t + 1])
582            if(all([x is None for x in newcontent])):
583                raise Exception('m_eff is undefined at all timeslices')
584
585            return np.log(Corr(newcontent, padding=[0, 1]))
586
587        elif variant in ['periodic', 'cosh', 'sinh']:
588            if variant in ['periodic', 'cosh']:
589                func = anp.cosh
590            else:
591                func = anp.sinh
592
593            def root_function(x, d):
594                return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d
595
596            newcontent = []
597            for t in range(self.T - 1):
598                if (self.content[t] is None) or (self.content[t + 1] is None):
599                    newcontent.append(None)
600                # Fill the two timeslices in the middle of the lattice with their predecessors
601                elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]:
602                    newcontent.append(newcontent[-1])
603                else:
604                    newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess)))
605            if(all([x is None for x in newcontent])):
606                raise Exception('m_eff is undefined at all timeslices')
607
608            return Corr(newcontent, padding=[0, 1])
609
610        elif variant == 'arccosh':
611            newcontent = []
612            for t in range(1, self.T - 1):
613                if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t - 1] is None):
614                    newcontent.append(None)
615                else:
616                    newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t]))
617            if(all([x is None for x in newcontent])):
618                raise Exception("m_eff is undefined at all timeslices")
619            return np.arccosh(Corr(newcontent, padding=[1, 1]))
620
621        else:
622            raise Exception('Unknown variant.')

Returns the effective mass of the correlator as correlator object

Parameters
  • variant (str): log : uses the standard effective mass log(C(t) / C(t+1)) cosh, periodic : Use periodicitiy of the correlator by solving C(t) / C(t+1) = cosh(m * (t - T/2)) / cosh(m * (t + 1 - T/2)) for m. sinh : Use anti-periodicitiy of the correlator by solving C(t) / C(t+1) = sinh(m * (t - T/2)) / sinh(m * (t + 1 - T/2)) for m. See, e.g., arXiv:1205.5380 arccosh : Uses the explicit form of the symmetrized correlator (not recommended)
  • guess (float): guess for the root finder, only relevant for the root variant
#   def fit(self, function, fitrange=None, silent=False, **kwargs):
View Source
624    def fit(self, function, fitrange=None, silent=False, **kwargs):
625        r'''Fits function to the data
626
627        Parameters
628        ----------
629        function : obj
630            function to fit to the data. See fits.least_squares for details.
631        fitrange : list
632            Two element list containing the timeslices on which the fit is supposed to start and stop.
633            Caution: This range is inclusive as opposed to standard python indexing.
634            `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6.
635            If not specified, self.prange or all timeslices are used.
636        silent : bool
637            Decides whether output is printed to the standard output.
638        '''
639        if self.N != 1:
640            raise Exception("Correlator must be projected before fitting")
641
642        if fitrange is None:
643            if self.prange:
644                fitrange = self.prange
645            else:
646                fitrange = [0, self.T - 1]
647        else:
648            if not isinstance(fitrange, list):
649                raise Exception("fitrange has to be a list with two elements")
650            if len(fitrange) != 2:
651                raise Exception("fitrange has to have exactly two elements [fit_start, fit_stop]")
652
653        xs = [x for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None]
654        ys = [self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None]
655        result = least_squares(xs, ys, function, silent=silent, **kwargs)
656        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):
View Source
658    def plateau(self, plateau_range=None, method="fit", auto_gamma=False):
659        """ Extract a plateau value from a Corr object
660
661        Parameters
662        ----------
663        plateau_range : list
664            list with two entries, indicating the first and the last timeslice
665            of the plateau region.
666        method : str
667            method to extract the plateau.
668                'fit' fits a constant to the plateau region
669                'avg', 'average' or 'mean' just average over the given timeslices.
670        auto_gamma : bool
671            apply gamma_method with default parameters to the Corr. Defaults to None
672        """
673        if not plateau_range:
674            if self.prange:
675                plateau_range = self.prange
676            else:
677                raise Exception("no plateau range provided")
678        if self.N != 1:
679            raise Exception("Correlator must be projected before getting a plateau.")
680        if(all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])):
681            raise Exception("plateau is undefined at all timeslices in plateaurange.")
682        if auto_gamma:
683            self.gamma_method()
684        if method == "fit":
685            def const_func(a, t):
686                return a[0]
687            return self.fit(const_func, plateau_range)[0]
688        elif method in ["avg", "average", "mean"]:
689            returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None])
690            return returnvalue
691
692        else:
693            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):
View Source
695    def set_prange(self, prange):
696        """Sets the attribute prange of the Corr object."""
697        if not len(prange) == 2:
698            raise Exception("prange must be a list or array with two values")
699        if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))):
700            raise Exception("Start and end point must be integers")
701        if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] < prange[1]):
702            raise Exception("Start and end point must define a range in the interval 0,T")
703
704        self.prange = prange
705        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 ):
View Source
707    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):
708        """Plots the correlator using the tag of the correlator as label if available.
709
710        Parameters
711        ----------
712        x_range : list
713            list of two values, determining the range of the x-axis e.g. [4, 8]
714        comp : Corr or list of Corr
715            Correlator or list of correlators which are plotted for comparison.
716            The tags of these correlators are used as labels if available.
717        logscale : bool
718            Sets y-axis to logscale
719        plateau : Obs
720            Plateau value to be visualized in the figure
721        fit_res : Fit_result
722            Fit_result object to be visualized
723        ylabel : str
724            Label for the y-axis
725        save : str
726            path to file in which the figure should be saved
727        auto_gamma : bool
728            Apply the gamma method with standard parameters to all correlators and plateau values before plotting.
729        hide_sigma : float
730            Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors.
731        references : list
732            List of floating point values that are displayed as horizontal lines for reference.
733        """
734        if self.N != 1:
735            raise Exception("Correlator must be projected before plotting")
736
737        if auto_gamma:
738            self.gamma_method()
739
740        if x_range is None:
741            x_range = [0, self.T - 1]
742
743        fig = plt.figure()
744        ax1 = fig.add_subplot(111)
745
746        x, y, y_err = self.plottable()
747        if hide_sigma:
748            hide_from = np.argmax((hide_sigma * np.array(y_err)) > np.abs(y)) - 1
749        else:
750            hide_from = None
751        ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag)
752        if logscale:
753            ax1.set_yscale('log')
754        else:
755            if y_range is None:
756                try:
757                    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)])
758                    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)])
759                    ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)])
760                except Exception:
761                    pass
762            else:
763                ax1.set_ylim(y_range)
764        if comp:
765            if isinstance(comp, (Corr, list)):
766                for corr in comp if isinstance(comp, list) else [comp]:
767                    if auto_gamma:
768                        corr.gamma_method()
769                    x, y, y_err = corr.plottable()
770                    if hide_sigma:
771                        hide_from = np.argmax((hide_sigma * np.array(y_err)) > np.abs(y)) - 1
772                    else:
773                        hide_from = None
774                    plt.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor'])
775            else:
776                raise Exception("'comp' must be a correlator or a list of correlators.")
777
778        if plateau:
779            if isinstance(plateau, Obs):
780                if auto_gamma:
781                    plateau.gamma_method()
782                ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau))
783                ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-')
784            else:
785                raise Exception("'plateau' must be an Obs")
786
787        if references:
788            if isinstance(references, list):
789                for ref in references:
790                    ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--')
791            else:
792                raise Exception("'references' must be a list of floating pint values.")
793
794        if self.prange:
795            ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',')
796            ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',')
797
798        if fit_res:
799            x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05)
800            ax1.plot(x_samples,
801                     fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples),
802                     ls='-', marker=',', lw=2)
803
804        ax1.set_xlabel(r'$x_0 / a$')
805        if ylabel:
806            ax1.set_ylabel(ylabel)
807        ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5])
808
809        handles, labels = ax1.get_legend_handles_labels()
810        if labels:
811            ax1.legend()
812        plt.draw()
813
814        if save:
815            if isinstance(save, str):
816                fig.savefig(save)
817            else:
818                raise Exception("'save' has to be a string.")

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

Parameters
  • x_range (list): list of two values, determining the range of the x-axis e.g. [4, 8]
  • comp (Corr or list of Corr): Correlator or list of correlators which are plotted for comparison. The tags of these correlators are used as labels if available.
  • logscale (bool): Sets y-axis to logscale
  • plateau (Obs): Plateau value to be visualized in the figure
  • fit_res (Fit_result): Fit_result object to be visualized
  • ylabel (str): Label for the y-axis
  • save (str): path to file in which the figure should be saved
  • auto_gamma (bool): Apply the gamma method with standard parameters to all correlators and plateau values before plotting.
  • hide_sigma (float): Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors.
  • references (list): List of floating point values that are displayed as horizontal lines for reference.
#   def spaghetti_plot(self, logscale=True):
View Source
820    def spaghetti_plot(self, logscale=True):
821        """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations.
822
823        Parameters
824        ----------
825        logscale : bool
826            Determines whether the scale of the y-axis is logarithmic or standard.
827        """
828        if self.N != 1:
829            raise Exception("Correlator needs to be projected first.")
830
831        mc_names = list(set([item for sublist in [o[0].mc_names for o in self.content if o is not None] for item in sublist]))
832        x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None]
833
834        for name in mc_names:
835            data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T
836
837            fig = plt.figure()
838            ax = fig.add_subplot(111)
839            for dat in data:
840                ax.plot(x0_vals, dat, ls='-', marker='')
841
842            if logscale is True:
843                ax.set_yscale('log')
844
845            ax.set_xlabel(r'$x_0 / a$')
846            plt.title(name)
847            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):
View Source
849    def dump(self, filename, datatype="json.gz", **kwargs):
850        """Dumps the Corr into a file of chosen type
851        Parameters
852        ----------
853        filename : str
854            Name of the file to be saved.
855        datatype : str
856            Format of the exported file. Supported formats include
857            "json.gz" and "pickle"
858        path : str
859            specifies a custom path for the file (default '.')
860        """
861        if datatype == "json.gz":
862            from .input.json import dump_to_json
863            if 'path' in kwargs:
864                file_name = kwargs.get('path') + '/' + filename
865            else:
866                file_name = filename
867            dump_to_json(self, file_name)
868        elif datatype == "pickle":
869            dump_object(self, filename, **kwargs)
870        else:
871            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, range=[0, None]):
View Source
873    def print(self, range=[0, None]):
874        print(self.__repr__(range))
#   def sqrt(self):
View Source
1036    def sqrt(self):
1037        return self**0.5
#   def log(self):
View Source
1039    def log(self):
1040        newcontent = [None if (item is None) else np.log(item) for item in self.content]
1041        return Corr(newcontent, prange=self.prange)
#   def exp(self):
View Source
1043    def exp(self):
1044        newcontent = [None if (item is None) else np.exp(item) for item in self.content]
1045        return Corr(newcontent, prange=self.prange)
#   def sin(self):
View Source
1058    def sin(self):
1059        return self._apply_func_to_corr(np.sin)
#   def cos(self):
View Source
1061    def cos(self):
1062        return self._apply_func_to_corr(np.cos)
#   def tan(self):
View Source
1064    def tan(self):
1065        return self._apply_func_to_corr(np.tan)
#   def sinh(self):
View Source
1067    def sinh(self):
1068        return self._apply_func_to_corr(np.sinh)
#   def cosh(self):
View Source
1070    def cosh(self):
1071        return self._apply_func_to_corr(np.cosh)
#   def tanh(self):
View Source
1073    def tanh(self):
1074        return self._apply_func_to_corr(np.tanh)
#   def arcsin(self):
View Source
1076    def arcsin(self):
1077        return self._apply_func_to_corr(np.arcsin)
#   def arccos(self):
View Source
1079    def arccos(self):
1080        return self._apply_func_to_corr(np.arccos)
#   def arctan(self):
View Source
1082    def arctan(self):
1083        return self._apply_func_to_corr(np.arctan)
#   def arcsinh(self):
View Source
1085    def arcsinh(self):
1086        return self._apply_func_to_corr(np.arcsinh)
#   def arccosh(self):
View Source
1088    def arccosh(self):
1089        return self._apply_func_to_corr(np.arccosh)
#   def arctanh(self):
View Source
1091    def arctanh(self):
1092        return self._apply_func_to_corr(np.arctanh)
#   real
#   imag
#   def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
View Source
1127    def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
1128        r''' Project large correlation matrix to lowest states
1129
1130        This method can be used to reduce the size of an (N x N) correlation matrix
1131        to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise
1132        is still small.
1133
1134        Parameters
1135        ----------
1136        Ntrunc: int
1137            Rank of the target matrix.
1138        tproj: int
1139            Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method.
1140            The default value is 3.
1141        t0proj: int
1142            Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly
1143            discouraged for O(a) improved theories, since the correctness of the procedure
1144            cannot be granted in this case. The default value is 2.
1145        basematrix : Corr
1146            Correlation matrix that is used to determine the eigenvectors of the
1147            lowest states based on a GEVP. basematrix is taken to be the Corr itself if
1148            is is not specified.
1149
1150        Notes
1151        -----
1152        We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving
1153        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}$
1154        and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the
1155        resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via
1156        $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large
1157        correlation matrix and to remove some noise that is added by irrelevant operators.
1158        This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated
1159        bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$.
1160        '''
1161
1162        if self.N == 1:
1163            raise Exception('Method cannot be applied to one-dimensional correlators.')
1164        if basematrix is None:
1165            basematrix = self
1166        if Ntrunc >= basematrix.N:
1167            raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N))
1168        if basematrix.N != self.N:
1169            raise Exception('basematrix and targetmatrix have to be of the same size.')
1170
1171        evecs = []
1172        for i in range(Ntrunc):
1173            evecs.append(basematrix.GEVP(t0proj, tproj, state=i, sorted_list=None))
1174
1175        tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object)
1176        rmat = []
1177        for t in range(basematrix.T):
1178            for i in range(Ntrunc):
1179                for j in range(Ntrunc):
1180                    tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j]
1181            rmat.append(np.copy(tmpmat))
1182
1183        newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)]
1184        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)$.