From 5da20365bf750892a54cf8b4ab9b0f2fdadfca46 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 22 Nov 2021 16:06:58 +0000 Subject: [PATCH] refactor!: npr module removed, Npr_matrix class moved to input.hadrons --- pyerrors/__init__.py | 1 - pyerrors/input/hadrons.py | 41 ++++++++++++++- pyerrors/npr.py | 108 -------------------------------------- 3 files changed, 40 insertions(+), 110 deletions(-) delete mode 100644 pyerrors/npr.py diff --git a/pyerrors/__init__.py b/pyerrors/__init__.py index c56af765..694f9017 100644 --- a/pyerrors/__init__.py +++ b/pyerrors/__init__.py @@ -251,7 +251,6 @@ from . import input from . import linalg from . import misc from . import mpm -from . import npr from . import roots from .version import __version__ diff --git a/pyerrors/input/hadrons.py b/pyerrors/input/hadrons.py index 7be3fd9d..b006d81b 100644 --- a/pyerrors/input/hadrons.py +++ b/pyerrors/input/hadrons.py @@ -6,7 +6,6 @@ import h5py import numpy as np from ..obs import Obs, CObs from ..correlators import Corr -from ..npr import Npr_matrix def _get_files(path, filestem): @@ -84,6 +83,46 @@ def read_meson_hd5(path, filestem, ens_id, meson='meson_0', tree='meson'): return corr +class Npr_matrix(np.ndarray): + + def __new__(cls, input_array, mom_in=None, mom_out=None): + obj = np.asarray(input_array).view(cls) + obj.mom_in = mom_in + obj.mom_out = mom_out + return obj + + @property + def g5H(self): + """Gamma_5 hermitean conjugate + + Uses the fact that the propagator is gamma5 hermitean, so just the + in and out momenta of the propagator are exchanged. + """ + return Npr_matrix(self, + mom_in=self.mom_out, + mom_out=self.mom_in) + + def _propagate_mom(self, other, name): + s_mom = getattr(self, name, None) + o_mom = getattr(other, name, None) + if s_mom is not None and o_mom is not None: + if not np.allclose(s_mom, o_mom): + raise Exception(name + ' does not match.') + return o_mom if o_mom is not None else s_mom + + def __matmul__(self, other): + return self.__new__(Npr_matrix, + super().__matmul__(other), + self._propagate_mom(other, 'mom_in'), + self._propagate_mom(other, 'mom_out')) + + def __array_finalize__(self, obj): + if obj is None: + return + self.mom_in = getattr(obj, 'mom_in', None) + self.mom_out = getattr(obj, 'mom_out', None) + + def read_ExternalLeg_hd5(path, filestem, ens_id, order='F'): """Read hadrons ExternalLeg hdf5 file and output an array of CObs diff --git a/pyerrors/npr.py b/pyerrors/npr.py deleted file mode 100644 index 5e67ce9e..00000000 --- a/pyerrors/npr.py +++ /dev/null @@ -1,108 +0,0 @@ -import warnings -import numpy as np -from .linalg import inv, matmul -from .dirac import gamma - - -L = None -T = None - - -class Npr_matrix(np.ndarray): - - def __new__(cls, input_array, mom_in=None, mom_out=None): - obj = np.asarray(input_array).view(cls) - obj.mom_in = mom_in - obj.mom_out = mom_out - return obj - - @property - def g5H(self): - """Gamma_5 hermitean conjugate - - Uses the fact that the propagator is gamma5 hermitean, so just the - in and out momenta of the propagator are exchanged. - """ - return Npr_matrix(self, - mom_in=self.mom_out, - mom_out=self.mom_in) - - def _propagate_mom(self, other, name): - s_mom = getattr(self, name, None) - o_mom = getattr(other, name, None) - if s_mom is not None and o_mom is not None: - if not np.allclose(s_mom, o_mom): - raise Exception(name + ' does not match.') - return o_mom if o_mom is not None else s_mom - - def __matmul__(self, other): - return self.__new__(Npr_matrix, - super().__matmul__(other), - self._propagate_mom(other, 'mom_in'), - self._propagate_mom(other, 'mom_out')) - - def __array_finalize__(self, obj): - if obj is None: - return - self.mom_in = getattr(obj, 'mom_in', None) - self.mom_out = getattr(obj, 'mom_out', None) - - -def _check_geometry(): - if L is None: - raise Exception("Spatial extent 'L' not set.") - else: - if not isinstance(L, int): - raise Exception("Spatial extent 'L' must be an integer.") - if T is None: - raise Exception("Temporal extent 'T' not set.") - if not isinstance(T, int): - raise Exception("Temporal extent 'T' must be an integer.") - - -def inv_propagator(prop): - """ Inverts a 12x12 quark propagator""" - if prop.shape != (12, 12): - raise Exception("Only 12x12 propagators can be inverted.") - return Npr_matrix(inv(prop), prop.mom_in) - - -def Zq(inv_prop, fermion='Wilson'): - """ Calculates the quark field renormalization constant Zq - - Parameters - ---------- - inv_prop : array - Inverted 12x12 quark propagator - fermion : str - Fermion type for which the tree-level propagator is used - in the calculation of Zq. Default Wilson. - """ - _check_geometry() - mom = np.copy(inv_prop.mom_in) - mom[3] /= T / L - sin_mom = np.sin(2 * np.pi / L * mom) - - if fermion == 'Wilson': - p_slash = -1j * (sin_mom[0] * gamma[0] + sin_mom[1] * gamma[1] + sin_mom[2] * gamma[2] + sin_mom[3] * gamma[3]) / np.sum(sin_mom ** 2) - elif fermion == 'Continuum': - p_mom = 2 * np.pi / L * mom - p_slash = -1j * (p_mom[0] * gamma[0] + p_mom[1] * gamma[1] + p_mom[2] * gamma[2] + p_mom[3] * gamma[3]) / np.sum(p_mom ** 2) - elif fermion == 'DWF': - W = np.sum(1 - np.cos(2 * np.pi / L * mom)) - s2 = np.sum(sin_mom ** 2) - p_slash = -1j * (sin_mom[0] * gamma[0] + sin_mom[1] * gamma[1] + sin_mom[2] * gamma[2] + sin_mom[3] * gamma[3]) - p_slash /= 2 * (W - 1 + np.sqrt((1 - W) ** 2 + s2)) - else: - raise Exception("Fermion type '" + fermion + "' not implemented") - - res = 1 / 12. * np.trace(matmul(inv_prop, np.kron(np.eye(3, dtype=int), p_slash))) - res.gamma_method() - - if not res.imag.is_zero_within_error(5): - warnings.warn("Imaginary part of Zq is not zero within 5 sigma") - return res - - res.real.tag = "Zq '" + fermion + "', p=" + str(inv_prop.mom_in) - - return res.real