pyerrors.input.hadrons

View Source
  0import os
  1import warnings
  2from collections import Counter
  3import h5py
  4import numpy as np
  5from ..obs import Obs, CObs
  6from ..correlators import Corr
  7from ..dirac import epsilon_tensor_rank4
  8
  9
 10def _get_files(path, filestem, idl):
 11    ls = os.listdir(path)
 12
 13    # Clean up file list
 14    files = list(filter(lambda x: x.startswith(filestem), ls))
 15
 16    if not files:
 17        raise Exception('No files starting with', filestem, 'in folder', path)
 18
 19    def get_cnfg_number(n):
 20        return int(n.replace(".h5", "")[len(filestem) + 1:])  # From python 3.9 onward the safer 'removesuffix' method can be used.
 21
 22    # Sort according to configuration number
 23    files.sort(key=get_cnfg_number)
 24
 25    cnfg_numbers = []
 26    filtered_files = []
 27    for line in files:
 28        no = get_cnfg_number(line)
 29        if idl:
 30            if no in list(idl):
 31                filtered_files.append(line)
 32                cnfg_numbers.append(no)
 33        else:
 34            filtered_files.append(line)
 35            cnfg_numbers.append(no)
 36
 37    if idl:
 38        if Counter(list(idl)) != Counter(cnfg_numbers):
 39            raise Exception("Not all configurations specified in idl found, configurations " + str(list(Counter(list(idl)) - Counter(cnfg_numbers))) + " are missing.")
 40
 41    # Check that configurations are evenly spaced
 42    dc = np.unique(np.diff(cnfg_numbers))
 43    if np.any(dc < 0):
 44        raise Exception("Unsorted files")
 45    if len(dc) == 1:
 46        idx = range(cnfg_numbers[0], cnfg_numbers[-1] + dc[0], dc[0])
 47    elif idl:
 48        idx = idl
 49        warnings.warn("Configurations are not evenly spaced.", RuntimeWarning)
 50    else:
 51        raise Exception("Configurations are not evenly spaced.")
 52
 53    return filtered_files, idx
 54
 55
 56def read_meson_hd5(path, filestem, ens_id, meson='meson_0', idl=None):
 57    """Read hadrons meson hdf5 file and extract the meson labeled 'meson'
 58
 59    Parameters
 60    -----------------
 61    path : str
 62        path to the files to read
 63    filestem : str
 64        namestem of the files to read
 65    ens_id : str
 66        name of the ensemble, required for internal bookkeeping
 67    meson : str
 68        label of the meson to be extracted, standard value meson_0 which
 69        corresponds to the pseudoscalar pseudoscalar two-point function.
 70    idl : range
 71        If specified only configurations in the given range are read in.
 72    """
 73
 74    files, idx = _get_files(path, filestem, idl)
 75
 76    tree = meson.rsplit('_')[0]
 77    corr_data = []
 78    infos = []
 79    for hd5_file in files:
 80        h5file = h5py.File(path + '/' + hd5_file, "r")
 81        if not tree + '/' + meson in h5file:
 82            raise Exception("Entry '" + meson + "' not contained in the files.")
 83        raw_data = list(h5file[tree + '/' + meson + '/corr'])
 84        real_data = [o[0] for o in raw_data]
 85        corr_data.append(real_data)
 86        if not infos:
 87            for k, i in h5file[tree + '/' + meson].attrs.items():
 88                infos.append(k + ': ' + i[0].decode())
 89        h5file.close()
 90    corr_data = np.array(corr_data)
 91
 92    l_obs = []
 93    for c in corr_data.T:
 94        l_obs.append(Obs([c], [ens_id], idl=[idx]))
 95
 96    corr = Corr(l_obs)
 97    corr.tag = r", ".join(infos)
 98    return corr
 99
