pyerrors.input.hadrons

  1import os
  2import warnings
  3from collections import Counter
  4import h5py
  5from pathlib import Path
  6import numpy as np
  7from ..obs import Obs, CObs
  8from ..correlators import Corr
  9from ..dirac import epsilon_tensor_rank4
 10
 11
 12def _get_files(path, filestem, idl):
 13    ls = os.listdir(path)
 14
 15    # Clean up file list
 16    files = list(filter(lambda x: x.startswith(filestem + "."), ls))
 17
 18    if not files:
 19        raise Exception('No files starting with', filestem, 'in folder', path)
 20
 21    def get_cnfg_number(n):
 22        return int(n.replace(".h5", "")[len(filestem) + 1:])  # From python 3.9 onward the safer 'removesuffix' method can be used.
 23
 24    # Sort according to configuration number
 25    files.sort(key=get_cnfg_number)
 26
 27    cnfg_numbers = []
 28    filtered_files = []
 29    for line in files:
 30        no = get_cnfg_number(line)
 31        if idl:
 32            if no in list(idl):
 33                filtered_files.append(line)
 34                cnfg_numbers.append(no)
 35        else:
 36            filtered_files.append(line)
 37            cnfg_numbers.append(no)
 38
 39    if idl:
 40        if Counter(list(idl)) != Counter(cnfg_numbers):
 41            raise Exception("Not all configurations specified in idl found, configurations " + str(list(Counter(list(idl)) - Counter(cnfg_numbers))) + " are missing.")
 42
 43    # Check that configurations are evenly spaced
 44    dc = np.unique(np.diff(cnfg_numbers))
 45    if np.any(dc < 0):
 46        raise Exception("Unsorted files")
 47    if len(dc) == 1:
 48        idx = range(cnfg_numbers[0], cnfg_numbers[-1] + dc[0], dc[0])
 49    elif idl:
 50        idx = idl
 51        warnings.warn("Configurations are not evenly spaced.", RuntimeWarning)
 52    else:
 53        raise Exception("Configurations are not evenly spaced. Provide an idl if you want to proceed with this set of configurations.")
 54
 55    return filtered_files, idx
 56
 57
 58def read_meson_hd5(path, filestem, ens_id, meson='meson_0', idl=None, gammas=None):
 59    r'''Read hadrons meson hdf5 file and extract the meson labeled 'meson'
 60
 61    Parameters
 62    -----------------
 63    path : str
 64        path to the files to read
 65    filestem : str
 66        namestem of the files to read
 67    ens_id : str
 68        name of the ensemble, required for internal bookkeeping
 69    meson : str
 70        label of the meson to be extracted, standard value meson_0 which
 71        corresponds to the pseudoscalar pseudoscalar two-point function.
 72    gammas : tuple of strings
 73        Instrad of a meson label one can also provide a tuple of two strings
 74        indicating the gamma matrices at source and sink.
 75        ("Gamma5", "Gamma5") corresponds to the pseudoscalar pseudoscalar
 76        two-point function. The gammas argument dominateds over meson.
 77    idl : range
 78        If specified only configurations in the given range are read in.
 79    '''
 80
 81    files, idx = _get_files(path, filestem, idl)
 82
 83    tree = meson.rsplit('_')[0]
 84    if gammas is not None:
 85        h5file = h5py.File(path + '/' + files[0], "r")
 86        found_meson = None
 87        for key in h5file[tree].keys():
 88            if gammas[0] == h5file[tree][key].attrs["gamma_snk"][0].decode() and h5file[tree][key].attrs["gamma_src"][0].decode() == gammas[1]:
 89                found_meson = key
 90                break
 91        h5file.close()
 92        if found_meson:
 93            meson = found_meson
 94        else:
 95            raise Exception("Source Sink combination " + str(gammas) + " not found.")
 96
 97    corr_data = []
 98    infos = []
 99    for hd5_file in files:
