From ba054fa11caed05f8b12e94368b6987492c113b9 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 16 May 2022 11:35:23 +0100 Subject: [PATCH 01/13] refactor: correlators._solve_GEVP simplified and optimized, test added. --- pyerrors/correlators.py | 8 +++----- tests/correlators_test.py | 15 +++++++++++++++ 2 files changed, 18 insertions(+), 5 deletions(-) diff --git a/pyerrors/correlators.py b/pyerrors/correlators.py index 233c7499..19d0c334 100644 --- a/pyerrors/correlators.py +++ b/pyerrors/correlators.py @@ -1219,8 +1219,6 @@ def _sort_vectors(vec_set, ts): return sorted_vec_set -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.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 +def _GEVP_solver(Gt, G0): + """Helper function for solving the GEVP and sorting the eigenvectors.""" + return scipy.linalg.eigh(Gt, G0)[1].T[::-1] diff --git a/tests/correlators_test.py b/tests/correlators_test.py index 5242339f..c520917a 100644 --- a/tests/correlators_test.py +++ b/tests/correlators_test.py @@ -1,5 +1,6 @@ import os import numpy as np +import scipy import pyerrors as pe import pytest @@ -290,6 +291,20 @@ def test_matrix_symmetric(): assert np.all([np.all(o == o.T) for o in sym_corr_mat]) +def test_GEVP_solver(): + + mat1 = np.random.rand(15, 15) + mat2 = np.random.rand(15, 15) + mat1 = mat1 @ mat1.T + mat2 = mat2 @ mat2.T + + sp_val, sp_vecs = scipy.linalg.eigh(mat1, mat2) + 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 @ mat2 @ v)) for v in sp_vecs] + + assert np.allclose(sp_vecs, pe.correlators._GEVP_solver(mat1, mat2), atol=1e-14) + + def test_hankel(): corr_content = [] for t in range(8): From 92b19cba9f852c88e797c3df44c983809401ad0c Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 16 May 2022 11:44:02 +0100 Subject: [PATCH 02/13] !feat: GEVP now returns all eigenvectors instead of just the ones for the specified state. --- pyerrors/correlators.py | 18 +++++++----------- tests/correlators_test.py | 7 +++---- 2 files changed, 10 insertions(+), 15 deletions(-) diff --git a/pyerrors/correlators.py b/pyerrors/correlators.py index 19d0c334..cf236d8d 100644 --- a/pyerrors/correlators.py +++ b/pyerrors/correlators.py @@ -241,8 +241,8 @@ class Corr: if self.N == 1: raise Exception("Trying to symmetrize a correlator matrix, that already has N=1.") - def GEVP(self, t0, ts=None, state=0, sorted_list="Eigenvalue"): - """Solve the general eigenvalue problem on the current correlator + def GEVP(self, t0, ts=None, sorted_list="Eigenvalue"): + """Solve the generalized eigenvalue problem on the current correlator and returns the corresponding eigenvectors. Parameters ---------- @@ -251,8 +251,6 @@ class Corr: ts : int fixed time G(t_s)v= lambda G(t_0)v if return_list=False If return_list=True and sorting=Eigenvector it gives a reference point for the sorting method. - state : int - The state one is interested in, ordered by energy. The lowest state is zero. sorted_list : string if this argument is set, a list of vectors (len=self.T) is returned. If it is left as None, only one vector is returned. "Eigenvalue" - The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice. @@ -278,7 +276,7 @@ class Corr: Gt[i, j] = symmetric_corr[ts][i, j].value sp_vecs = _GEVP_solver(Gt, G0) - sp_vec = sp_vecs[state] + sp_vec = [sp_vecs[s] for s in range(self.N)] return sp_vec elif sorted_list in ["Eigenvalue", "Eigenvector"]: if sorted_list == "Eigenvalue" and ts is not None: @@ -294,7 +292,7 @@ class Corr: sp_vecs = _GEVP_solver(Gt, G0) if sorted_list == "Eigenvalue": - sp_vec = sp_vecs[state] + sp_vec = [sp_vecs[s] for s in range(self.N)] all_vecs.append(sp_vec) else: all_vecs.append(sp_vecs) @@ -304,7 +302,7 @@ class Corr: if (ts is None): raise Exception("ts is required for the Eigenvector sorting method.") all_vecs = _sort_vectors(all_vecs, ts) - all_vecs = [a[state] if a is not None else None for a in all_vecs] + all_vecs = [[a[s] if a is not None else None for a in all_vecs] for s in range(self.N)] else: raise Exception("Unkown value for 'sorted_list'.") @@ -328,7 +326,7 @@ class Corr: "Eigenvector" - Use the method described in arXiv:2004.10472 [hep-lat] to find the set of v(t) belonging to the state. The reference state is identified by its eigenvalue at t=ts """ - vec = self.GEVP(t0, ts=ts, state=state, sorted_list=sorted_list) + vec = self.GEVP(t0, ts=ts, sorted_list=sorted_list)[state] return self.projected(vec) def Hankel(self, N, periodic=False): @@ -1176,9 +1174,7 @@ class Corr: if basematrix.N != self.N: raise Exception('basematrix and targetmatrix have to be of the same size.') - evecs = [] - for i in range(Ntrunc): - evecs.append(basematrix.GEVP(t0proj, tproj, state=i, sorted_list=None)) + evecs = basematrix.GEVP(t0proj, tproj, sorted_list=None)[:Ntrunc] tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object) rmat = [] diff --git a/tests/correlators_test.py b/tests/correlators_test.py index c520917a..23c99db6 100644 --- a/tests/correlators_test.py +++ b/tests/correlators_test.py @@ -231,11 +231,10 @@ def test_matrix_corr(): corr_mat = pe.Corr(np.array([[corr_aa, corr_ab], [corr_ab, corr_aa]])) corr_mat.item(0, 0) - vec_0 = corr_mat.GEVP(0, 1, sorted_list=None) - vec_1 = corr_mat.GEVP(0, 1, state=1, sorted_list=None) + vecs = corr_mat.GEVP(0, 1, sorted_list=None) - corr_0 = corr_mat.projected(vec_0) - corr_1 = corr_mat.projected(vec_1) + corr_0 = corr_mat.projected(vecs[0]) + corr_1 = corr_mat.projected(vecs[0]) assert np.all([o == 0 for o in corr_0 - corr_aa]) assert np.all([o == 0 for o in corr_1 - corr_aa]) From c00c21ee8629a9b444d905a0a37c2dfb5079a466 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 16 May 2022 11:45:10 +0100 Subject: [PATCH 03/13] !refactor: argument sorted_list of Corr.GEVP renamed to sort. --- pyerrors/correlators.py | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/pyerrors/correlators.py b/pyerrors/correlators.py index cf236d8d..23c8b57b 100644 --- a/pyerrors/correlators.py +++ b/pyerrors/correlators.py @@ -241,7 +241,7 @@ class Corr: if self.N == 1: raise Exception("Trying to symmetrize a correlator matrix, that already has N=1.") - def GEVP(self, t0, ts=None, sorted_list="Eigenvalue"): + def GEVP(self, t0, ts=None, sort="Eigenvalue"): """Solve the generalized eigenvalue problem on the current correlator and returns the corresponding eigenvectors. Parameters @@ -251,7 +251,7 @@ class Corr: ts : int fixed time G(t_s)v= lambda G(t_0)v if return_list=False If return_list=True and sorting=Eigenvector it gives a reference point for the sorting method. - sorted_list : string + sort : string if this argument is set, a list of vectors (len=self.T) is returned. If it is left as None, only one vector is returned. "Eigenvalue" - The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice. "Eigenvector" - Use the method described in arXiv:2004.10472 [hep-lat] to find the set of v(t) belonging to the state. @@ -262,9 +262,9 @@ class Corr: raise Exception("GEVP methods only works on correlator matrices and not single correlators.") symmetric_corr = self.matrix_symmetric() - if sorted_list is None: + if sort is None: if (ts is None): - raise Exception("ts is required if sorted_list=None.") + raise Exception("ts is required if sort=None.") if (ts <= t0): raise Exception("ts has to be larger than t0.") if (self.content[t0] is None) or (self.content[ts] is None): @@ -278,8 +278,8 @@ class Corr: sp_vecs = _GEVP_solver(Gt, G0) sp_vec = [sp_vecs[s] for s in range(self.N)] return sp_vec - elif sorted_list in ["Eigenvalue", "Eigenvector"]: - if sorted_list == "Eigenvalue" and ts is not None: + elif sort in ["Eigenvalue", "Eigenvector"]: + if sort == "Eigenvalue" and ts is not None: warnings.warn("ts has no effect when sorting by eigenvalue is chosen.", RuntimeWarning) all_vecs = [None] * (t0 + 1) for t in range(t0 + 1, self.T): @@ -291,24 +291,24 @@ class Corr: Gt[i, j] = symmetric_corr[t][i, j].value sp_vecs = _GEVP_solver(Gt, G0) - if sorted_list == "Eigenvalue": + if sort == "Eigenvalue": sp_vec = [sp_vecs[s] for s in range(self.N)] all_vecs.append(sp_vec) else: all_vecs.append(sp_vecs) except Exception: all_vecs.append(None) - if sorted_list == "Eigenvector": + if sort == "Eigenvector": if (ts is None): raise Exception("ts is required for the Eigenvector sorting method.") all_vecs = _sort_vectors(all_vecs, ts) all_vecs = [[a[s] if a is not None else None for a in all_vecs] for s in range(self.N)] else: - raise Exception("Unkown value for 'sorted_list'.") + raise Exception("Unkown value for 'sort'.") return all_vecs - def Eigenvalue(self, t0, ts=None, state=0, sorted_list=None): + def Eigenvalue(self, t0, ts=None, state=0, sort=None): """Determines the eigenvalue of the GEVP by solving and projecting the correlator Parameters @@ -320,13 +320,13 @@ class Corr: If return_list=True and sorting=Eigenvector it gives a reference point for the sorting method. state : int The state one is interested in ordered by energy. The lowest state is zero. - sorted_list : string + sort : string if this argument is set, a list of vectors (len=self.T) is returned. If it is left as None, only one vector is returned. "Eigenvalue" - The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice. "Eigenvector" - Use the method described in arXiv:2004.10472 [hep-lat] to find the set of v(t) belonging to the state. The reference state is identified by its eigenvalue at t=ts """ - vec = self.GEVP(t0, ts=ts, sorted_list=sorted_list)[state] + vec = self.GEVP(t0, ts=ts, sort=sort)[state] return self.projected(vec) def Hankel(self, N, periodic=False): @@ -1174,7 +1174,7 @@ class Corr: if basematrix.N != self.N: raise Exception('basematrix and targetmatrix have to be of the same size.') - evecs = basematrix.GEVP(t0proj, tproj, sorted_list=None)[:Ntrunc] + evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc] tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object) rmat = [] From cfc21a2a5e9d0662e1092ee1fc39a9462d0824a5 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 16 May 2022 12:05:33 +0100 Subject: [PATCH 04/13] docs: GEVP example updated to new syntax. --- examples/06_gevp.ipynb | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/examples/06_gevp.ipynb b/examples/06_gevp.ipynb index 9bb9ed49..e969008c 100644 --- a/examples/06_gevp.ipynb +++ b/examples/06_gevp.ipynb @@ -168,9 +168,9 @@ } ], "source": [ - "vec = matrix_V1V1.GEVP(t0=3, ts=6 ,state=0, sorted_list=None)\n", - "assert len(vec) == matrix_V1V1.N\n", - "print(vec)" + "vecs = matrix_V1V1.GEVP(t0=3, ts=6, sort=None)\n", + "assert len(vecs[0]) == matrix_V1V1.N\n", + "print(vecs[0])" ] }, { @@ -202,7 +202,7 @@ } ], "source": [ - "matrix_V1V1.projected(vec).m_eff().show(comp=single_smearing.m_eff(), auto_gamma=True)" + "matrix_V1V1.projected(vecs[0]).m_eff().show(comp=single_smearing.m_eff(), auto_gamma=True)" ] }, { @@ -266,7 +266,7 @@ "\n", "# We then solve the GEVP to get eigenvectors corresponding to the ground state.\n", "\n", - "vec_ground = matrix_VnVn.GEVP(t0=3, ts=6, state=0, sorted_list=None)\n", + "vec_ground = matrix_VnVn.GEVP(t0=3, ts=6, sort=None)[0]\n", "\n", "# Now we project the matrix-correlators to get new correlators belonging to the ground state.\n", "\n", From 410d7618a0f343d60b465216c0e777a6b39847d3 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 16 May 2022 12:07:36 +0100 Subject: [PATCH 05/13] tests: GEVP tests updated to new syntax. --- tests/correlators_test.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/correlators_test.py b/tests/correlators_test.py index 23c99db6..4d6e9b61 100644 --- a/tests/correlators_test.py +++ b/tests/correlators_test.py @@ -231,7 +231,7 @@ def test_matrix_corr(): corr_mat = pe.Corr(np.array([[corr_aa, corr_ab], [corr_ab, corr_aa]])) corr_mat.item(0, 0) - vecs = corr_mat.GEVP(0, 1, sorted_list=None) + vecs = corr_mat.GEVP(0, 1, sort=None) corr_0 = corr_mat.projected(vecs[0]) corr_1 = corr_mat.projected(vecs[0]) @@ -239,13 +239,13 @@ def test_matrix_corr(): assert np.all([o == 0 for o in corr_0 - corr_aa]) assert np.all([o == 0 for o in corr_1 - corr_aa]) - corr_mat.GEVP(0, sorted_list="Eigenvalue") - corr_mat.GEVP(0, 1, sorted_list="Eigenvector") + corr_mat.GEVP(0, sort="Eigenvalue") + corr_mat.GEVP(0, 1, sort="Eigenvector") corr_mat.matrix_symmetric() with pytest.warns(RuntimeWarning): - corr_mat.GEVP(0, 1, sorted_list="Eigenvalue") + corr_mat.GEVP(0, 1, sort="Eigenvalue") with pytest.raises(Exception): corr_mat.plottable() From 4c06f9886d9f724310c6d75137702141b34ce287 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 16 May 2022 13:22:42 +0100 Subject: [PATCH 06/13] fix: GEVP sorted vectors fixed and simplified. --- pyerrors/correlators.py | 14 ++++---------- 1 file changed, 4 insertions(+), 10 deletions(-) diff --git a/pyerrors/correlators.py b/pyerrors/correlators.py index 23c8b57b..60c72ec0 100644 --- a/pyerrors/correlators.py +++ b/pyerrors/correlators.py @@ -276,8 +276,8 @@ class Corr: Gt[i, j] = symmetric_corr[ts][i, j].value sp_vecs = _GEVP_solver(Gt, G0) - sp_vec = [sp_vecs[s] for s in range(self.N)] - return sp_vec + return sp_vecs + elif sort in ["Eigenvalue", "Eigenvector"]: if sort == "Eigenvalue" and ts is not None: warnings.warn("ts has no effect when sorting by eigenvalue is chosen.", RuntimeWarning) @@ -290,23 +290,17 @@ class Corr: G0[i, j] = symmetric_corr[t0][i, j].value Gt[i, j] = symmetric_corr[t][i, j].value - sp_vecs = _GEVP_solver(Gt, G0) - if sort == "Eigenvalue": - sp_vec = [sp_vecs[s] for s in range(self.N)] - all_vecs.append(sp_vec) - else: - all_vecs.append(sp_vecs) + all_vecs.append(_GEVP_solver(Gt, G0)) except Exception: all_vecs.append(None) if sort == "Eigenvector": if (ts is None): raise Exception("ts is required for the Eigenvector sorting method.") all_vecs = _sort_vectors(all_vecs, ts) - all_vecs = [[a[s] if a is not None else None for a in all_vecs] for s in range(self.N)] else: raise Exception("Unkown value for 'sort'.") - return all_vecs + return [[v[s] if v is not None else None for v in all_vecs] for s in range(self.N)] def Eigenvalue(self, t0, ts=None, state=0, sort=None): """Determines the eigenvalue of the GEVP by solving and projecting the correlator From eb8090a90cc2390433a75fbb03961c5707c98276 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 16 May 2022 13:28:10 +0100 Subject: [PATCH 07/13] tests: GEVP tests extended to all three sorting variants. --- tests/correlators_test.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/tests/correlators_test.py b/tests/correlators_test.py index 4d6e9b61..18c578f4 100644 --- a/tests/correlators_test.py +++ b/tests/correlators_test.py @@ -231,22 +231,23 @@ def test_matrix_corr(): corr_mat = pe.Corr(np.array([[corr_aa, corr_ab], [corr_ab, corr_aa]])) corr_mat.item(0, 0) - vecs = corr_mat.GEVP(0, 1, sort=None) + for (ts, sort) in zip([None, 1, 1], ["Eigenvalue", "Eigenvector", None]): + vecs = corr_mat.GEVP(0, ts=ts, sort=sort) - corr_0 = corr_mat.projected(vecs[0]) - corr_1 = corr_mat.projected(vecs[0]) + corr_0 = corr_mat.projected(vecs[0]) + corr_1 = corr_mat.projected(vecs[1]) assert np.all([o == 0 for o in corr_0 - corr_aa]) assert np.all([o == 0 for o in corr_1 - corr_aa]) - corr_mat.GEVP(0, sort="Eigenvalue") - corr_mat.GEVP(0, 1, sort="Eigenvector") - corr_mat.matrix_symmetric() with pytest.warns(RuntimeWarning): corr_mat.GEVP(0, 1, sort="Eigenvalue") + with pytest.raises(Exception): + corr_mat.GEVP(0, 1, sort="This sorting method does not exist.") + with pytest.raises(Exception): corr_mat.plottable() From 2136958fbc54426fb14bfc97d3c60a8deb31d872 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 16 May 2022 13:42:12 +0100 Subject: [PATCH 08/13] feat: sorted_list argument readded to Corr.GEVP with a deprecation warning. --- pyerrors/correlators.py | 6 +++++- tests/correlators_test.py | 3 +++ 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/pyerrors/correlators.py b/pyerrors/correlators.py index 60c72ec0..ef50b743 100644 --- a/pyerrors/correlators.py +++ b/pyerrors/correlators.py @@ -241,7 +241,7 @@ class Corr: if self.N == 1: raise Exception("Trying to symmetrize a correlator matrix, that already has N=1.") - def GEVP(self, t0, ts=None, sort="Eigenvalue"): + def GEVP(self, t0, ts=None, sort="Eigenvalue", **kwargs): """Solve the generalized eigenvalue problem on the current correlator and returns the corresponding eigenvectors. Parameters @@ -261,6 +261,10 @@ class Corr: if self.N == 1: raise Exception("GEVP methods only works on correlator matrices and not single correlators.") + if "sorted_list" in kwargs: + sort = kwargs.get("sorted_list") + warnings.warn("Argument 'sorted_list' is deprecated, use 'sort' instead.", DeprecationWarning) + symmetric_corr = self.matrix_symmetric() if sort is None: if (ts is None): diff --git a/tests/correlators_test.py b/tests/correlators_test.py index 18c578f4..8cc6e520 100644 --- a/tests/correlators_test.py +++ b/tests/correlators_test.py @@ -245,6 +245,9 @@ def test_matrix_corr(): with pytest.warns(RuntimeWarning): corr_mat.GEVP(0, 1, sort="Eigenvalue") + with pytest.warns(DeprecationWarning): + corr_mat.GEVP(0, sorted_list="Eigenvalue") + with pytest.raises(Exception): corr_mat.GEVP(0, 1, sort="This sorting method does not exist.") From 9c17b8e7192ad9d63794ce53ea0f60cb445ab2a2 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 16 May 2022 14:13:33 +0100 Subject: [PATCH 09/13] feat: argument state readded to Corr.GEVP as deprecated kwarg. Documentation of GEVP and _solve_GEVP extended. --- pyerrors/correlators.py | 42 ++++++++++++++++++++++++++--------------- 1 file changed, 27 insertions(+), 15 deletions(-) diff --git a/pyerrors/correlators.py b/pyerrors/correlators.py index ef50b743..5784f186 100644 --- a/pyerrors/correlators.py +++ b/pyerrors/correlators.py @@ -242,28 +242,31 @@ class Corr: raise Exception("Trying to symmetrize a correlator matrix, that already has N=1.") def GEVP(self, t0, ts=None, sort="Eigenvalue", **kwargs): - """Solve the generalized eigenvalue problem on the current correlator and returns the corresponding eigenvectors. + r'''Solve the generalized eigenvalue problem on the correlator matrix and returns the corresponding eigenvectors. + + The eigenvectors are sorted according to the descending eigenvalues, the zeroth eigenvector(s) correspond to the + largest eigenvalue(s). Parameters ---------- t0 : int - The time t0 for G(t)v= lambda G(t_0)v + The time t0 for the right hand side of the GEVP according to $G(t)v_i=\lambda_i G(t_0)v_i$ ts : int - fixed time G(t_s)v= lambda G(t_0)v if return_list=False - If return_list=True and sorting=Eigenvector it gives a reference point for the sorting method. + fixed time $G(t_s)v_i=\lambda_i G(t_0)v_i$ if sort=None. + If sort="Eigenvector" it gives a reference point for the sorting method. sort : string - if this argument is set, a list of vectors (len=self.T) is returned. If it is left as None, only one vector is returned. - "Eigenvalue" - The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice. - "Eigenvector" - Use the method described in arXiv:2004.10472 [hep-lat] to find the set of v(t) belonging to the state. - The reference state is identified by its eigenvalue at t=ts - """ + If this argument is set, a list of self.T vectors per state is returned. If it is set to None, only one vector is returned. + - "Eigenvalue": The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice. + - "Eigenvector": Use the method described in arXiv:2004.10472 to find the set of v(t) belonging to the state. + The reference state is identified by its eigenvalue at $t=t_s$. + ''' if self.N == 1: raise Exception("GEVP methods only works on correlator matrices and not single correlators.") if "sorted_list" in kwargs: - sort = kwargs.get("sorted_list") warnings.warn("Argument 'sorted_list' is deprecated, use 'sort' instead.", DeprecationWarning) + sort = kwargs.get("sorted_list") symmetric_corr = self.matrix_symmetric() if sort is None: @@ -279,8 +282,7 @@ class Corr: G0[i, j] = symmetric_corr[t0][i, j].value Gt[i, j] = symmetric_corr[ts][i, j].value - sp_vecs = _GEVP_solver(Gt, G0) - return sp_vecs + reordered_vecs = _GEVP_solver(Gt, G0) elif sort in ["Eigenvalue", "Eigenvector"]: if sort == "Eigenvalue" and ts is not None: @@ -301,10 +303,16 @@ class Corr: if (ts is None): raise Exception("ts is required for the Eigenvector sorting method.") all_vecs = _sort_vectors(all_vecs, ts) + + reordered_vecs = [[v[s] if v is not None else None for v in all_vecs] for s in range(self.N)] else: raise Exception("Unkown value for 'sort'.") - return [[v[s] if v is not None else None for v in all_vecs] for s in range(self.N)] + if "state" in kwargs: + warnings.warn("Argument 'state' is deprecated, use slicing instead.", DeprecationWarning) + return reordered_vecs[kwargs.get("state")] + else: + return reordered_vecs def Eigenvalue(self, t0, ts=None, state=0, sort=None): """Determines the eigenvalue of the GEVP by solving and projecting the correlator @@ -1214,5 +1222,9 @@ def _sort_vectors(vec_set, ts): def _GEVP_solver(Gt, G0): - """Helper function for solving the GEVP and sorting the eigenvectors.""" - return scipy.linalg.eigh(Gt, G0)[1].T[::-1] + """Helper function for solving the GEVP and sorting the eigenvectors. + + The helper function assumes that both provided matrices are symmetric and + only processes the lower triangular part of both matrices. In case the matrices + are not symmetric the upper triangular parts are effectively discarded.""" + return scipy.linalg.eigh(Gt, G0, lower=True)[1].T[::-1] From 43118bb67da9d788347f51d07958f82e53c9d067 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 16 May 2022 14:19:15 +0100 Subject: [PATCH 10/13] docs: docstring of Corr.Eigenvalue simplified. --- pyerrors/correlators.py | 14 +++----------- 1 file changed, 3 insertions(+), 11 deletions(-) diff --git a/pyerrors/correlators.py b/pyerrors/correlators.py index 5784f186..f6de50e5 100644 --- a/pyerrors/correlators.py +++ b/pyerrors/correlators.py @@ -314,23 +314,15 @@ class Corr: else: return reordered_vecs - def Eigenvalue(self, t0, ts=None, state=0, sort=None): + def Eigenvalue(self, t0, ts=None, state=0, sort="Eigenvalue"): """Determines the eigenvalue of the GEVP by solving and projecting the correlator Parameters ---------- - t0 : int - The time t0 for G(t)v= lambda G(t_0)v - ts : int - fixed time G(t_s)v= lambda G(t_0)v if return_list=False - If return_list=True and sorting=Eigenvector it gives a reference point for the sorting method. state : int The state one is interested in ordered by energy. The lowest state is zero. - sort : string - if this argument is set, a list of vectors (len=self.T) is returned. If it is left as None, only one vector is returned. - "Eigenvalue" - The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice. - "Eigenvector" - Use the method described in arXiv:2004.10472 [hep-lat] to find the set of v(t) belonging to the state. - The reference state is identified by its eigenvalue at t=ts + + All other parameters are identical to the ones of Corr.GEVP. """ vec = self.GEVP(t0, ts=ts, sort=sort)[state] return self.projected(vec) From 98ce553521794e9d8bc5ad2b80112a6b80c4b343 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Wed, 18 May 2022 10:02:56 +0100 Subject: [PATCH 11/13] fix: exception for ts<=t0 generalized, tests added. --- pyerrors/correlators.py | 5 +++-- tests/correlators_test.py | 24 ++++++++++++++++++++++++ 2 files changed, 27 insertions(+), 2 deletions(-) diff --git a/pyerrors/correlators.py b/pyerrors/correlators.py index f6de50e5..a8f4def6 100644 --- a/pyerrors/correlators.py +++ b/pyerrors/correlators.py @@ -263,6 +263,9 @@ class Corr: if self.N == 1: raise Exception("GEVP methods only works on correlator matrices and not single correlators.") + if ts is not None: + if (ts <= t0): + raise Exception("ts has to be larger than t0.") if "sorted_list" in kwargs: warnings.warn("Argument 'sorted_list' is deprecated, use 'sort' instead.", DeprecationWarning) @@ -272,8 +275,6 @@ class Corr: if sort is None: if (ts is None): raise Exception("ts is required if sort=None.") - if (ts <= t0): - raise Exception("ts has to be larger than t0.") 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") diff --git a/tests/correlators_test.py b/tests/correlators_test.py index 8cc6e520..0c1d089b 100644 --- a/tests/correlators_test.py +++ b/tests/correlators_test.py @@ -242,12 +242,36 @@ def test_matrix_corr(): corr_mat.matrix_symmetric() + +def test_GEVP_warnings(): + corr_aa = _gen_corr(1) + corr_ab = 0.5 * corr_aa + + corr_mat = pe.Corr(np.array([[corr_aa, corr_ab], [corr_ab, corr_aa]])) + corr_mat.item(0, 0) + with pytest.warns(RuntimeWarning): corr_mat.GEVP(0, 1, sort="Eigenvalue") with pytest.warns(DeprecationWarning): corr_mat.GEVP(0, sorted_list="Eigenvalue") + with pytest.warns(DeprecationWarning): + corr_mat.GEVP(0, state=0) + +def test_GEVP_exceptions(): + corr_aa = _gen_corr(1) + corr_ab = 0.5 * corr_aa + + corr_mat = pe.Corr(np.array([[corr_aa, corr_ab], [corr_ab, corr_aa]])) + corr_mat.item(0, 0) + + with pytest.raises(Exception): + corr_mat.GEVP(0, 0, sort=None) + + with pytest.raises(Exception): + corr_mat.GEVP(1, 0, sort="Eigenvector") + with pytest.raises(Exception): corr_mat.GEVP(0, 1, sort="This sorting method does not exist.") From 088a7da7e60814c3c2aa84a1635a0aed3c800ff9 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Wed, 18 May 2022 10:09:23 +0100 Subject: [PATCH 12/13] feat: Deprecation warning for Corr.GEVP kwarg state removed, documentation extended. --- pyerrors/correlators.py | 12 ++++++++++-- tests/correlators_test.py | 3 --- 2 files changed, 10 insertions(+), 5 deletions(-) diff --git a/pyerrors/correlators.py b/pyerrors/correlators.py index a8f4def6..ab4f6f50 100644 --- a/pyerrors/correlators.py +++ b/pyerrors/correlators.py @@ -245,7 +245,11 @@ class Corr: r'''Solve the generalized eigenvalue problem on the correlator matrix and returns the corresponding eigenvectors. The eigenvectors are sorted according to the descending eigenvalues, the zeroth eigenvector(s) correspond to the - largest eigenvalue(s). + largest eigenvalue(s). The eigenvector(s) for the individual states can be accessed via slicing + ```python + C.GEVP(t0=2)[0] # Ground state vector(s) + C.GEVP(t0=2)[:3] # Vectors for the lowest three states + ``` Parameters ---------- @@ -259,6 +263,11 @@ class Corr: - "Eigenvalue": The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice. - "Eigenvector": Use the method described in arXiv:2004.10472 to find the set of v(t) belonging to the state. The reference state is identified by its eigenvalue at $t=t_s$. + + Other Parameters + ---------------- + state : int + Returns only the vector(s) for a specified state. The lowest state is zero. ''' if self.N == 1: @@ -310,7 +319,6 @@ class Corr: raise Exception("Unkown value for 'sort'.") if "state" in kwargs: - warnings.warn("Argument 'state' is deprecated, use slicing instead.", DeprecationWarning) return reordered_vecs[kwargs.get("state")] else: return reordered_vecs diff --git a/tests/correlators_test.py b/tests/correlators_test.py index 0c1d089b..90dac90a 100644 --- a/tests/correlators_test.py +++ b/tests/correlators_test.py @@ -256,9 +256,6 @@ def test_GEVP_warnings(): with pytest.warns(DeprecationWarning): corr_mat.GEVP(0, sorted_list="Eigenvalue") - with pytest.warns(DeprecationWarning): - corr_mat.GEVP(0, state=0) - def test_GEVP_exceptions(): corr_aa = _gen_corr(1) corr_ab = 0.5 * corr_aa From f45b0aa6ecb9067323fde164630e02470b754e55 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Wed, 18 May 2022 10:54:09 +0100 Subject: [PATCH 13/13] tests: coverage of Corr class improved. --- tests/correlators_test.py | 20 +++++++++++++++++++- 1 file changed, 19 insertions(+), 1 deletion(-) diff --git a/tests/correlators_test.py b/tests/correlators_test.py index 90dac90a..97bbebb1 100644 --- a/tests/correlators_test.py +++ b/tests/correlators_test.py @@ -229,6 +229,7 @@ def test_matrix_corr(): corr_ab = 0.5 * corr_aa corr_mat = pe.Corr(np.array([[corr_aa, corr_ab], [corr_ab, corr_aa]])) + corr_mat.gamma_method() corr_mat.item(0, 0) for (ts, sort) in zip([None, 1, 1], ["Eigenvalue", "Eigenvector", None]): @@ -241,6 +242,8 @@ def test_matrix_corr(): assert np.all([o == 0 for o in corr_1 - corr_aa]) corr_mat.matrix_symmetric() + corr_mat.GEVP(0, state=0) + corr_mat.Eigenvalue(2, state=0) def test_GEVP_warnings(): @@ -263,9 +266,21 @@ def test_GEVP_exceptions(): corr_mat = pe.Corr(np.array([[corr_aa, corr_ab], [corr_ab, corr_aa]])) corr_mat.item(0, 0) + with pytest.raises(Exception): + corr_mat.item(0, 0).projected() + + with pytest.raises(Exception): + corr_mat.item(0, 0).GEVP(2) + + with pytest.raises(Exception): + corr_mat.item(0, 0).matrix_symmetric() + with pytest.raises(Exception): corr_mat.GEVP(0, 0, sort=None) + with pytest.raises(Exception): + corr_mat.GEVP(0, sort=None) + with pytest.raises(Exception): corr_mat.GEVP(1, 0, sort="Eigenvector") @@ -275,6 +290,9 @@ def test_GEVP_exceptions(): with pytest.raises(Exception): corr_mat.plottable() + with pytest.raises(Exception): + corr_mat.spaghetti_plot() + with pytest.raises(Exception): corr_mat.show() @@ -282,7 +300,7 @@ def test_GEVP_exceptions(): corr_mat.m_eff() with pytest.raises(Exception): - corr_mat.Hankel() + corr_mat.Hankel(2) with pytest.raises(Exception): corr_mat.plateau()