diff --git a/pyerrors/correlators.py b/pyerrors/correlators.py index b49256b0..feb91147 100644 --- a/pyerrors/correlators.py +++ b/pyerrors/correlators.py @@ -1,15 +1,12 @@ import numpy as np import autograd.numpy as anp -#from scipy.special.orthogonal import _IntegerType -from .pyerrors import * +import scipy.linalg +from .pyerrors import Obs, dump_object from .fits import standard_fit -from .linalg import * +from .linalg import eigh, mat_mat_op from .roots import find_root -from matplotlib import pyplot as plt -from matplotlib.ticker import NullFormatter # useful for `logit` scale -from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg -#import PySimpleGUI as sg -import matplotlib +import matplotlib.pyplot as plt + class Corr: """The class for a correlator (time dependent sequence of pe.Obs). @@ -25,41 +22,41 @@ class Corr: """ def __init__(self, data_input, padding_front=0, padding_back=0, prange=None): - #All data_input should be a list of things at different timeslices. This needs to be verified + # All data_input should be a list of things at different timeslices. This needs to be verified if not (isinstance(data_input, list)): raise TypeError('Corr__init__ expects a list of timeslices.') # data_input can have multiple shapes. The simplest one is a list of Obs. - #We check, if this is the case + # We check, if this is the case if all([isinstance(item, Obs) for item in data_input]): - self.content=[np.asarray([item]) for item in data_input] - #Wrapping the Obs in an array ensures that the data structure is consistent with smearing matrices. - self.N = 1 # number of smearings + self.content = [np.asarray([item]) for item in data_input] + # Wrapping the Obs in an array ensures that the data structure is consistent with smearing matrices. + self.N = 1 # number of smearings - #data_input in the form [np.array(Obs,NxN)] - elif all([isinstance(item,np.ndarray) or item==None for item in data_input]) and any([isinstance(item,np.ndarray)for item in data_input]): + # data_input in the form [np.array(Obs,NxN)] + 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]): self.content = data_input - noNull=[a for a in self.content if not (a is None)] #To check if the matrices are correct for all undefined elements + noNull = [a for a in self.content if not (a is None)] # To check if the matrices are correct for all undefined elements self.N = noNull[0].shape[0] # The checks are now identical to the case above if self.N > 1 and noNull[0].shape[0] != noNull[0].shape[1]: raise Exception("Smearing matrices are not NxN") - if (not all([item.shape == noNull[0].shape for item in noNull])): + if (not all([item.shape == noNull[0].shape for item in noNull])): raise Exception("Items in data_input are not of identical shape." + str(noNull)) - else: # In case its a list of something else. - raise Exception ("data_input contains item of wrong type") + else: # In case its a list of something else. + raise Exception("data_input contains item of wrong type") self.tag = None - #We now apply some padding to our list. In case that our list represents a correlator of length T but is not defined at every value. - #An undefined timeslice is represented by the None object + # We now apply some padding to our list. In case that our list represents a correlator of length T but is not defined at every value. + # An undefined timeslice is represented by the None object self.content = [None] * padding_front + self.content + [None] * padding_back - self.T = len(self.content) #for convenience: will be used a lot + self.T = len(self.content) # for convenience: will be used a lot - #The attribute "range" [start,end] marks a range of two timeslices. - #This is useful for keeping track of plateaus and fitranges. - #The range can be inherited from other Corrs, if the operation should not alter a chosen range eg. multiplication with a constant. + # The attribute "range" [start,end] marks a range of two timeslices. + # This is useful for keeping track of plateaus and fitranges. + # The range can be inherited from other Corrs, if the operation should not alter a chosen range eg. multiplication with a constant. self.prange = prange self.gamma_method() @@ -72,64 +69,64 @@ class Corr: else: for i in range(self.N): for j in range(self.N): - item[i,j].gamma_method() + item[i, j].gamma_method() - #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 + # 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 projected(self, vector_l=None, vector_r=None): if self.N == 1: raise Exception("Trying to project a Corr, that already has N=1.") - #This Exception is in no way necessary. One could just return self - #But there is no scenario, where a user would want that to happen and the error message might be more informative. + # This Exception is in no way necessary. One could just return self + # But there is no scenario, where a user would want that to happen and the error message might be more informative. self.gamma_method() if vector_l is None: - vector_l,vector_r=np.asarray([1.] + (self.N - 1) * [0.]),np.asarray([1.] + (self.N - 1) * [0.]) + vector_l, vector_r = np.asarray([1.] + (self.N - 1) * [0.]), np.asarray([1.] + (self.N - 1) * [0.]) elif(vector_r is None): - vector_r=vector_l + vector_r = vector_l if not vector_l.shape == vector_r.shape == (self.N,): raise Exception("Vectors are of wrong shape!") - #We always normalize before projecting! But we only raise a warning, when it is clear, they where not meant to be normalized. - if (not(0.95 < vector_r@vector_r < 1.05)) or (not(0.95 < vector_l@vector_l < 1.05)): + # We always normalize before projecting! But we only raise a warning, when it is clear, they where not meant to be normalized. + if (not (0.95 < vector_r @ vector_r < 1.05)) or (not (0.95 < vector_l @ vector_l < 1.05)): print("Vectors are normalized before projection!") - vector_l,vector_r = vector_l / np.sqrt((vector_l@vector_l)), vector_r / np.sqrt(vector_r@vector_r) + vector_l, vector_r = vector_l / np.sqrt((vector_l @ vector_l)), vector_r / np.sqrt(vector_r @ vector_r) - newcontent = [None if (item is None) else np.asarray([vector_l.T@item@vector_r]) for item in self.content] + newcontent = [None if (item is None) else np.asarray([vector_l.T @ item @ vector_r]) for item in self.content] return Corr(newcontent) def sum(self): - return np.sqrt(self.N)*self.projected(np.ones(self.N)) - #For purposes of debugging and verification, one might want to see a single smearing level. smearing will return a Corr at the specified i,j. where both are integers 0<=i,j 1: transposed = [None if (G is None) else G.T for G in self.content] - return 0.5 * (Corr(transposed)+self) + return 0.5 * (Corr(transposed) + self) if self.N == 1: raise Exception("Trying to symmetrize a smearing matrix, that already has N=1.") - - #We also include a simple GEVP method based on Scipy.linalg - def GEVP(self, t0, ts,state=1): + # We also include a simple GEVP method based on Scipy.linalg + def GEVP(self, t0, ts, state=1): if (self.content[t0] is None) or (self.content[ts] is None): raise Exception("Corr not defined at t0/ts") G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double") @@ -178,39 +172,35 @@ class Corr: G0[i, j] = self.content[t0][i, j].value Gt[i, j] = self.content[ts][i, j].value - sp_val,sp_vec = scipy.linalg.eig(Gt, G0) - sp_vec=sp_vec[:,np.argsort(sp_val)[-state]] #we only want the eigenvector belonging to the selected state - sp_vec = sp_vec/np.sqrt(sp_vec@sp_vec) + sp_val, sp_vec = scipy.linalg.eig(Gt, G0) + sp_vec = sp_vec[:, np.argsort(sp_val)[-state]] # We only want the eigenvector belonging to the selected state + sp_vec = sp_vec / np.sqrt(sp_vec @ sp_vec) return sp_vec - - def Eigenvalue(self,t0,state=1): - G=self.smearing_symmetric() - G0=G.content[t0] + def Eigenvalue(self, t0, state=1): + G = self.smearing_symmetric() + G0 = G.content[t0] L = mat_mat_op(anp.linalg.cholesky, G0) Li = mat_mat_op(anp.linalg.inv, L) - LT=L.T - LTi=mat_mat_op(anp.linalg.inv, LT) - newcontent=[] + LT = L.T + LTi = mat_mat_op(anp.linalg.inv, LT) + newcontent = [] for t in range(self.T): - Gt=G.content[t] - M=Li@Gt@LTi + Gt = G.content[t] + M = Li @ Gt @ LTi eigenvalues = eigh(M)[0] - #print(eigenvalues) - eigenvalue=eigenvalues[-state] + eigenvalue = eigenvalues[-state] newcontent.append(eigenvalue) return Corr(newcontent) - def roll(self, dt): return Corr(list(np.roll(np.array(self.content, dtype=object), dt))) - - def deriv(self, symmetric=True): #Defaults to symmetric derivative + def deriv(self, symmetric=True): # Defaults to symmetric derivative if not symmetric: newcontent = [] for t in range(self.T - 1): - if (self.content[t] is None) or (self.content[t+1] is None): + if (self.content[t] is None) or (self.content[t + 1] is None): newcontent.append(None) else: newcontent.append(self.content[t + 1] - self.content[t]) @@ -219,8 +209,8 @@ class Corr: return Corr(newcontent, padding_back=1) if symmetric: newcontent = [] - for t in range(1, self.T-1): - if (self.content[t-1] is None) or (self.content[t+1] is None): + for t in range(1, self.T - 1): + if (self.content[t - 1] is None) or (self.content[t + 1] is None): newcontent.append(None) else: newcontent.append(0.5 * (self.content[t + 1] - self.content[t - 1])) @@ -228,11 +218,10 @@ class Corr: raise Exception('Derivative is undefined at all timeslices') return Corr(newcontent, padding_back=1, padding_front=1) - def second_deriv(self): newcontent = [] - for t in range(1, self.T-1): - if (self.content[t-1] is None) or (self.content[t+1] is None): + for t in range(1, self.T - 1): + if (self.content[t - 1] is None) or (self.content[t + 1] is None): newcontent.append(None) else: newcontent.append((self.content[t + 1] - 2 * self.content[t] + self.content[t - 1])) @@ -240,7 +229,6 @@ class Corr: raise Exception("Derivative is undefined at all timeslices") return Corr(newcontent, padding_back=1, padding_front=1) - def m_eff(self, variant='log', guess=1.0): """Returns the effective mass of the correlator as correlator object @@ -252,7 +240,7 @@ class Corr: """ if self.N != 1: raise Exception('Correlator must be projected before getting m_eff') - if variant is 'log': + if variant == 'log': newcontent = [] for t in range(self.T - 1): if (self.content[t] is None) or (self.content[t + 1] is None): @@ -264,42 +252,42 @@ class Corr: return np.log(Corr(newcontent, padding_back=1)) - elif variant is 'periodic': + elif variant == 'periodic': newcontent = [] for t in range(self.T - 1): if (self.content[t] is None) or (self.content[t + 1] is None): newcontent.append(None) else: - func = lambda x, d : anp.cosh(x * (t - self.T / 2)) / anp.cosh(x * (t + 1 - self.T / 2)) - d + func = lambda x, d: anp.cosh(x * (t - self.T / 2)) / anp.cosh(x * (t + 1 - self.T / 2)) - d newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], func, guess=guess))) if(all([x is None for x in newcontent])): raise Exception('m_eff is undefined at all timeslices') return Corr(newcontent, padding_back=1) - elif variant is 'arccosh': - newcontent=[] - for t in range(1,self.T-1): - if (self.content[t] is None) or (self.content[t+1] is None)or (self.content[t-1] is None): + + elif variant == 'arccosh': + newcontent = [] + for t in range(1, self.T - 1): + if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t - 1] is None): newcontent.append(None) else: - newcontent.append((self.content[t+1]+self.content[t-1])/(2*self.content[t])) + newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t])) if(all([x is None for x in newcontent])): raise Exception("m_eff is undefined at all timeslices") - return np.arccosh(Corr(newcontent,padding_back=1,padding_front=1)) + return np.arccosh(Corr(newcontent, padding_back=1, padding_front=1)) else: raise Exception('Unkown variant.') - #We want to apply a pe.standard_fit directly to the Corr using an arbitrary function and range. + # We want to apply a pe.standard_fit directly to the Corr using an arbitrary function and range. def fit(self, function, fitrange=None, silent=False, **kwargs): if self.N != 1: raise Exception("Correlator must be projected before fitting") - #The default behaviour is: - #1 use explicit fitrange - # if none is provided, use the range of the corr - # if this is also not set, use the whole length of the corr (This could come with a warning!) - + # The default behaviour is: + # 1 use explicit fitrange + # if none is provided, use the range of the corr + # if this is also not set, use the whole length of the corr (This could come with a warning!) if fitrange is None: if self.prange: @@ -311,15 +299,13 @@ class Corr: ys = [self.content[x][0] for x in range(fitrange[0], fitrange[1]) if not self.content[x] is None] result = standard_fit(xs, ys, function, silent=silent, **kwargs) if isinstance(result, list): - [item.gamma_method() for item in result if isinstance(item,Obs)] + [item.gamma_method() for item in result if isinstance(item, Obs)] elif isinstance(result, dict): - [item.gamma_method() for item in result['fit_parameters'] if isinstance(item,Obs)] + [item.gamma_method() for item in result['fit_parameters'] if isinstance(item, Obs)] else: raise Exception('Unexpected fit result.') return result - - #we want to quickly get a plateau def plateau(self, plateau_range=None, method="fit"): if not plateau_range: if self.prange: @@ -329,13 +315,13 @@ class Corr: if self.N != 1: raise Exception("Correlator must be projected before getting a plateau.") if(all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1])])): - raise Exception("plateau is undefined at all timeslices in plateaurange.") + raise Exception("plateau is undefined at all timeslices in plateaurange.") if method == "fit": def const_func(a, t): - return a[0] # At some point pe.standard fit had an issue with single parameter fits. Being careful does not hurt - return self.fit(const_func,plateau_range)[0] - elif method in ["avg","average","mean"]: - returnvalue= np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1]+1] if not item is None]) + return a[0] # At some point pe.standard fit had an issue with single parameter fits. Being careful does not hurt + return self.fit(const_func, plateau_range)[0] + elif method in ["avg", "average", "mean"]: + returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None]) returnvalue.gamma_method() return returnvalue @@ -343,11 +329,11 @@ class Corr: raise Exception("Unsupported plateau method: " + method) def set_prange(self, prange): - if not len(prange)==2: + if not len(prange) == 2: raise Exception("prange must be a list or array with two values") - if not ((isinstance(prange[0],int)) and (isinstance(prange[1],int))): + if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))): raise Exception("Start and end point must be integers") - if not (0<=prange[0]<=self.T and 0<=prange[1]<=self.T and prange[0]