100        h5file = h5py.File(path + '/' + hd5_file, "r")
101        if not tree + '/' + meson in h5file:
102            raise Exception("Entry '" + meson + "' not contained in the files.")
103        raw_data = h5file[tree + '/' + meson + '/corr']
104        real_data = raw_data[:]["re"].astype(np.double)
105        corr_data.append(real_data)
106        if not infos:
107            for k, i in h5file[tree + '/' + meson].attrs.items():
108                infos.append(k + ': ' + i[0].decode())
109        h5file.close()
110    corr_data = np.array(corr_data)
111
112    l_obs = []
113    for c in corr_data.T:
114        l_obs.append(Obs([c], [ens_id], idl=[idx]))
115
116    corr = Corr(l_obs)
117    corr.tag = r", ".join(infos)
118    return corr
119
120
121def read_DistillationContraction_hd5(path, ens_id, diagrams=["direct"], idl=None):
122    """Read hadrons DistillationContraction hdf5 files in given directory structure
123
124    Parameters
125    -----------------
126    path : str
127        path to the directories to read
128    ens_id : str
129        name of the ensemble, required for internal bookkeeping
130    diagrams : list
131        List of strings of the diagrams to extract, e.g. ["direct", "box", "cross"].
132    idl : range
133        If specified only configurations in the given range are read in.
134    """
135
136    res_dict = {}
137
138    directories, idx = _get_files(path, "data", idl)
139
140    explore_path = Path(path + "/" + directories[0])
141
142    for explore_file in explore_path.iterdir():
143        if explore_file.is_file():
144            stem = explore_file.with_suffix("").with_suffix("").as_posix().split("/")[-1]
145        else:
146            continue
147
148        file_list = []
149        for dir in directories:
150            tmp_path = Path(path + "/" + dir)
151            file_list.append((tmp_path / stem).as_posix() + tmp_path.suffix + ".h5")
152
153        corr_data = {}
154
155        for diagram in diagrams:
156            corr_data[diagram] = []
157
158        for n_file, (hd5_file, n_traj) in enumerate(zip(file_list, list(idx))):
159            h5file = h5py.File(hd5_file)
160
161            if n_file == 0:
162                if h5file["DistillationContraction/Metadata"].attrs.get("TimeSources")[0].decode() != "0...":
163                    raise Exception("Routine is only implemented for files containing inversions on all timeslices.")
164
165                Nt = h5file["DistillationContraction/Metadata"].attrs.get("Nt")[0]
166
167                identifier = []
168                for in_file in range(len(h5file["DistillationContraction/Metadata/DmfInputFiles"].attrs.keys()) - 1):
169                    encoded_info = h5file["DistillationContraction/Metadata/DmfInputFiles"].attrs.get("DmfInputFiles_" + str(in_file))
170                    full_info = encoded_info[0].decode().split("/")[-1].replace(".h5", "").split("_")
171                    my_tuple = (full_info[0], full_info[1][1:], full_info[2], full_info[3])
172                    identifier.append(my_tuple)
173                identifier = tuple(identifier)
174                # "DistillationContraction/Metadata/DmfSuffix" contains info about different quarks, irrelevant in the SU(3) case.
175
176            for diagram in diagrams:
177                real_data = np.zeros(Nt)
178                for x0 in range(Nt):
179                    raw_data = h5file["DistillationContraction/Correlators/" + diagram + "/" + str(x0)][:]["re"].astype(np.double)
180                    real_data += np.roll(raw_data, -x0)
181                real_data /= Nt
182
183                corr_data[diagram].append(real_data)
184            h5file.close()
185
186        res_dict[str(identifier)] = {}
187
188        for diagram in diagrams:
189
190            tmp_data = np.array(corr_data[diagram])
191
192            l_obs = []
193            for c in tmp_data.T:
194                l_obs.append(Obs([c], [ens_id], idl=[idx]))
195
196            corr = Corr(l_obs)
197            corr.tag = str(identifier)
198
199            res_dict[str(identifier)][diagram] = corr
200
201    return res_dict
202
203
204class Npr_matrix(np.ndarray):
205
206    def __new__(cls, input_array, mom_in=None, mom_out=None):
207        obj = np.asarray(input_array).view(cls)
208        obj.mom_in = mom_in
209        obj.mom_out = mom_out
210        return obj
211
212    @property
213    def g5H(self):
214        """Gamma_5 hermitean conjugate
215
216        Uses the fact that the propagator is gamma5 hermitean, so just the
217        in and out momenta of the propagator are exchanged.
218        """
219        return Npr_matrix(self,
220                          mom_in=self.mom_out,
221                          mom_out=self.mom_in)
222
223    def _propagate_mom(self, other, name):
224        s_mom = getattr(self, name, None)
225        o_mom = getattr(other, name, None)
226        if s_mom is not None and o_mom is not None:
227            if not np.allclose(s_mom, o_mom):
228                raise Exception(name + ' does not match.')
229        return o_mom if o_mom is not None else s_mom
230
231    def __matmul__(self, other):
232        return self.__new__(Npr_matrix,
233                            super().__matmul__(other),
234                            self._propagate_mom(other, 'mom_in'),
235                            self._propagate_mom(other, 'mom_out'))
236
237    def __array_finalize__(self, obj):
238        if obj is None:
239            return
240        self.mom_in = getattr(obj, 'mom_in', None)
241        self.mom_out = getattr(obj, 'mom_out', None)
242
243
244def read_ExternalLeg_hd5(path, filestem, ens_id, idl=None):
245    """Read hadrons ExternalLeg hdf5 file and output an array of CObs
246
247    Parameters
248    ----------
249    path : str
250        path to the files to read
251    filestem : str
252        namestem of the files to read
253    ens_id : str
254        name of the ensemble, required for internal bookkeeping
255    idl : range
256        If specified only configurations in the given range are read in.
257    """
258
259    files, idx = _get_files(path, filestem, idl)
260
261    mom = None
262
263    corr_data = []
264    for hd5_file in files:
265        file = h5py.File(path + '/' + hd5_file, "r")
266        raw_data = file['ExternalLeg/corr'][0][0].view('complex')
267        corr_data.append(raw_data)
268        if mom is None:
269            mom = np.array(str(file['ExternalLeg/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float)
270        file.close()
271    corr_data = np.array(corr_data)
272
273    rolled_array = np.rollaxis(corr_data, 0, 5)
274
275    matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
276    for si, sj, ci, cj in np.ndindex(rolled_array.shape[:-1]):
277        real = Obs([rolled_array[si, sj, ci, cj].real], [ens_id], idl=[idx])
278        imag = Obs([rolled_array[si, sj, ci, cj].imag], [ens_id], idl=[idx])
279        matrix[si, sj, ci, cj] = CObs(real, imag)
280
281    return Npr_matrix(matrix, mom_in=mom)
282
283
284def read_Bilinear_hd5(path, filestem, ens_id, idl=None):
285    """Read hadrons Bilinear hdf5 file and output an array of CObs
286
287    Parameters
288    ----------
289    path : str
290        path to the files to read
291    filestem : str
292        namestem of the files to read
293    ens_id : str
294        name of the ensemble, required for internal bookkeeping
295    idl : range
296        If specified only configurations in the given range are read in.
297    """
298
299    files, idx = _get_files(path, filestem, idl)
300
301    mom_in = None
302    mom_out = None
303
304    corr_data = {}
305    for hd5_file in files:
306        file = h5py.File(path + '/' + hd5_file, "r")
307        for i in range(16):
308            name = file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['gamma'][0].decode('UTF-8')
309            if name not in corr_data:
310                corr_data[name] = []
311            raw_data = file['Bilinear/Bilinear_' + str(i) + '/corr'][0][0].view('complex')
312            corr_data[name].append(raw_data)
313            if mom_in is None:
314                mom_in = np.array(str(file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float)
315            if mom_out is None:
316                mom_out = np.array(str(file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['pOut'])[3:-2].strip().split(), dtype=float)
317
318        file.close()
319
320    result_dict = {}
321
322    for key, data in corr_data.items():
323        local_data = np.array(data)
324
325        rolled_array = np.rollaxis(local_data, 0, 5)
326
327        matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
328        for si, sj, ci, cj in np.ndindex(rolled_array.shape[:-1]):
329            real = Obs([rolled_array[si, sj, ci, cj].real], [ens_id], idl=[idx])
330            imag = Obs([rolled_array[si, sj, ci, cj].imag], [ens_id], idl=[idx])
331            matrix[si, sj, ci, cj] = CObs(real, imag)
332
333        result_dict[key] = Npr_matrix(matrix, mom_in=mom_in, mom_out=mom_out)
334
335    return result_dict
336
337
338def read_Fourquark_hd5(path, filestem, ens_id, idl=None, vertices=["VA", "AV"]):
339    """Read hadrons FourquarkFullyConnected hdf5 file and output an array of CObs
340
341    Parameters
342    ----------
343    path : str
344        path to the files to read
345    filestem : str
346        namestem of the files to read
347    ens_id : str
348        name of the ensemble, required for internal bookkeeping
349    idl : range
350        If specified only configurations in the given range are read in.
351    vertices : list
352        Vertex functions to be extracted.
353    """
354
355    files, idx = _get_files(path, filestem, idl)
356
357    mom_in = None
358    mom_out = None
359
360    vertex_names = []
361    for vertex in vertices:
362        vertex_names += _get_lorentz_names(vertex)
363
364    corr_data = {}
365
366    tree = 'FourQuarkFullyConnected/FourQuarkFullyConnected_'
367
368    for hd5_file in files:
369        file = h5py.File(path + '/' + hd5_file, "r")
370
371        for i in range(32):
372            name = (file[tree + str(i) + '/info'].attrs['gammaA'][0].decode('UTF-8'), file[tree + str(i) + '/info'].attrs['gammaB'][0].decode('UTF-8'))
373            if name in vertex_names:
374                if name not in corr_data:
375                    corr_data[name] = []
376                raw_data = file[tree + str(i) + '/corr'][0][0].view('complex')
377                corr_data[name].append(raw_data)
378                if mom_in is None:
379                    mom_in = np.array(str(file[tree + str(i) + '/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float)
380                if mom_out is None:
381                    mom_out = np.array(str(file[tree + str(i) + '/info'].attrs['pOut'])[3:-2].strip().split(), dtype=float)
382
383        file.close()
384
385    intermediate_dict = {}
386
387    for vertex in vertices:
388        lorentz_names = _get_lorentz_names(vertex)
389        for v_name in lorentz_names:
390            if v_name in [('SigmaXY', 'SigmaZT'),
391                          ('SigmaXT', 'SigmaYZ'),
392                          ('SigmaYZ', 'SigmaXT'),
393                          ('SigmaZT', 'SigmaXY')]:
394                sign = -1
395            else:
396                sign = 1
397            if vertex not in intermediate_dict:
398                intermediate_dict[vertex] = sign * np.array(corr_data[v_name])
399            else:
400                intermediate_dict[vertex] += sign * np.array(corr_data[v_name])
401
402    result_dict = {}
403
404    for key, data in intermediate_dict.items():
405
406        rolled_array = np.moveaxis(data, 0, 8)
407
408        matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
409        for index in np.ndindex(rolled_array.shape[:-1]):
410            real = Obs([rolled_array[index].real], [ens_id], idl=[idx])
411            imag = Obs([rolled_array[index].imag], [ens_id], idl=[idx])
412            matrix[index] = CObs(real, imag)
413
414        result_dict[key] = Npr_matrix(matrix, mom_in=mom_in, mom_out=mom_out)
415
416    return result_dict
417
418
419def _get_lorentz_names(name):
420    lorentz_index = ['X', 'Y', 'Z', 'T']
421
422    res = []
423
424    if name == "TT":
425        for i in range(4):
426            for j in range(i + 1, 4):
427                res.append(("Sigma" + lorentz_index[i] + lorentz_index[j], "Sigma" + lorentz_index[i] + lorentz_index[j]))
428        return res
429
430    if name == "TTtilde":
431        for i in range(4):
432            for j in range(i + 1, 4):
433                for k in range(4):
434                    for o in range(k + 1, 4):
435                        fac = epsilon_tensor_rank4(i, j, k, o)
436                        if not np.isclose(fac, 0.0):
437                            res.append(("Sigma" + lorentz_index[i] + lorentz_index[j], "Sigma" + lorentz_index[k] + lorentz_index[o]))
438        return res
439
440    assert len(name) == 2
441
442    if 'S' in name or 'P' in name:
443        if not set(name) <= set(['S', 'P']):
444            raise Exception("'" + name + "' is not a Lorentz scalar")
445
446        g_names = {'S': 'Identity',
447                   'P': 'Gamma5'}
448
449        res.append((g_names[name[0]], g_names[name[1]]))
450
451    else:
452        if not set(name) <= set(['V', 'A']):
453            raise Exception("'" + name + "' is not a Lorentz scalar")
454
455        for ind in lorentz_index:
456            res.append(('Gamma' + ind + (name[0] == 'A') * 'Gamma5',
457                        'Gamma' + ind + (name[1] == 'A') * 'Gamma5'))
458
459    return res
def read_meson_hd5(path, filestem, ens_id, meson='meson_0', idl=None, gammas=None):
 59def read_meson_hd5(path, filestem, ens_id, meson='meson_0', idl=None, gammas=None):
 60    r'''Read hadrons meson hdf5 file and extract the meson labeled 'meson'
 61
 62    Parameters
 63    -----------------
 64    path : str
 65        path to the files to read
 66    filestem : str
 67        namestem of the files to read
 68    ens_id : str
 69        name of the ensemble, required for internal bookkeeping
 70    meson : str
 71        label of the meson to be extracted, standard value meson_0 which
 72        corresponds to the pseudoscalar pseudoscalar two-point function.
 73    gammas : tuple of strings
 74        Instrad of a meson label one can also provide a tuple of two strings
 75        indicating the gamma matrices at source and sink.
 76        ("Gamma5", "Gamma5") corresponds to the pseudoscalar pseudoscalar
 77        two-point function. The gammas argument dominateds over meson.
 78    idl : range
 79        If specified only configurations in the given range are read in.
 80    '''
 81
 82    files, idx = _get_files(path, filestem, idl)
 83
 84    tree = meson.rsplit('_')[0]
 85    if gammas is not None:
 86        h5file = h5py.File(path + '/' + files[0], "r")
 87        found_meson = None
 88        for key in h5file[tree].keys():
 89            if gammas[0] == h5file[tree][key].attrs["gamma_snk"][0].decode() and h5file[tree][key].attrs["gamma_src"][0].decode() == gammas[1]:
 90                found_meson = key
 91                break
 92        h5file.close()
 93        if found_meson:
 94            meson = found_meson
 95        else:
 96            raise Exception("Source Sink combination " + str(gammas) + " not found.")
 97
 98    corr_data = []
 99    infos = []
100    for hd5_file in files:
101        h5file = h5py.File(path + '/' + hd5_file, "r")
102        if not tree + '/' + meson in h5file:
103            raise Exception("Entry '" + meson + "' not contained in the files.")
104        raw_data = h5file[tree + '/' + meson + '/corr']
105        real_data = raw_data[:]["re"].astype(np.double)
106        corr_data.append(real_data)
107        if not infos:
108            for k, i in h5file[tree + '/' + meson].attrs.items():
109                infos.append(k + ': ' + i[0].decode())
110        h5file.close()
111    corr_data = np.array(corr_data)
112
113    l_obs = []
114    for c in corr_data.T:
115        l_obs.append(Obs([c], [ens_id], idl=[idx]))
116
117    corr = Corr(l_obs)
118    corr.tag = r", ".join(infos)
119    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.
  • gammas (tuple of strings): Instrad of a meson label one can also provide a tuple of two strings indicating the gamma matrices at source and sink. ("Gamma5", "Gamma5") corresponds to the pseudoscalar pseudoscalar two-point function. The gammas argument dominateds over meson.
  • idl (range): If specified only configurations in the given range are read in.
def read_DistillationContraction_hd5(path, ens_id, diagrams=['direct'], idl=None):
122def read_DistillationContraction_hd5(path, ens_id, diagrams=["direct"], idl=None):
123    """Read hadrons DistillationContraction hdf5 files in given directory structure
124
125    Parameters
126    -----------------
127    path : str
128        path to the directories to read
129    ens_id : str
130        name of the ensemble, required for internal bookkeeping
131    diagrams : list
132        List of strings of the diagrams to extract, e.g. ["direct", "box", "cross"].
133    idl : range
134        If specified only configurations in the given range are read in.
135    """
136
137    res_dict = {}
138
139    directories, idx = _get_files(path, "data", idl)
140
141    explore_path = Path(path + "/" + directories[0])
142
143    for explore_file in explore_path.iterdir():
144        if explore_file.is_file():
145            stem = explore_file.with_suffix("").with_suffix("").as_posix().split("/")[-1]
146        else:
147            continue
148
149        file_list = []
150        for dir in directories:
151            tmp_path = Path(path + "/" + dir)
152            file_list.append((tmp_path / stem).as_posix() + tmp_path.suffix + ".h5")
153
154        corr_data = {}
155
156        for diagram in diagrams:
157            corr_data[diagram] = []
158
159        for n_file, (hd5_file, n_traj) in enumerate(zip(file_list, list(idx))):
160            h5file = h5py.File(hd5_file)
161
162            if n_file == 0:
163                if h5file["DistillationContraction/Metadata"].attrs.get("TimeSources")[0].decode() != "0...":
164                    raise Exception("Routine is only implemented for files containing inversions on all timeslices.")
165
166                Nt = h5file["DistillationContraction/Metadata"].attrs.get("Nt")[0]
167
168                identifier = []
169                for in_file in range(len(h5file["DistillationContraction/Metadata/DmfInputFiles"].attrs.keys()) - 1):
170                    encoded_info = h5file["DistillationContraction/Metadata/DmfInputFiles"].attrs.get("DmfInputFiles_" + str(in_file))
171                    full_info = encoded_info[0].decode().split("/")[-1].replace(".h5", "").split("_")
172                    my_tuple = (full_info[0], full_info[1][1:], full_info[2], full_info[3])
173                    identifier.append(my_tuple)
174                identifier = tuple(identifier)
175                # "DistillationContraction/Metadata/DmfSuffix" contains info about different quarks, irrelevant in the SU(3) case.
176
177            for diagram in diagrams:
178                real_data = np.zeros(Nt)
179                for x0 in range(Nt):
180                    raw_data = h5file["DistillationContraction/Correlators/" + diagram + "/" + str(x0)][:]["re"].astype(np.double)
181                    real_data += np.roll(raw_data, -x0)
182                real_data /= Nt
183
184                corr_data[diagram].append(real_data)
185            h5file.close()
186
187        res_dict[str(identifier)] = {}
188
189        for diagram in diagrams:
190
191            tmp_data = np.array(corr_data[diagram])
192
193            l_obs = []
194            for c in tmp_data.T:
195                l_obs.append(Obs([c], [ens_id], idl=[idx]))
196
197            corr = Corr(l_obs)
198            corr.tag = str(identifier)
199
200            res_dict[str(identifier)][diagram] = corr
201
202    return res_dict

Read hadrons DistillationContraction hdf5 files in given directory structure

Parameters
  • path (str): path to the directories to read
  • ens_id (str): name of the ensemble, required for internal bookkeeping
  • diagrams (list): List of strings of the diagrams to extract, e.g. ["direct", "box", "cross"].
  • idl (range): If specified only configurations in the given range are read in.
class Npr_matrix(numpy.ndarray):
205class Npr_matrix(np.ndarray):
206
207    def __new__(cls, input_array, mom_in=None, mom_out=None):
208        obj = np.asarray(input_array).view(cls)
209        obj.mom_in = mom_in
210        obj.mom_out = mom_out
211        return obj
212
213    @property
214    def g5H(self):
215        """Gamma_5 hermitean conjugate
216
217        Uses the fact that the propagator is gamma5 hermitean, so just the
218        in and out momenta of the propagator are exchanged.
219        """
220        return Npr_matrix(self,
221                          mom_in=self.mom_out,
222                          mom_out=self.mom_in)
223
224    def _propagate_mom(self, other, name):
225        s_mom = getattr(self, name, None)
226        o_mom = getattr(other, name, None)
227        if s_mom is not None and o_mom is not None:
228            if not np.allclose(s_mom, o_mom):
229                raise Exception(name + ' does not match.')
230        return o_mom if o_mom is not None else s_mom
231
232    def __matmul__(self, other):
233        return self.__new__(Npr_matrix,
234                            super().__matmul__(other),
235                            self._propagate_mom(other, 'mom_in'),
236                            self._propagate_mom(other, 'mom_out'))
237
238    def __array_finalize__(self, obj):
239        if obj is None:
240            return
241        self.mom_in = getattr(obj, 'mom_in', None)
242        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):
245def read_ExternalLeg_hd5(path, filestem, ens_id, idl=None):
246    """Read hadrons ExternalLeg hdf5 file and output an array of CObs
247
248    Parameters
249    ----------
250    path : str
251        path to the files to read
252    filestem : str
253        namestem of the files to read
254    ens_id : str
255        name of the ensemble, required for internal bookkeeping
256    idl : range
257        If specified only configurations in the given range are read in.
258    """
259
260    files, idx = _get_files(path, filestem, idl)
261
262    mom = None
263
264    corr_data = []
265    for hd5_file in files:
266        file = h5py.File(path + '/' + hd5_file, "r")
267        raw_data = file['ExternalLeg/corr'][0][0].view('complex')
268        corr_data.append(raw_data)
269        if mom is None:
270            mom = np.array(str(file['ExternalLeg/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float)
271        file.close()
272    corr_data = np.array(corr_data)
273
274    rolled_array = np.rollaxis(corr_data, 0, 5)
275
276    matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
277    for si, sj, ci, cj in np.ndindex(rolled_array.shape[:-1]):
278        real = Obs([rolled_array[si, sj, ci, cj].real], [ens_id], idl=[idx])
279        imag = Obs([rolled_array[si, sj, ci, cj].imag], [ens_id], idl=[idx])
280        matrix[si, sj, ci, cj] = CObs(real, imag)
281
282    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):
285def read_Bilinear_hd5(path, filestem, ens_id, idl=None):
286    """Read hadrons Bilinear hdf5 file and output an array of CObs
287
288    Parameters
289    ----------
290    path : str
291        path to the files to read
292    filestem : str
293        namestem of the files to read
294    ens_id : str
295        name of the ensemble, required for internal bookkeeping
296    idl : range
297        If specified only configurations in the given range are read in.
298    """
299
300    files, idx = _get_files(path, filestem, idl)
301
302    mom_in = None
303    mom_out = None
304
305    corr_data = {}
306    for hd5_file in files:
307        file = h5py.File(path + '/' + hd5_file, "r")
308        for i in range(16):
309            name = file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['gamma'][0].decode('UTF-8')
310            if name not in corr_data:
311                corr_data[name] = []
312            raw_data = file['Bilinear/Bilinear_' + str(i) + '/corr'][0][0].view('complex')
313            corr_data[name].append(raw_data)
314            if mom_in is None:
315                mom_in = np.array(str(file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float)
316            if mom_out is None:
317                mom_out = np.array(str(file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['pOut'])[3:-2].strip().split(), dtype=float)
318
319        file.close()
320
321    result_dict = {}
322
323    for key, data in corr_data.items():
324        local_data = np.array(data)
325
326        rolled_array = np.rollaxis(local_data, 0, 5)
327
328        matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
329        for si, sj, ci, cj in np.ndindex(rolled_array.shape[:-1]):
330            real = Obs([rolled_array[si, sj, ci, cj].real], [ens_id], idl=[idx])
331            imag = Obs([rolled_array[si, sj, ci, cj].imag], [ens_id], idl=[idx])
332            matrix[si, sj, ci, cj] = CObs(real, imag)
333
334        result_dict[key] = Npr_matrix(matrix, mom_in=mom_in, mom_out=mom_out)
335
336    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']):
339def read_Fourquark_hd5(path, filestem, ens_id, idl=None, vertices=["VA", "AV"]):
340    """Read hadrons FourquarkFullyConnected hdf5 file and output an array of CObs
341
342    Parameters
343    ----------
344    path : str
345        path to the files to read
346    filestem : str
347        namestem of the files to read
348    ens_id : str
349        name of the ensemble, required for internal bookkeeping
350    idl : range
351        If specified only configurations in the given range are read in.
352    vertices : list
353        Vertex functions to be extracted.
354    """
355
356    files, idx = _get_files(path, filestem, idl)
357
358    mom_in = None
359    mom_out = None
360
361    vertex_names = []
362    for vertex in vertices:
363        vertex_names += _get_lorentz_names(vertex)
364
365    corr_data = {}
366
367    tree = 'FourQuarkFullyConnected/FourQuarkFullyConnected_'
368
369    for hd5_file in files:
370        file = h5py.File(path + '/' + hd5_file, "r")
371
372        for i in range(32):
373            name = (file[tree + str(i) + '/info'].attrs['gammaA'][0].decode('UTF-8'), file[tree + str(i) + '/info'].attrs['gammaB'][0].decode('UTF-8'))
374            if name in vertex_names:
375                if name not in corr_data:
376                    corr_data[name] = []
377                raw_data = file[tree + str(i) + '/corr'][0][0].view('complex')
378                corr_data[name].append(raw_data)
379                if mom_in is None:
380                    mom_in = np.array(str(file[tree + str(i) + '/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float)
381                if mom_out is None:
382                    mom_out = np.array(str(file[tree + str(i) + '/info'].attrs['pOut'])[3:-2].strip().split(), dtype=float)
383
384        file.close()
385
386    intermediate_dict = {}
387
388    for vertex in vertices:
389        lorentz_names = _get_lorentz_names(vertex)
390        for v_name in lorentz_names:
391            if v_name in [('SigmaXY', 'SigmaZT'),
392                          ('SigmaXT', 'SigmaYZ'),
393                          ('SigmaYZ', 'SigmaXT'),
394                          ('SigmaZT', 'SigmaXY')]:
395                sign = -1
396            else:
397                sign = 1
398            if vertex not in intermediate_dict:
399                intermediate_dict[vertex] = sign * np.array(corr_data[v_name])
400            else:
401                intermediate_dict[vertex] += sign * np.array(corr_data[v_name])
402
403    result_dict = {}
404
405    for key, data in intermediate_dict.items():
406
407        rolled_array = np.moveaxis(data, 0, 8)
408
409        matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
410        for index in np.ndindex(rolled_array.shape[:-1]):
411            real = Obs([rolled_array[index].real], [ens_id], idl=[idx])
412            imag = Obs([rolled_array[index].imag], [ens_id], idl=[idx])
413            matrix[index] = CObs(real, imag)
414
415        result_dict[key] = Npr_matrix(matrix, mom_in=mom_in, mom_out=mom_out)
416
417    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.