pyerrors.correlators

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