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.")
 54
 55    return filtered_files, idx
 56
 57
 58def read_meson_hd5(path, filestem, ens_id, meson='meson_0', idl=None):
 59    """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    idl : range
 73        If specified only configurations in the given range are read in.
 74    """
 75
 76    files, idx = _get_files(path, filestem, idl)
 77
 78    tree = meson.rsplit('_')[0]
 79    corr_data = []
 80    infos = []
 81    for hd5_file in files:
 82        h5file = h5py.File(path + '/' + hd5_file, "r")
 83        if not tree + '/' + meson in h5file:
 84            raise Exception("Entry '" + meson + "' not contained in the files.")
 85        raw_data = h5file[tree + '/' + meson + '/corr']
 86        real_data = raw_data[:]["re"].astype(np.double)
 87        corr_data.append(real_data)
 88        if not infos:
 89            for k, i in h5file[tree + '/' + meson].attrs.items():
 90                infos.append(k + ': ' + i[0].decode())
 91        h5file.close()
 92    corr_data = np.array(corr_data)
 93
 94    l_obs = []
 95    for c in corr_data.T:
 96        l_obs.append(Obs([c], [ens_id], idl=[idx]))
 97
 98    corr = Corr(l_obs)
 99    corr.tag = r", ".join(infos)
100    return corr
101
102
103def read_DistillationContraction_hd5(path, ens_id, diagrams=["direct"], idl=None):
104    """Read hadrons DistillationContraction hdf5 files in given directory structure
105
106    Parameters
107    -----------------
108    path : str
109        path to the directories to read
110    ens_id : str
111        name of the ensemble, required for internal bookkeeping
112    diagrams : list
113        List of strings of the diagrams to extract, e.g. ["direct", "box", "cross"].
114    idl : range
115        If specified only configurations in the given range are read in.
116    """
117
118    res_dict = {}
119
120    directories, idx = _get_files(path, "data", idl)
121
122    explore_path = Path(path + "/" + directories[0])
123
124    for explore_file in explore_path.iterdir():
125        if explore_file.is_file():
126            stem = explore_file.with_suffix("").with_suffix("").as_posix().split("/")[-1]
127
128        file_list = []
129        for dir in directories:
130            tmp_path = Path(path + "/" + dir)
131            file_list.append((tmp_path / stem).as_posix() + tmp_path.suffix + ".h5")
132
133        corr_data = {}
134
135        for diagram in diagrams:
136            corr_data[diagram] = []
137
138        for n_file, (hd5_file, n_traj) in enumerate(zip(file_list, list(idx))):
139            h5file = h5py.File(hd5_file)
140
141            if n_file == 0:
142                if h5file["DistillationContraction/Metadata"].attrs.get("TimeSources")[0].decode() != "0...":
143                    raise Exception("Routine is only implemented for files containing inversions on all timeslices.")
144
145                Nt = h5file["DistillationContraction/Metadata"].attrs.get("Nt")[0]
146
147                identifier = []
148                for in_file in range(4):
149                    encoded_info = h5file["DistillationContraction/Metadata/DmfInputFiles"].attrs.get("DmfInputFiles_" + str(in_file))
150                    full_info = encoded_info[0].decode().split("/")[-1].replace(".h5", "").split("_")
151                    my_tuple = (full_info[0], full_info[1][1:], full_info[2], full_info[3])
152                    identifier.append(my_tuple)
153                identifier = tuple(identifier)
154                # "DistillationContraction/Metadata/DmfSuffix" contains info about different quarks, irrelevant in the SU(3) case.
155
156            check_traj = h5file["DistillationContraction/Metadata"].attrs.get("Traj")[0]
157
158            assert check_traj == n_traj
159
160            for diagram in diagrams:
161                real_data = np.zeros(Nt)
162                for x0 in range(Nt):
163                    raw_data = h5file["DistillationContraction/Correlators/" + diagram + "/" + str(x0)]
164                    real_data += np.roll(raw_data[:]["re"].astype(np.double), -x0)
165                real_data /= Nt
166
167                corr_data[diagram].append(real_data)
168            h5file.close()
169
170        res_dict[str(identifier)] = {}
171
172        for diagram in diagrams:
173
174            tmp_data = np.array(corr_data[diagram])
175
176            l_obs = []
177            for c in tmp_data.T:
178                l_obs.append(Obs([c], [ens_id], idl=[idx]))
179
180            corr = Corr(l_obs)
181            corr.tag = str(identifier)
182
183            res_dict[str(identifier)][diagram] = corr
184
185    return res_dict
186
187
188class Npr_matrix(np.ndarray):
189
190    def __new__(cls, input_array, mom_in=None, mom_out=None):
191        obj = np.asarray(input_array).view(cls)
192        obj.mom_in = mom_in
193        obj.mom_out = mom_out
194        return obj
195
196    @property
197    def g5H(self):
198        """Gamma_5 hermitean conjugate
199
200        Uses the fact that the propagator is gamma5 hermitean, so just the
201        in and out momenta of the propagator are exchanged.
202        """
203        return Npr_matrix(self,
204                          mom_in=self.mom_out,
205                          mom_out=self.mom_in)
206
207    def _propagate_mom(self, other, name):
208        s_mom = getattr(self, name, None)
209        o_mom = getattr(other, name, None)
210        if s_mom is not None and o_mom is not None:
211            if not np.allclose(s_mom, o_mom):
212                raise Exception(name + ' does not match.')
213        return o_mom if o_mom is not None else s_mom
214
215    def __matmul__(self, other):
216        return self.__new__(Npr_matrix,
217                            super().__matmul__(other),
218                            self._propagate_mom(other, 'mom_in'),
219                            self._propagate_mom(other, 'mom_out'))
220
221    def __array_finalize__(self, obj):
222        if obj is None:
223            return
224        self.mom_in = getattr(obj, 'mom_in', None)
225        self.mom_out = getattr(obj, 'mom_out', None)
226
227
228def read_ExternalLeg_hd5(path, filestem, ens_id, idl=None):
229    """Read hadrons ExternalLeg hdf5 file and output an array of CObs
230
231    Parameters
232    ----------
233    path : str
234        path to the files to read
235    filestem : str
236        namestem of the files to read
237    ens_id : str
238        name of the ensemble, required for internal bookkeeping
239    idl : range
240        If specified only configurations in the given range are read in.
241    """
242
243    files, idx = _get_files(path, filestem, idl)
244
245    mom = None
246
247    corr_data = []
248    for hd5_file in files:
249        file = h5py.File(path + '/' + hd5_file, "r")
250        raw_data = file['ExternalLeg/corr'][0][0].view('complex')
251        corr_data.append(raw_data)
252        if mom is None:
253            mom = np.array(str(file['ExternalLeg/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float)
254        file.close()
255    corr_data = np.array(corr_data)
256
257    rolled_array = np.rollaxis(corr_data, 0, 5)
258
259    matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
260    for si, sj, ci, cj in np.ndindex(rolled_array.shape[:-1]):
261        real = Obs([rolled_array[si, sj, ci, cj].real], [ens_id], idl=[idx])
262        imag = Obs([rolled_array[si, sj, ci, cj].imag], [ens_id], idl=[idx])
263        matrix[si, sj, ci, cj] = CObs(real, imag)
264
265    return Npr_matrix(matrix, mom_in=mom)
266
267
268def read_Bilinear_hd5(path, filestem, ens_id, idl=None):
269    """Read hadrons Bilinear hdf5 file and output an array of CObs
270
271    Parameters
272    ----------
273    path : str
274        path to the files to read
275    filestem : str
276        namestem of the files to read
277    ens_id : str
278        name of the ensemble, required for internal bookkeeping
279    idl : range
280        If specified only configurations in the given range are read in.
281    """
282
283    files, idx = _get_files(path, filestem, idl)
284
285    mom_in = None
286    mom_out = None
287
288    corr_data = {}
289    for hd5_file in files:
290        file = h5py.File(path + '/' + hd5_file, "r")
291        for i in range(16):
292            name = file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['gamma'][0].decode('UTF-8')
293            if name not in corr_data:
294                corr_data[name] = []
295            raw_data = file['Bilinear/Bilinear_' + str(i) + '/corr'][0][0].view('complex')
296            corr_data[name].append(raw_data)
297            if mom_in is None:
298                mom_in = np.array(str(file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float)
299            if mom_out is None:
300                mom_out = np.array(str(file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['pOut'])[3:-2].strip().split(), dtype=float)
301
302        file.close()
303
304    result_dict = {}
305
306    for key, data in corr_data.items():
307        local_data = np.array(data)
308
309        rolled_array = np.rollaxis(local_data, 0, 5)
310
311        matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
312        for si, sj, ci, cj in np.ndindex(rolled_array.shape[:-1]):
313            real = Obs([rolled_array[si, sj, ci, cj].real], [ens_id], idl=[idx])
314            imag = Obs([rolled_array[si, sj, ci, cj].imag], [ens_id], idl=[idx])
315            matrix[si, sj, ci, cj] = CObs(real, imag)
316
317        result_dict[key] = Npr_matrix(matrix, mom_in=mom_in, mom_out=mom_out)
318
319    return result_dict
320
321
322def read_Fourquark_hd5(path, filestem, ens_id, idl=None, vertices=["VA", "AV"]):
323    """Read hadrons FourquarkFullyConnected hdf5 file and output an array of CObs
324
325    Parameters
326    ----------
327    path : str
328        path to the files to read
329    filestem : str
330        namestem of the files to read
331    ens_id : str
332        name of the ensemble, required for internal bookkeeping
333    idl : range
334        If specified only configurations in the given range are read in.
335    vertices : list
336        Vertex functions to be extracted.
337    """
338
339    files, idx = _get_files(path, filestem, idl)
340
341    mom_in = None
342    mom_out = None
343
344    vertex_names = []
345    for vertex in vertices:
346        vertex_names += _get_lorentz_names(vertex)
347
348    corr_data = {}
349
350    tree = 'FourQuarkFullyConnected/FourQuarkFullyConnected_'
351
352    for hd5_file in files:
353        file = h5py.File(path + '/' + hd5_file, "r")
354
355        for i in range(32):
356            name = (file[tree + str(i) + '/info'].attrs['gammaA'][0].decode('UTF-8'), file[tree + str(i) + '/info'].attrs['gammaB'][0].decode('UTF-8'))
357            if name in vertex_names:
358                if name not in corr_data:
359                    corr_data[name] = []
360                raw_data = file[tree + str(i) + '/corr'][0][0].view('complex')
361                corr_data[name].append(raw_data)
362                if mom_in is None:
363                    mom_in = np.array(str(file[tree + str(i) + '/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float)
364                if mom_out is None:
365                    mom_out = np.array(str(file[tree + str(i) + '/info'].attrs['pOut'])[3:-2].strip().split(), dtype=float)
366
367        file.close()
368
369    intermediate_dict = {}
370
371    for vertex in vertices:
372        lorentz_names = _get_lorentz_names(vertex)
373        for v_name in lorentz_names:
374            if v_name in [('SigmaXY', 'SigmaZT'),
375                          ('SigmaXT', 'SigmaYZ'),
376                          ('SigmaYZ', 'SigmaXT'),
377                          ('SigmaZT', 'SigmaXY')]:
378                sign = -1
379            else:
380                sign = 1
381            if vertex not in intermediate_dict:
382                intermediate_dict[vertex] = sign * np.array(corr_data[v_name])
383            else:
384                intermediate_dict[vertex] += sign * np.array(corr_data[v_name])
385
386    result_dict = {}
387
388    for key, data in intermediate_dict.items():
389
390        rolled_array = np.moveaxis(data, 0, 8)
391
392        matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
393        for index in np.ndindex(rolled_array.shape[:-1]):
394            real = Obs([rolled_array[index].real], [ens_id], idl=[idx])
395            imag = Obs([rolled_array[index].imag], [ens_id], idl=[idx])
396            matrix[index] = CObs(real, imag)
397
398        result_dict[key] = Npr_matrix(matrix, mom_in=mom_in, mom_out=mom_out)
399
400    return result_dict
401
402
403def _get_lorentz_names(name):
404    lorentz_index = ['X', 'Y', 'Z', 'T']
405
406    res = []
407
408    if name == "TT":
409        for i in range(4):
410            for j in range(i + 1, 4):
411                res.append(("Sigma" + lorentz_index[i] + lorentz_index[j], "Sigma" + lorentz_index[i] + lorentz_index[j]))
412        return res
413
414    if name == "TTtilde":
415        for i in range(4):
416            for j in range(i + 1, 4):
417                for k in range(4):
418                    for o in range(k + 1, 4):
419                        fac = epsilon_tensor_rank4(i, j, k, o)
420                        if not np.isclose(fac, 0.0):
421                            res.append(("Sigma" + lorentz_index[i] + lorentz_index[j], "Sigma" + lorentz_index[k] + lorentz_index[o]))
422        return res
423
424    assert len(name) == 2
425
426    if 'S' in name or 'P' in name:
427        if not set(name) <= set(['S', 'P']):
428            raise Exception("'" + name + "' is not a Lorentz scalar")
429
430        g_names = {'S': 'Identity',
431                   'P': 'Gamma5'}
432
433        res.append((g_names[name[0]], g_names[name[1]]))
434
435    else:
436        if not set(name) <= set(['V', 'A']):
437            raise Exception("'" + name + "' is not a Lorentz scalar")
438
439        for ind in lorentz_index:
440            res.append(('Gamma' + ind + (name[0] == 'A') * 'Gamma5',
441                        'Gamma' + ind + (name[1] == 'A') * 'Gamma5'))
442
443    return res