lattice/pyerrors
1
0
Fork 0
mirror of https://github.com/fjosw/pyerrors.git synced 2025-05-15 12:03:42 +02:00
Code Issues Projects Releases Packages Wiki Activity

The changes i tried to push before

Browse source
This commit is contained in:
JanNeuendorf 2022-01-18 16:16:54 +01:00
parent 39f176585e
commit 2f11d0d30b
1 changed files with 41 additions and 58 deletions
Download patch file Download diff file Expand all files Collapse all files

95
pyerrors/correlators.py
View file

@ -30,7 +30,7 @@ class Corr:
raise TypeError('Corr__init__ expects a list of timeslices.') raise TypeError('Corr__init__ expects a list of timeslices.')
# data_input can have multiple shapes. The simplest one is a list of Obs. # 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) 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] 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. # Wrapping the Obs in an array ensures that the data structure is consistent with smearing matrices.
self.N = 1 # number of smearings self.N = 1 # number of smearings
@ -97,7 +97,7 @@ class Corr:
# The method can use one or two vectors. # The method can use one or two vectors.
# If two are specified it returns v1@G@v2 (the order might be very important.) # 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 # 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: if self.N == 1:
raise Exception("Trying to project a Corr, that already has 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 # 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.]) vector_l, vector_r = np.asarray([1.] + (self.N - 1) * [0.]), np.asarray([1.] + (self.N - 1) * [0.])
elif(vector_r is None): elif(vector_r is None):
vector_r = vector_l vector_r = vector_l
if isinstance(vector_l, list) and not isinstance(vector_r, list):
if isinstance(vector_l,list) and not isinstance(vector_r,list): if len(vector_l) != self.T:
if len(vector_l)!=self.T:
raise Exception("Length of vector list must be equal to T") raise Exception("Length of vector list must be equal to T")
vector_r=[vector_r]*self.T vector_r = [vector_r] * self.T
if isinstance(vector_r,list) and not isinstance(vector_l,list): if isinstance(vector_r, list) and not isinstance(vector_l, list):
if len(vector_r)!=self.T: if len(vector_r) != self.T:
raise Exception("Length of vector list must be equal to 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,): if not vector_l.shape == vector_r.shape == (self.N,):
raise Exception("Vectors are of wrong shape!") raise Exception("Vectors are of wrong shape!")
if normalize: if normalize:
@ -215,9 +213,7 @@ class Corr:
# There are two ways, the GEVP metod can be called. # 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 # 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] # 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 not return_list:
if (self.content[t0] is None) or (self.content[ts] is None): if (self.content[t0] is None) or (self.content[ts] is None):
raise Exception("Corr not defined at t0/ts") raise Exception("Corr not defined at t0/ts")
@ -227,11 +223,11 @@ class Corr:
G0[i, j] = self.content[t0][i, j].value G0[i, j] = self.content[t0][i, j].value
Gt[i, j] = self.content[ts][i, j].value Gt[i, j] = self.content[ts][i, j].value
sp_vecs=GEVP_solver(Gt,G0) sp_vecs = GEVP_solver(Gt, G0)
sp_vec=sp_vecs[state] sp_vec = sp_vecs[state]
return sp_vec return sp_vec
if return_list: if return_list:
all_vecs=[] all_vecs = []
for t in range(self.T): for t in range(self.T):
try: try:
G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double") G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double")
@ -240,25 +236,20 @@ class Corr:
G0[i, j] = self.content[t0][i, j].value G0[i, j] = self.content[t0][i, j].value
Gt[i, j] = self.content[t][i, j].value Gt[i, j] = self.content[t][i, j].value
sp_vecs = GEVP_solver(Gt,G0) sp_vecs = GEVP_solver(Gt, G0)
if sorting=="Eigenvalue": if sorting == "Eigenvalue":
sp_vec = sp_vecs[state] sp_vec = sp_vecs[state]
all_vecs.append(sp_vec) all_vecs.append(sp_vec)
else: else:
all_vecs.append(sp_vecs) all_vecs.append(sp_vecs)
except "Failure to solve for one timeslice": # This could contain a check for real eigenvectors
except: #This could contain a check for real eigenvectors
all_vecs.append(None) all_vecs.append(None)
if sorting=="Eigenvector": if sorting == "Eigenvector":
all_vecs=sort_vectors(all_vecs,ts) all_vecs = sort_vectors(all_vecs, ts)
all_vecs=[a[state] for a in all_vecs] all_vecs = [a[state] for a in all_vecs]
return all_vecs return all_vecs
def Eigenvalue(self, t0, state=1): def Eigenvalue(self, t0, state=1):
G = self.smearing_symmetric() G = self.smearing_symmetric()
G0 = G.content[t0] G0 = G.content[t0]
@ -278,27 +269,23 @@ class Corr:
newcontent.append(eigenvalue) newcontent.append(eigenvalue)
return Corr(newcontent) 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))
if self.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!") raise Exception("Multi-operator Prony not implemented!")
array = np.empty([N, N], dtype="object")
array=np.empty([N,N],dtype="object") new_content = []
new_content=[]
for t in range(self.T): for t in range(self.T):
new_content.append(array.copy()) new_content.append(array.copy())
def wrap(i): def wrap(i):
if i>=self.T: if i >= self.T:
return i-self.T return i-self.T
return i return i
@ -306,18 +293,14 @@ class Corr:
for i in range(N): for i in range(N):
for j in range(N): for j in range(N):
if periodic: if periodic:
new_content[t][i,j]=self.content[wrap(t+i+j)][0] new_content[t][i, j] = self.content[wrap(t+i+j)][0]
elif (t+i+j)>=self.T: elif (t+i+j) >= self.T:
new_content[t]=None new_content[t]=None
else: 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) return Corr(new_content)
def roll(self, dt): def roll(self, dt):
"""Periodically shift the correlator by dt timeslices """Periodically shift the correlator by dt timeslices
@ -701,8 +684,8 @@ class Corr:
content_string += "Description: " + self.tag + "\n" content_string += "Description: " + self.tag + "\n"
if self.N!=1: if self.N!=1:
return content_string 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 # 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. # be a wall of numbers.
if range[1]: if range[1]:
@ -929,7 +912,7 @@ class Corr:
if isinstance(obs_OR_cobs, CObs): if isinstance(obs_OR_cobs, CObs):
return obs_OR_cobs.imag return obs_OR_cobs.imag
else: 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) 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]) reference_sorting=np.array(vec_set[ts])
N=reference_sorting.shape[0] N=reference_sorting.shape[0]
sorted_vec_set=[] 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): for k in range(N):
new_sorting=reference_sorting.copy() new_sorting=reference_sorting.copy()
new_sorting[perm[k],:]=vec_set[t][k] 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: if current_score>best_score:
best_score=current_score best_score=current_score
best_perm=perm 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: if len(lst) == 1:
return [lst] return [lst]
l = [] l = []
@ -993,7 +976,7 @@ def permutation(lst): #Shamelessly copied
return l 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_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=[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] sp_vecs=[v/np.sqrt((v.T@G0@v)) for v in sp_vecs]