diff --git a/pyerrors/correlators.py b/pyerrors/correlators.py index b78b13d6..e6ea9caf 100644 --- a/pyerrors/correlators.py +++ b/pyerrors/correlators.py @@ -30,7 +30,7 @@ class Corr: 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 - if all([ (isinstance(item, Obs) or isinstance(item, CObs)) for item in data_input]): + if all([(isinstance(item, Obs) or isinstance(item, CObs)) 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 @@ -97,7 +97,7 @@ class Corr: # 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,normalize=False): + def projected(self, vector_l=None, vector_r=None, normalize=False): 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 @@ -109,18 +109,16 @@ class Corr: 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 - - if isinstance(vector_l,list) and not isinstance(vector_r,list): - if len(vector_l)!=self.T: + if isinstance(vector_l, list) and not isinstance(vector_r, list): + if len(vector_l) != self.T: raise Exception("Length of vector list must be equal to T") - vector_r=[vector_r]*self.T - if isinstance(vector_r,list) and not isinstance(vector_l,list): - if len(vector_r)!=self.T: + vector_r = [vector_r] * self.T + if isinstance(vector_r, list) and not isinstance(vector_l, list): + if len(vector_r) != self.T: raise Exception("Length of vector list must be equal to T") - vector_l=[vector_l]*self.T + vector_l = [vector_l] * self.T - - if not isinstance(vector_l,list): + if not isinstance(vector_l, list): if not vector_l.shape == vector_r.shape == (self.N,): raise Exception("Vectors are of wrong shape!") if normalize: @@ -215,9 +213,7 @@ class Corr: # There are two ways, the GEVP metod can be called. # 1. return_list=False will return a single eigenvector, normalized according to V*C(t_0)*V=1 # 2. return_list=True will return a new eigenvector for every timeslice. The time t_s is used to order the vectors according to. arXiv:2004.10472 [hep-lat] - - - def GEVP(self, t0, ts, state=0, sorting="Eigenvalue",return_list=False): + def GEVP(self, t0, ts, state=0, sorting="Eigenvalue", return_list=False): if not return_list: if (self.content[t0] is None) or (self.content[ts] is None): raise Exception("Corr not defined at t0/ts") @@ -227,11 +223,11 @@ class Corr: G0[i, j] = self.content[t0][i, j].value Gt[i, j] = self.content[ts][i, j].value - sp_vecs=GEVP_solver(Gt,G0) - sp_vec=sp_vecs[state] + sp_vecs = GEVP_solver(Gt, G0) + sp_vec = sp_vecs[state] return sp_vec if return_list: - all_vecs=[] + all_vecs = [] for t in range(self.T): try: G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double") @@ -240,24 +236,19 @@ class Corr: G0[i, j] = self.content[t0][i, j].value Gt[i, j] = self.content[t][i, j].value - sp_vecs = GEVP_solver(Gt,G0) - if sorting=="Eigenvalue": + sp_vecs = GEVP_solver(Gt, G0) + if sorting == "Eigenvalue": sp_vec = sp_vecs[state] all_vecs.append(sp_vec) else: all_vecs.append(sp_vecs) - - - - except: #This could contain a check for real eigenvectors + except "Failure to solve for one timeslice": # This could contain a check for real eigenvectors all_vecs.append(None) - if sorting=="Eigenvector": - all_vecs=sort_vectors(all_vecs,ts) - all_vecs=[a[state] for a in all_vecs] + if sorting == "Eigenvector": + all_vecs = sort_vectors(all_vecs, ts) + all_vecs = [a[state] for a in all_vecs] return all_vecs - - def Eigenvalue(self, t0, state=1): G = self.smearing_symmetric() @@ -278,27 +269,23 @@ class Corr: newcontent.append(eigenvalue) return Corr(newcontent) - - - def Hankel(self,N,periodic=False): - #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)) + def Hankel(self, N, periodic=False): + # 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)) - if self.N!=1: - raise Exception("Multi-operator Prony not implemented!") + if self.N != 1: + raise Exception("Multi-operator Prony not implemented!") - - array=np.empty([N,N],dtype="object") - new_content=[] + array = np.empty([N, N], dtype="object") + new_content = [] for t in range(self.T): new_content.append(array.copy()) - def wrap(i): - if i>=self.T: + if i >= self.T: return i-self.T return i @@ -306,18 +293,14 @@ class Corr: for i in range(N): for j in range(N): if periodic: - new_content[t][i,j]=self.content[wrap(t+i+j)][0] - elif (t+i+j)>=self.T: + new_content[t][i, j] = self.content[wrap(t+i+j)][0] + elif (t+i+j) >= self.T: new_content[t]=None else: - new_content[t][i,j]=self.content[t+i+j][0] + new_content[t][i, j] = self.content[t+i+j][0] - - return Corr(new_content) - - def roll(self, dt): """Periodically shift the correlator by dt timeslices @@ -701,8 +684,8 @@ class Corr: content_string += "Description: " + self.tag + "\n" if self.N!=1: return content_string - #This avoids a crash for N>1. I do not know, what else to do here. I like the list representation for N==1. We could print only one "smearing" or one matrix. Printing everything will just - #be a wall of numbers. + # This avoids a crash for N>1. I do not know, what else to do here. I like the list representation for N==1. We could print only one "smearing" or one matrix. Printing everything will just + # be a wall of numbers. if range[1]: @@ -929,7 +912,7 @@ class Corr: if isinstance(obs_OR_cobs, CObs): return obs_OR_cobs.imag else: - return obs_OR_cobs*0 #So it stays the right type + return obs_OR_cobs*0 # So it stays the right type return self._apply_func_to_corr(return_imag) @@ -948,7 +931,7 @@ class Corr: -def sort_vectors(vec_set, ts): #Helper function used to find a set of Eigenvectors consistent over all timeslices +def sort_vectors(vec_set, ts): # Helper function used to find a set of Eigenvectors consistent over all timeslices reference_sorting=np.array(vec_set[ts]) N=reference_sorting.shape[0] sorted_vec_set=[] @@ -963,7 +946,7 @@ def sort_vectors(vec_set, ts): #Helper function used to find a set of Eigenvecto for k in range(N): new_sorting=reference_sorting.copy() new_sorting[perm[k],:]=vec_set[t][k] - current_score*=abs(np.linalg.det(new_sorting)) + current_score *= abs(np.linalg.det(new_sorting)) if current_score>best_score: best_score=current_score best_perm=perm @@ -980,7 +963,7 @@ def sort_vectors(vec_set, ts): #Helper function used to find a set of Eigenvecto -def permutation(lst): #Shamelessly copied +def permutation(lst): # Shamelessly copied if len(lst) == 1: return [lst] l = [] @@ -993,8 +976,8 @@ def permutation(lst): #Shamelessly copied return l -def GEVP_solver(Gt,G0): #Just so normalization an sorting does not need to be repeated. Here we could later put in some checks +def GEVP_solver(Gt,G0): # Just so normalization an sorting does not need to be repeated. Here we could later put in some checks sp_val, sp_vecs = scipy.linalg.eig(Gt, G0) sp_vecs=[sp_vecs[:, np.argsort(sp_val)[-i]]for i in range(1,sp_vecs.shape[0]+1) ] sp_vecs=[v/np.sqrt((v.T@G0@v)) for v in sp_vecs] - return sp_vecs \ No newline at end of file + return sp_vecs