lattice/pyerrors
1
0
Fork 0
mirror of https://github.com/fjosw/pyerrors.git synced 2025-03-15 14:50:25 +01:00
Code Issues Projects Releases Packages Wiki Activity

linting tested again

Browse source
This commit is contained in:
JanNeuendorf 2022-01-18 16:46:14 +01:00
parent 2f11d0d30b
commit 909ef85ff8
1 changed files with 65 additions and 87 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

152
pyerrors/correlators.py
View file

@ -116,24 +116,24 @@ class Corr:
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 vector_l.shape == vector_r.shape == (self.N,):
raise Exception("Vectors are of wrong shape!")
if normalize:
vector_l, vector_r = vector_l / np.sqrt((vector_l @ vector_l)), vector_r / np.sqrt(vector_r @ vector_r)
#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!")
# 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!")
newcontent = [None if (item is None) else np.asarray([vector_l.T @ item @ vector_r]) for item in self.content]
else:
#There are no checks here yet. There are so many possible scenarios, where this can go wrong.
# There are no checks here yet. There are so many possible scenarios, where this can go wrong.
if normalize:
for t in range(self.T):
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])
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)]
return Corr(newcontent)
@ -210,7 +210,7 @@ class Corr:
if self.N == 1:
raise Exception("Trying to symmetrize a smearing matrix, that already has N=1.")
# 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
# 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):
@ -226,7 +226,7 @@ class Corr:
sp_vecs = GEVP_solver(Gt, G0)
sp_vec = sp_vecs[state]
return sp_vec
if return_list:
if return_list:
all_vecs = []
for t in range(self.T):
try:
@ -242,12 +242,12 @@ class Corr:
all_vecs.append(sp_vec)
else:
all_vecs.append(sp_vecs)
except "Failure to solve for one timeslice": # 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]
return all_vecs
def Eigenvalue(self, t0, state=1):
@ -269,36 +269,36 @@ class Corr:
newcontent.append(eigenvalue)
return Corr(newcontent)
def Hankel(self, N, periodic=False):
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!")
array = np.empty([N, N], dtype="object")
new_content = []
new_content = []
for t in range(self.T):
new_content.append(array.copy())
def wrap(i):
if i >= self.T:
return i-self.T
return i
new_content.append(array.copy())
for t in range(self.T):
def wrap(i):
if i >= self.T:
return i - self.T
return i
for t in range(self.T):
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]=None
else:
new_content[t][i, j] = self.content[t+i+j][0]
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]
return Corr(new_content)
def roll(self, dt):
@ -313,7 +313,7 @@ class Corr:
def reverse(self):
"""Reverse the time ordering of the Corr"""
return Corr(self.content[::-1])
return Corr(self.content[:: -1])
def correlate(self, partner):
"""Correlate the correlator with another correlator or Obs
@ -335,7 +335,7 @@ class Corr:
new_content.append(None)
else:
new_content.append(np.array([correlate(o, partner.content[x0][0]) for o in t_slice]))
elif isinstance(partner, Obs): # Should this include CObs?
elif isinstance(partner, Obs): # Should this include CObs?
new_content.append(np.array([correlate(o, partner) for o in t_slice]))
else:
raise Exception("Can only correlate with an Obs or a Corr.")
@ -676,18 +676,15 @@ class Corr:
def __repr__(self, range=[0, None]):
content_string = ""
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
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
if self.tag is not None:
content_string += "Description: " + self.tag + "\n"
if self.N!=1:
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
# 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]:
range[1] += 1
content_string += 'x0/a\tCorr(x0/a)\n------------------\n'
@ -782,7 +779,7 @@ class Corr:
if y.value == 0:
raise Exception('Division by zero will return undefined correlator')
if isinstance(y, CObs):
if y.is_zero():
if y.is_zero():
raise Exception('Division by zero will return undefined correlator')
newcontent = []
@ -837,11 +834,11 @@ class Corr:
return Corr(newcontent, prange=self.prange)
def _apply_func_to_corr(self, func):
newcontent = [None if (item is None ) else func(item) for item in self.content]
newcontent = [None if (item is None) else func(item) for item in self.content]
for t in range(self.T):
if newcontent[t] is None:
continue
if np.isnan(np.sum(newcontent[t]).value):
if np.isnan(np.sum(newcontent[t]).value):
newcontent[t] = None
if all([item is None for item in newcontent]):
raise Exception('Operation returns undefined correlator')
@ -897,87 +894,68 @@ class Corr:
return (self / y) ** (-1)
@property
def real(self):
def real(self):
def return_real(obs_OR_cobs):
if isinstance(obs_OR_cobs, CObs):
return obs_OR_cobs.real
return obs_OR_cobs.real
else:
return obs_OR_cobs
return obs_OR_cobs
return self._apply_func_to_corr(return_real)
@property
def imag(self):
def imag(self):
def return_imag(obs_OR_cobs):
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)
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=[]
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 = []
for t in range(len(vec_set)):
if vec_set[t] is None:
sorted_vec_set.append(None)
elif not t==ts:
perms=permutation([i for i in range(N)])
best_score=0
elif not t == ts:
perms = permutation([i for i in range(N)])
best_score = 0
for perm in perms:
current_score=1
current_score = 1
for k in range(N):
new_sorting=reference_sorting.copy()
new_sorting[perm[k],:]=vec_set[t][k]
new_sorting = reference_sorting.copy()
new_sorting[perm[k], :] = vec_set[t][k]
current_score *= abs(np.linalg.det(new_sorting))
if current_score>best_score:
best_score=current_score
best_perm=perm
#print("best perm", best_perm)
if current_score > best_score:
best_score = current_score
best_perm = perm
# print("best perm", best_perm)
sorted_vec_set.append([vec_set[t][k] for k in best_perm])
else:
sorted_vec_set.append(vec_set[t])
sorted_vec_set.append(vec_set[t])
return sorted_vec_set
def permutation(lst): # Shamelessly copied
def permutation(lst): # Shamelessly copied
if len(lst) == 1:
return [lst]
l = []
ll = []
for i in range(len(lst)):
m = lst[i]
remLst = lst[:i] + lst[i+1:]
remLst = lst[:i] + lst[i + 1:]
# Generating all permutations where m is first
for p in permutation(remLst):
l.append([m] + p)
return l
ll.append([m] + p)
return ll
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]
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