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

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
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)

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
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)))

Periodically shift the correlator by dt timeslices

Parameters
  • dt (int): number of timeslices
#   def reverse(self):
View Source
376    def reverse(self):
377        """Reverse the time ordering of the Corr"""
378        return Corr(self.content[:: -1])

Reverse the time ordering of the Corr

#   def thin(self, spacing=2, offset=0):
View Source
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)

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
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)

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
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)

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
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

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
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.")

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
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.")

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
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.')

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
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

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
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)

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
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

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 ):
View Source
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):
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        """
731        if self.N != 1:
732            raise Exception("Correlator must be projected before plotting")
733
734        if auto_gamma:
735            self.gamma_method()
736
737        if x_range is None:
738            x_range = [0, self.T - 1]
739
740        fig = plt.figure()
741        ax1 = fig.add_subplot(111)
742
743        x, y, y_err = self.plottable()
744        if hide_sigma:
745            hide_from = np.argmax((hide_sigma * np.array(y_err)) > np.abs(y)) - 1
746        else:
747            hide_from = None
748        ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag)
749        if logscale:
750            ax1.set_yscale('log')
751        else:
752            if y_range is None:
753                try:
754                    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)])
755                    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)])
756                    ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)])
757                except Exception:
758                    pass
759            else:
760                ax1.set_ylim(y_range)
761        if comp:
762            if isinstance(comp, (Corr, list)):
763                for corr in comp if isinstance(comp, list) else [comp]:
764                    if auto_gamma:
765                        corr.gamma_method()
766                    x, y, y_err = corr.plottable()
767                    if hide_sigma:
768                        hide_from = np.argmax((hide_sigma * np.array(y_err)) > np.abs(y)) - 1
769                    else:
770                        hide_from = None
771                    plt.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor'])
772            else:
773                raise Exception("'comp' must be a correlator or a list of correlators.")
774
775        if plateau:
776            if isinstance(plateau, Obs):
777                if auto_gamma:
778                    plateau.gamma_method()
779                ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau))
780                ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-')
781            else:
782                raise Exception("'plateau' must be an Obs")
783        if self.prange:
784            ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',')
785            ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',')
786
787        if fit_res:
788            x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05)
789            ax1.plot(x_samples,
790                     fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples),
791                     ls='-', marker=',', lw=2)
792
793        ax1.set_xlabel(r'$x_0 / a$')
794        if ylabel:
795            ax1.set_ylabel(ylabel)
796        ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5])
797
798        handles, labels = ax1.get_legend_handles_labels()
799        if labels:
800            ax1.legend()
801        plt.draw()
802
803        if save:
804            if isinstance(save, str):
805                fig.savefig(save)
806            else:
807                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.
#   def spaghetti_plot(self, logscale=True):
View Source
809    def spaghetti_plot(self, logscale=True):
810        """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations.
811
812        Parameters
813        ----------
814        logscale : bool
815            Determines whether the scale of the y-axis is logarithmic or standard.
816        """
817        if self.N != 1:
818            raise Exception("Correlator needs to be projected first.")
819
820        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]))
821        x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None]
822
823        for name in mc_names:
824            data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T
825
826            fig = plt.figure()
827            ax = fig.add_subplot(111)
828            for dat in data:
829                ax.plot(x0_vals, dat, ls='-', marker='')
830
831            if logscale is True:
832                ax.set_yscale('log')
833
834            ax.set_xlabel(r'$x_0 / a$')
835            plt.title(name)
836            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
838    def dump(self, filename, datatype="json.gz", **kwargs):
839        """Dumps the Corr into a file of chosen type
840        Parameters
841        ----------
842        filename : str
843            Name of the file to be saved.
844        datatype : str
845            Format of the exported file. Supported formats include
846            "json.gz" and "pickle"
847        path : str
848            specifies a custom path for the file (default '.')
849        """
850        if datatype == "json.gz":
851            from .input.json import dump_to_json
852            if 'path' in kwargs:
853                file_name = kwargs.get('path') + '/' + filename
854            else:
855                file_name = filename
856            dump_to_json(self, file_name)
857        elif datatype == "pickle":
858            dump_object(self, filename, **kwargs)
859        else:
860            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
862    def print(self, range=[0, None]):
863        print(self.__repr__(range))
#   def sqrt(self):
View Source
1025    def sqrt(self):
1026        return self**0.5
#   def log(self):
View Source
1028    def log(self):
1029        newcontent = [None if (item is None) else np.log(item) for item in self.content]
1030        return Corr(newcontent, prange=self.prange)
#   def exp(self):
View Source
1032    def exp(self):
1033        newcontent = [None if (item is None) else np.exp(item) for item in self.content]
1034        return Corr(newcontent, prange=self.prange)
#   def sin(self):
View Source
1047    def sin(self):
1048        return self._apply_func_to_corr(np.sin)
#   def cos(self):
View Source
1050    def cos(self):
1051        return self._apply_func_to_corr(np.cos)
#   def tan(self):
View Source
1053    def tan(self):
1054        return self._apply_func_to_corr(np.tan)
#   def sinh(self):
View Source
1056    def sinh(self):
1057        return self._apply_func_to_corr(np.sinh)
#   def cosh(self):
View Source
1059    def cosh(self):
1060        return self._apply_func_to_corr(np.cosh)
#   def tanh(self):
View Source
1062    def tanh(self):
1063        return self._apply_func_to_corr(np.tanh)
#   def arcsin(self):
View Source
1065    def arcsin(self):
1066        return self._apply_func_to_corr(np.arcsin)
#   def arccos(self):
View Source
1068    def arccos(self):
1069        return self._apply_func_to_corr(np.arccos)
#   def arctan(self):
View Source
1071    def arctan(self):
1072        return self._apply_func_to_corr(np.arctan)
#   def arcsinh(self):
View Source
1074    def arcsinh(self):
1075        return self._apply_func_to_corr(np.arcsinh)
#   def arccosh(self):
View Source
1077    def arccosh(self):
1078        return self._apply_func_to_corr(np.arccosh)
#   def arctanh(self):
View Source
1080    def arctanh(self):
1081        return self._apply_func_to_corr(np.arctanh)
#   real
#   imag
#   def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
View Source
1116    def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
1117        r''' Project large correlation matrix to lowest states
1118
1119        This method can be used to reduce the size of an (N x N) correlation matrix
1120        to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise
1121        is still small.
1122
1123        Parameters
1124        ----------
1125        Ntrunc: int
1126            Rank of the target matrix.
1127        tproj: int
1128            Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method.
1129            The default value is 3.
1130        t0proj: int
1131            Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly
1132            discouraged for O(a) improved theories, since the correctness of the procedure
1133            cannot be granted in this case. The default value is 2.
1134        basematrix : Corr
1135            Correlation matrix that is used to determine the eigenvectors of the
1136            lowest states based on a GEVP. basematrix is taken to be the Corr itself if
1137            is is not specified.
1138
1139        Notes
1140        -----
1141        We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving
1142        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}$
1143        and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the
1144        resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via
1145        $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large
1146        correlation matrix and to remove some noise that is added by irrelevant operators.
1147        This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated
1148        bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$.
1149        '''
1150
1151        if self.N == 1:
1152            raise Exception('Method cannot be applied to one-dimensional correlators.')
1153        if basematrix is None:
1154            basematrix = self
1155        if Ntrunc >= basematrix.N:
1156            raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N))
1157        if basematrix.N != self.N:
1158            raise Exception('basematrix and targetmatrix have to be of the same size.')
1159
1160        evecs = []
1161        for i in range(Ntrunc):
1162            evecs.append(basematrix.GEVP(t0proj, tproj, state=i))
1163
1164        tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object)
1165        rmat = []
1166        for t in range(basematrix.T):
1167            for i in range(Ntrunc):
1168                for j in range(Ntrunc):
1169                    tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j]
1170            rmat.append(np.copy(tmpmat))
1171
1172        newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)]
1173        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)$.