100
101class Npr_matrix(np.ndarray):
102
103    def __new__(cls, input_array, mom_in=None, mom_out=None):
104        obj = np.asarray(input_array).view(cls)
105        obj.mom_in = mom_in
106        obj.mom_out = mom_out
107        return obj
108
109    @property
110    def g5H(self):
111        """Gamma_5 hermitean conjugate
112
113        Uses the fact that the propagator is gamma5 hermitean, so just the
114        in and out momenta of the propagator are exchanged.
115        """
116        return Npr_matrix(self,
117                          mom_in=self.mom_out,
118                          mom_out=self.mom_in)
119
120    def _propagate_mom(self, other, name):
121        s_mom = getattr(self, name, None)
122        o_mom = getattr(other, name, None)
123        if s_mom is not None and o_mom is not None:
124            if not np.allclose(s_mom, o_mom):
125                raise Exception(name + ' does not match.')
126        return o_mom if o_mom is not None else s_mom
127
128    def __matmul__(self, other):
129        return self.__new__(Npr_matrix,
130                            super().__matmul__(other),
131                            self._propagate_mom(other, 'mom_in'),
132                            self._propagate_mom(other, 'mom_out'))
133
134    def __array_finalize__(self, obj):
135        if obj is None:
136            return
137        self.mom_in = getattr(obj, 'mom_in', None)
138        self.mom_out = getattr(obj, 'mom_out', None)
139
140
141def read_ExternalLeg_hd5(path, filestem, ens_id, idl=None):
142    """Read hadrons ExternalLeg hdf5 file and output an array of CObs
143
144    Parameters
145    ----------
146    path : str
147        path to the files to read
148    filestem : str
149        namestem of the files to read
150    ens_id : str
151        name of the ensemble, required for internal bookkeeping
152    idl : range
153        If specified only configurations in the given range are read in.
154    """
155
156    files, idx = _get_files(path, filestem, idl)
157
158    mom = None
159
160    corr_data = []
161    for hd5_file in files:
162        file = h5py.File(path + '/' + hd5_file, "r")
163        raw_data = file['ExternalLeg/corr'][0][0].view('complex')
164        corr_data.append(raw_data)
165        if mom is None:
166            mom = np.array(str(file['ExternalLeg/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float)
167        file.close()
168    corr_data = np.array(corr_data)
169
170    rolled_array = np.rollaxis(corr_data, 0, 5)
171
172    matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
173    for si, sj, ci, cj in np.ndindex(rolled_array.shape[:-1]):
174        real = Obs([rolled_array[si, sj, ci, cj].real], [ens_id], idl=[idx])
175        imag = Obs([rolled_array[si, sj, ci, cj].imag], [ens_id], idl=[idx])
176        matrix[si, sj, ci, cj] = CObs(real, imag)
177
178    return Npr_matrix(matrix, mom_in=mom)
179
180
181def read_Bilinear_hd5(path, filestem, ens_id, idl=None):
182    """Read hadrons Bilinear hdf5 file and output an array of CObs
183
184    Parameters
185    ----------
186    path : str
187        path to the files to read
188    filestem : str
189        namestem of the files to read
190    ens_id : str
191        name of the ensemble, required for internal bookkeeping
192    idl : range
193        If specified only configurations in the given range are read in.
194    """
195
196    files, idx = _get_files(path, filestem, idl)
197
198    mom_in = None
199    mom_out = None
200
201    corr_data = {}
202    for hd5_file in files:
203        file = h5py.File(path + '/' + hd5_file, "r")
204        for i in range(16):
205            name = file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['gamma'][0].decode('UTF-8')
206            if name not in corr_data:
207                corr_data[name] = []
208            raw_data = file['Bilinear/Bilinear_' + str(i) + '/corr'][0][0].view('complex')
209            corr_data[name].append(raw_data)
210            if mom_in is None:
211                mom_in = np.array(str(file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float)
212            if mom_out is None:
213                mom_out = np.array(str(file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['pOut'])[3:-2].strip().split(), dtype=float)
214
215        file.close()
216
217    result_dict = {}
218
219    for key, data in corr_data.items():
220        local_data = np.array(data)
221
222        rolled_array = np.rollaxis(local_data, 0, 5)
223
224        matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
225        for si, sj, ci, cj in np.ndindex(rolled_array.shape[:-1]):
226            real = Obs([rolled_array[si, sj, ci, cj].real], [ens_id], idl=[idx])
227            imag = Obs([rolled_array[si, sj, ci, cj].imag], [ens_id], idl=[idx])
228            matrix[si, sj, ci, cj] = CObs(real, imag)
229
230        result_dict[key] = Npr_matrix(matrix, mom_in=mom_in, mom_out=mom_out)
231
232    return result_dict
233
234
235def read_Fourquark_hd5(path, filestem, ens_id, idl=None, vertices=["VA", "AV"]):
236    """Read hadrons FourquarkFullyConnected hdf5 file and output an array of CObs
237
238    Parameters
239    ----------
240    path : str
241        path to the files to read
242    filestem : str
243        namestem of the files to read
244    ens_id : str
245        name of the ensemble, required for internal bookkeeping
246    idl : range
247        If specified only configurations in the given range are read in.
248    vertices : list
249        Vertex functions to be extracted.
250    """
251
252    files, idx = _get_files(path, filestem, idl)
253
254    mom_in = None
255    mom_out = None
256
257    vertex_names = []
258    for vertex in vertices:
259        vertex_names += _get_lorentz_names(vertex)
260
261    corr_data = {}
262
263    tree = 'FourQuarkFullyConnected/FourQuarkFullyConnected_'
264
265    for hd5_file in files:
266        file = h5py.File(path + '/' + hd5_file, "r")
267
268        for i in range(32):
269            name = (file[tree + str(i) + '/info'].attrs['gammaA'][0].decode('UTF-8'), file[tree + str(i) + '/info'].attrs['gammaB'][0].decode('UTF-8'))
270            if name in vertex_names:
271                if name not in corr_data:
272                    corr_data[name] = []
273                raw_data = file[tree + str(i) + '/corr'][0][0].view('complex')
274                corr_data[name].append(raw_data)
275                if mom_in is None:
276                    mom_in = np.array(str(file[tree + str(i) + '/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float)
277                if mom_out is None:
278                    mom_out = np.array(str(file[tree + str(i) + '/info'].attrs['pOut'])[3:-2].strip().split(), dtype=float)
279
280        file.close()
281
282    intermediate_dict = {}
283
284    for vertex in vertices:
285        lorentz_names = _get_lorentz_names(vertex)
286        for v_name in lorentz_names:
287            if v_name in [('SigmaXY', 'SigmaZT'),
288                          ('SigmaXT', 'SigmaYZ'),
289                          ('SigmaYZ', 'SigmaXT'),
290                          ('SigmaZT', 'SigmaXY')]:
291                sign = -1
292            else:
293                sign = 1
294            if vertex not in intermediate_dict:
295                intermediate_dict[vertex] = sign * np.array(corr_data[v_name])
296            else:
297                intermediate_dict[vertex] += sign * np.array(corr_data[v_name])
298
299    result_dict = {}
300
301    for key, data in intermediate_dict.items():
302
303        rolled_array = np.moveaxis(data, 0, 8)
304
305        matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
306        for index in np.ndindex(rolled_array.shape[:-1]):
307            real = Obs([rolled_array[index].real], [ens_id], idl=[idx])
308            imag = Obs([rolled_array[index].imag], [ens_id], idl=[idx])
309            matrix[index] = CObs(real, imag)
310
311        result_dict[key] = Npr_matrix(matrix, mom_in=mom_in, mom_out=mom_out)
312
313    return result_dict
314
315
316def _get_lorentz_names(name):
317    lorentz_index = ['X', 'Y', 'Z', 'T']
318
319    res = []
320
321    if name == "TT":
322        for i in range(4):
323            for j in range(i + 1, 4):
324                res.append(("Sigma" + lorentz_index[i] + lorentz_index[j], "Sigma" + lorentz_index[i] + lorentz_index[j]))
325        return res
326
327    if name == "TTtilde":
328        for i in range(4):
329            for j in range(i + 1, 4):
330                for k in range(4):
331                    for o in range(k + 1, 4):
332                        fac = epsilon_tensor_rank4(i, j, k, o)
333                        if not np.isclose(fac, 0.0):
334                            res.append(("Sigma" + lorentz_index[i] + lorentz_index[j], "Sigma" + lorentz_index[k] + lorentz_index[o]))
335        return res
336
337    assert len(name) == 2
338
339    if 'S' in name or 'P' in name:
340        if not set(name) <= set(['S', 'P']):
341            raise Exception("'" + name + "' is not a Lorentz scalar")
342
343        g_names = {'S': 'Identity',
344                   'P': 'Gamma5'}
345
346        res.append((g_names[name[0]], g_names[name[1]]))
347
348    else:
349        if not set(name) <= set(['V', 'A']):
350            raise Exception("'" + name + "' is not a Lorentz scalar")
351
352        for ind in lorentz_index:
353            res.append(('Gamma' + ind + (name[0] == 'A') * 'Gamma5',
354                        'Gamma' + ind + (name[1] == 'A') * 'Gamma5'))
355
356    return res
#   def read_meson_hd5(path, filestem, ens_id, meson='meson_0', idl=None):
View Source
57def read_meson_hd5(path, filestem, ens_id, meson='meson_0', idl=None):
58    """Read hadrons meson hdf5 file and extract the meson labeled 'meson'
59
60    Parameters
61    -----------------
62    path : str
63        path to the files to read
64    filestem : str
65        namestem of the files to read
66    ens_id : str
67        name of the ensemble, required for internal bookkeeping
68    meson : str
69        label of the meson to be extracted, standard value meson_0 which
70        corresponds to the pseudoscalar pseudoscalar two-point function.
71    idl : range
72        If specified only configurations in the given range are read in.
73    """
74
75    files, idx = _get_files(path, filestem, idl)
76
77    tree = meson.rsplit('_')[0]
78    corr_data = []
79    infos = []
80    for hd5_file in files:
81        h5file = h5py.File(path + '/' + hd5_file, "r")
82        if not tree + '/' + meson in h5file:
83            raise Exception("Entry '" + meson + "' not contained in the files.")
84        raw_data = list(h5file[tree + '/' + meson + '/corr'])
85        real_data = [o[0] for o in raw_data]
86        corr_data.append(real_data)
87        if not infos:
88            for k, i in h5file[tree + '/' + meson].attrs.items():
89                infos.append(k + ': ' + i[0].decode())
90        h5file.close()
91    corr_data = np.array(corr_data)
92
93    l_obs = []
94    for c in corr_data.T:
95        l_obs.append(Obs([c], [ens_id], idl=[idx]))
96
97    corr = Corr(l_obs)
98    corr.tag = r", ".join(infos)
99    return corr

Read hadrons meson hdf5 file and extract the meson labeled 'meson'

Parameters
  • path (str): path to the files to read
  • filestem (str): namestem of the files to read
  • ens_id (str): name of the ensemble, required for internal bookkeeping
  • meson (str): label of the meson to be extracted, standard value meson_0 which corresponds to the pseudoscalar pseudoscalar two-point function.
  • idl (range): If specified only configurations in the given range are read in.
#   class Npr_matrix(numpy.ndarray):
View Source
102class Npr_matrix(np.ndarray):
103
104    def __new__(cls, input_array, mom_in=None, mom_out=None):
105        obj = np.asarray(input_array).view(cls)
106        obj.mom_in = mom_in
107        obj.mom_out = mom_out
108        return obj
109
110    @property
111    def g5H(self):
112        """Gamma_5 hermitean conjugate
113
114        Uses the fact that the propagator is gamma5 hermitean, so just the
115        in and out momenta of the propagator are exchanged.
116        """
117        return Npr_matrix(self,
118                          mom_in=self.mom_out,
119                          mom_out=self.mom_in)
120
121    def _propagate_mom(self, other, name):
122        s_mom = getattr(self, name, None)
123        o_mom = getattr(other, name, None)
124        if s_mom is not None and o_mom is not None:
125            if not np.allclose(s_mom, o_mom):
126                raise Exception(name + ' does not match.')
127        return o_mom if o_mom is not None else s_mom
128
129    def __matmul__(self, other):
130        return self.__new__(Npr_matrix,
131                            super().__matmul__(other),
132                            self._propagate_mom(other, 'mom_in'),
133                            self._propagate_mom(other, 'mom_out'))
134
135    def __array_finalize__(self, obj):
136        if obj is None:
137            return
138        self.mom_in = getattr(obj, 'mom_in', None)
139        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: An ndarray alias :term:generic <generic type> w.r.t. its dtype.type <numpy.dtype.type>.

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

Uses the fact that the propagator is gamma5 hermitean, so just the in and out momenta of the propagator are exchanged.

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 read_ExternalLeg_hd5(path, filestem, ens_id, idl=None):
View Source
142def read_ExternalLeg_hd5(path, filestem, ens_id, idl=None):
143    """Read hadrons ExternalLeg hdf5 file and output an array of CObs
144
145    Parameters
146    ----------
147    path : str
148        path to the files to read
149    filestem : str
150        namestem of the files to read
151    ens_id : str
152        name of the ensemble, required for internal bookkeeping
153    idl : range
154        If specified only configurations in the given range are read in.
155    """
156
157    files, idx = _get_files(path, filestem, idl)
158
159    mom = None
160
161    corr_data = []
162    for hd5_file in files:
163        file = h5py.File(path + '/' + hd5_file, "r")
164        raw_data = file['ExternalLeg/corr'][0][0].view('complex')
165        corr_data.append(raw_data)
166        if mom is None:
167            mom = np.array(str(file['ExternalLeg/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float)
168        file.close()
169    corr_data = np.array(corr_data)
170
171    rolled_array = np.rollaxis(corr_data, 0, 5)
172
173    matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
174    for si, sj, ci, cj in np.ndindex(rolled_array.shape[:-1]):
175        real = Obs([rolled_array[si, sj, ci, cj].real], [ens_id], idl=[idx])
176        imag = Obs([rolled_array[si, sj, ci, cj].imag], [ens_id], idl=[idx])
177        matrix[si, sj, ci, cj] = CObs(real, imag)
178
179    return Npr_matrix(matrix, mom_in=mom)

Read hadrons ExternalLeg hdf5 file and output an array of CObs

Parameters
  • path (str): path to the files to read
  • filestem (str): namestem of the files to read
  • ens_id (str): name of the ensemble, required for internal bookkeeping
  • idl (range): If specified only configurations in the given range are read in.
#   def read_Bilinear_hd5(path, filestem, ens_id, idl=None):
View Source
182def read_Bilinear_hd5(path, filestem, ens_id, idl=None):
183    """Read hadrons Bilinear hdf5 file and output an array of CObs
184
185    Parameters
186    ----------
187    path : str
188        path to the files to read
189    filestem : str
190        namestem of the files to read
191    ens_id : str
192        name of the ensemble, required for internal bookkeeping
193    idl : range
194        If specified only configurations in the given range are read in.
195    """
196
197    files, idx = _get_files(path, filestem, idl)
198
199    mom_in = None
200    mom_out = None
201
202    corr_data = {}
203    for hd5_file in files:
204        file = h5py.File(path + '/' + hd5_file, "r")
205        for i in range(16):
206            name = file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['gamma'][0].decode('UTF-8')
207            if name not in corr_data:
208                corr_data[name] = []
209            raw_data = file['Bilinear/Bilinear_' + str(i) + '/corr'][0][0].view('complex')
210            corr_data[name].append(raw_data)
211            if mom_in is None:
212                mom_in = np.array(str(file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float)
213            if mom_out is None:
214                mom_out = np.array(str(file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['pOut'])[3:-2].strip().split(), dtype=float)
215
216        file.close()
217
218    result_dict = {}
219
220    for key, data in corr_data.items():
221        local_data = np.array(data)
222
223        rolled_array = np.rollaxis(local_data, 0, 5)
224
225        matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
226        for si, sj, ci, cj in np.ndindex(rolled_array.shape[:-1]):
227            real = Obs([rolled_array[si, sj, ci, cj].real], [ens_id], idl=[idx])
228            imag = Obs([rolled_array[si, sj, ci, cj].imag], [ens_id], idl=[idx])
229            matrix[si, sj, ci, cj] = CObs(real, imag)
230
231        result_dict[key] = Npr_matrix(matrix, mom_in=mom_in, mom_out=mom_out)
232
233    return result_dict

Read hadrons Bilinear hdf5 file and output an array of CObs

Parameters
  • path (str): path to the files to read
  • filestem (str): namestem of the files to read
  • ens_id (str): name of the ensemble, required for internal bookkeeping
  • idl (range): If specified only configurations in the given range are read in.
#   def read_Fourquark_hd5(path, filestem, ens_id, idl=None, vertices=['VA', 'AV']):
View Source
236def read_Fourquark_hd5(path, filestem, ens_id, idl=None, vertices=["VA", "AV"]):
237    """Read hadrons FourquarkFullyConnected hdf5 file and output an array of CObs
238
239    Parameters
240    ----------
241    path : str
242        path to the files to read
243    filestem : str
244        namestem of the files to read
245    ens_id : str
246        name of the ensemble, required for internal bookkeeping
247    idl : range
248        If specified only configurations in the given range are read in.
249    vertices : list
250        Vertex functions to be extracted.
251    """
252
253    files, idx = _get_files(path, filestem, idl)
254
255    mom_in = None
256    mom_out = None
257
258    vertex_names = []
259    for vertex in vertices:
260        vertex_names += _get_lorentz_names(vertex)
261
262    corr_data = {}
263
264    tree = 'FourQuarkFullyConnected/FourQuarkFullyConnected_'
265
266    for hd5_file in files:
267        file = h5py.File(path + '/' + hd5_file, "r")
268
269        for i in range(32):
270            name = (file[tree + str(i) + '/info'].attrs['gammaA'][0].decode('UTF-8'), file[tree + str(i) + '/info'].attrs['gammaB'][0].decode('UTF-8'))
271            if name in vertex_names:
272                if name not in corr_data:
273                    corr_data[name] = []
274                raw_data = file[tree + str(i) + '/corr'][0][0].view('complex')
275                corr_data[name].append(raw_data)
276                if mom_in is None:
277                    mom_in = np.array(str(file[tree + str(i) + '/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float)
278                if mom_out is None:
279                    mom_out = np.array(str(file[tree + str(i) + '/info'].attrs['pOut'])[3:-2].strip().split(), dtype=float)
280
281        file.close()
282
283    intermediate_dict = {}
284
285    for vertex in vertices:
286        lorentz_names = _get_lorentz_names(vertex)
287        for v_name in lorentz_names:
288            if v_name in [('SigmaXY', 'SigmaZT'),
289                          ('SigmaXT', 'SigmaYZ'),
290                          ('SigmaYZ', 'SigmaXT'),
291                          ('SigmaZT', 'SigmaXY')]:
292                sign = -1
293            else:
294                sign = 1
295            if vertex not in intermediate_dict:
296                intermediate_dict[vertex] = sign * np.array(corr_data[v_name])
297            else:
298                intermediate_dict[vertex] += sign * np.array(corr_data[v_name])
299
300    result_dict = {}
301
302    for key, data in intermediate_dict.items():
303
304        rolled_array = np.moveaxis(data, 0, 8)
305
306        matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
307        for index in np.ndindex(rolled_array.shape[:-1]):
308            real = Obs([rolled_array[index].real], [ens_id], idl=[idx])
309            imag = Obs([rolled_array[index].imag], [ens_id], idl=[idx])
310            matrix[index] = CObs(real, imag)
311
312        result_dict[key] = Npr_matrix(matrix, mom_in=mom_in, mom_out=mom_out)
313
314    return result_dict

Read hadrons FourquarkFullyConnected hdf5 file and output an array of CObs

Parameters
  • path (str): path to the files to read
  • filestem (str): namestem of the files to read
  • ens_id (str): name of the ensemble, required for internal bookkeeping
  • idl (range): If specified only configurations in the given range are read in.
  • vertices (list): Vertex functions to be extracted.