pyerrors.npr

View Source 
import warnings
import numpy as np
from .linalg import inv, matmul
from .dirac import gamma, gamma5


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

        Returns gamma_5 @ M.T.conj() @ gamma_5 and exchanges in and out going
        momenta. Works only for 12x12 matrices.
        """
        if self.shape != (12, 12):
            raise Exception('g5H only works for 12x12 matrices.')
        extended_g5 = np.kron(np.eye(3, dtype=int), gamma5)
        return Npr_matrix(matmul(extended_g5, self.conj().T, extended_g5),
                          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
    return res.real
#   class Npr_matrix(numpy.ndarray):
View Source 
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

        Returns gamma_5 @ M.T.conj() @ gamma_5 and exchanges in and out going
        momenta. Works only for 12x12 matrices.
        """
        if self.shape != (12, 12):
            raise Exception('g5H only works for 12x12 matrices.')
        extended_g5 = np.kron(np.eye(3, dtype=int), gamma5)
        return Npr_matrix(matmul(extended_g5, self.conj().T, extended_g5),
                          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)

ndarray(shape, dtype=float, buffer=None, offset=0, strides=None, order=None)

An array object represents a multidimensional, homogeneous array of fixed-size items. An associated data-type object describes the format of each element in the array (its byte-order, how many bytes it occupies in memory, whether it is an integer, a floating point number, or something else, etc.)

Arrays should be constructed using array, zeros or empty (refer to the See Also section below). The parameters given here refer to a low-level method (ndarray(...)) for instantiating an array.

For more information, refer to the numpy module and examine the methods and attributes of an array.

Parameters
  • (for the __new__ method; see Notes below)
  • shape (tuple of ints): Shape of created array.
  • dtype (data-type, optional): Any object that can be interpreted as a numpy data type.
  • buffer (object exposing buffer interface, optional): Used to fill the array with data.
  • offset (int, optional): Offset of array data in buffer.
  • strides (tuple of ints, optional): Strides of data in memory.
  • order ({'C', 'F'}, optional): Row-major (C-style) or column-major (Fortran-style) order.
Attributes
  • T (ndarray): Transpose of the array.
  • data (buffer): The array's elements, in memory.
  • dtype (dtype object): Describes the format of the elements in the array.
  • flags (dict): Dictionary containing information related to memory use, e.g., 'C_CONTIGUOUS', 'OWNDATA', 'WRITEABLE', etc.
  • flat (numpy.flatiter object): Flattened version of the array as an iterator. The iterator allows assignments, e.g., x.flat = 3 (See ndarray.flat for assignment examples; TODO).
  • imag (ndarray): Imaginary part of the array.
  • real (ndarray): Real part of the array.
  • size (int): Number of elements in the array.
  • itemsize (int): The memory use of each array element in bytes.
  • nbytes (int): The total number of bytes required to store the array data, i.e., itemsize * size.
  • ndim (int): The array's number of dimensions.
  • shape (tuple of ints): Shape of the array.
  • strides (tuple of ints): The step-size required to move from one element to the next in memory. For example, a contiguous (3, 4) array of type int16 in C-order has strides (8, 2). This implies that to move from element to element in memory requires jumps of 2 bytes. To move from row-to-row, one needs to jump 8 bytes at a time (2 * 4).
  • ctypes (ctypes object): Class containing properties of the array needed for interaction with ctypes.
  • base (ndarray): If the array is a view into another array, that array is its base (unless that array is also a view). The base array is where the array data is actually stored.
See Also

array: Construct an array.
zeros: Create an array, each element of which is zero.
empty: Create an array, but leave its allocated memory unchanged (i.e., it contains "garbage").
dtype: Create a data-type.
numpy.typing.NDArray: A :term:generic <generic type> version of ndarray.

Notes

There are two modes of creating an array using __new__:

  1. If buffer is None, then only shape, dtype, and order are used.
  2. If buffer is an object exposing the buffer interface, then all keywords are interpreted.

No __init__ method is needed because the array is fully initialized after the __new__ method.

Examples

These examples illustrate the low-level ndarray constructor. Refer to the See Also section above for easier ways of constructing an ndarray.

First mode, buffer is None:

>>> np.ndarray(shape=(2,2), dtype=float, order='F')
array([[0.0e+000, 0.0e+000], # random
       [     nan, 2.5e-323]])

Second mode:

>>> np.ndarray((2,), buffer=np.array([1,2,3]),
...            offset=np.int_().itemsize,
...            dtype=int) # offset = 1*itemsize, i.e. skip first element
array([2, 3])
#   Npr_matrix()
#   g5H

Gamma_5 hermitean conjugate

Returns gamma_5 @ M.T.conj() @ gamma_5 and exchanges in and out going momenta. Works only for 12x12 matrices.

Inherited Members
numpy.ndarray
dumps
dump
all
any
argmax
argmin
argpartition
argsort
astype
byteswap
choose
clip
compress
conj
conjugate
copy
cumprod
cumsum
diagonal
dot
fill
flatten
getfield
item
itemset
max
mean
min
newbyteorder
nonzero
partition
prod
ptp
put
ravel
repeat
reshape
resize
round
searchsorted
setfield
setflags
sort
squeeze
std
sum
swapaxes
take
tobytes
tofile
tolist
tostring
trace
transpose
var
view
ndim
flags
shape
strides
data
itemsize
size
nbytes
base
dtype
real
imag
flat
ctypes
T
#   def inv_propagator(prop):
View Source 
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)

Inverts a 12x12 quark propagator

#   def Zq(inv_prop, fermion='Wilson'):
View Source 
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
    return res.real

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.