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

Merge pull request #77 from JanNeuendorf/develop

Browse source 
Changes to Eigenvalue
This commit is contained in:
Fabian Joswig 2022-02-22 21:26:21 +00:00 • committed by GitHub
commit 7830ba7b6c
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
1 changed files with 7 additions and 22 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

29
pyerrors/correlators.py
View file

@ -7,7 +7,6 @@ import scipy.linalg
from .obs import Obs, reweight, correlate, CObs
from .misc import dump_object, _assert_equal_properties
from .fits import least_squares
from .linalg import eigh, inv, cholesky
from .roots import find_root
@ -266,10 +265,11 @@ class Corr:
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")
symmetric_corr = self.matrix_symmetric()
for i in range(self.N):
for j in range(self.N):
G0[i, j] = self.content[t0][i, j].value
Gt[i, j] = self.content[ts][i, j].value
G0[i, j] = symmetric_corr.content[t0][i, j].value
Gt[i, j] = symmetric_corr[ts][i, j].value
sp_vecs = _GEVP_solver(Gt, G0)
sp_vec = sp_vecs[state]
@ -301,24 +301,9 @@ class Corr:
return all_vecs
def Eigenvalue(self, t0, state=1):
G = self.matrix_symmetric()
G0 = G.content[t0]
L = cholesky(G0)
Li = inv(L)
LT = L.T
LTi = inv(LT)
newcontent = []
for t in range(self.T):
if self.content[t] is None:
newcontent.append(None)
else:
Gt = G.content[t]
M = Li @ Gt @ LTi
eigenvalues = eigh(M)[0]
eigenvalue = eigenvalues[-state]
newcontent.append(eigenvalue)
return Corr(newcontent)
def Eigenvalue(self, t0, ts=None, state=0, sorted_list=None):
vec = self.GEVP(t0, ts=ts, state=state, sorted_list=sorted_list)
return self.projected(vec)
def Hankel(self, N, periodic=False):
"""Constructs an NxN Hankel matrix
@ -1096,7 +1081,7 @@ def _sort_vectors(vec_set, ts):
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.eigh(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