pyerrors.input.hadrons

  1import os
  2from collections import Counter
  3import h5py
  4from pathlib import Path
  5import numpy as np
  6from ..obs import Obs, CObs
  7from ..correlators import Corr
  8from ..dirac import epsilon_tensor_rank4
  9from .misc import fit_t0
 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    else:
 52        raise Exception("Configurations are not evenly spaced. Provide an idl if you want to proceed with this set of configurations.")
 53
 54    return filtered_files, idx
 55
 56
 57def read_meson_hd5(path, filestem, ens_id, meson='meson_0', idl=None, gammas=None):
 58    r'''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    gammas : tuple of strings
 72        Instrad of a meson label one can also provide a tuple of two strings
 73        indicating the gamma matrices at sink and source (gamma_snk, gamma_src).
 74        ("Gamma5", "Gamma5") corresponds to the pseudoscalar pseudoscalar
 75        two-point function. The gammas argument dominateds over meson.
 76    idl : range
 77        If specified only configurations in the given range are read in.
 78
 79    Returns
 80    -------
 81    corr : Corr
 82        Correlator of the source sink combination in question.
 83    '''
 84
 85    files, idx = _get_files(path, filestem, idl)
 86
 87    tree = meson.rsplit('_')[0]
 88    if gammas is not None:
 89        h5file = h5py.File(path + '/' + files[0], "r")
 90        found_meson = None
 91        for key in h5file[tree].keys():
 92            if gammas[0] == h5file[tree][key].attrs["gamma_snk"][0].decode() and h5file[tree][key].attrs["gamma_src"][0].decode() == gammas[1]:
 93                found_meson = key
 94                break
 95        h5file.close()
 96        if found_meson:
 97            meson = found_meson
 98        else:
 99            raise Exception("Source Sink combination " + str(gammas) + " not found.")
100
101    corr_data = []
102    infos = []
103    for hd5_file in files:
104        h5file = h5py.File(path + '/' + hd5_file, "r")
105        if not tree + '/' + meson in h5file:
106            raise Exception("Entry '" + meson + "' not contained in the files.")
107        raw_data = h5file[tree + '/' + meson + '/corr']
108        real_data = raw_data[:]["re"].astype(np.double)
109        corr_data.append(real_data)
110        if not infos:
111            for k, i in h5file[tree + '/' + meson].attrs.items():
112                infos.append(k + ': ' + i[0].decode())
113        h5file.close()
114    corr_data = np.array(corr_data)
115
116    l_obs = []
117    for c in corr_data.T:
118        l_obs.append(Obs([c], [ens_id], idl=[idx]))
119
120    corr = Corr(l_obs)
121    corr.tag = r", ".join(infos)
122    return corr
123
124
125def _extract_real_arrays(path, files, tree, keys):
126    corr_data = {}
127    for key in keys:
128        corr_data[key] = []
129    for hd5_file in files:
130        h5file = h5py.File(path + '/' + hd5_file, "r")
131        for key in keys:
132            if not tree + '/' + key in h5file:
133                raise Exception("Entry '" + key + "' not contained in the files.")
134            raw_data = h5file[tree + '/' + key + '/data']
135            real_data = raw_data[:].astype(np.double)
136            corr_data[key].append(real_data)
137        h5file.close()
138    for key in keys:
139        corr_data[key] = np.array(corr_data[key])
140    return corr_data
141
142
143def extract_t0_hd5(path, filestem, ens_id, obs='Clover energy density', fit_range=5, idl=None, **kwargs):
144    r'''Read hadrons FlowObservables hdf5 file and extract t0
145
146    Parameters
147    -----------------
148    path : str
149        path to the files to read
150    filestem : str
151        namestem of the files to read
152    ens_id : str
153        name of the ensemble, required for internal bookkeeping
154    obs : str
155        label of the observable from which t0 should be extracted.
156        Options: 'Clover energy density' and 'Plaquette energy density'
157    fit_range : int
158        Number of data points left and right of the zero
159        crossing to be included in the linear fit. (Default: 5)
160    idl : range
161        If specified only configurations in the given range are read in.
162    plot_fit : bool
163        If true, the fit for the extraction of t0 is shown together with the data.
164    '''
165
166    files, idx = _get_files(path, filestem, idl)
167    tree = "FlowObservables"
168
169    h5file = h5py.File(path + '/' + files[0], "r")
170    obs_key = None
171    for key in h5file[tree].keys():
172        if obs == h5file[tree][key].attrs["description"][0].decode():
173            obs_key = key
174            break
175    h5file.close()
176    if obs_key is None:
177        raise Exception(f"Observable {obs} not found.")
178
179    corr_data = _extract_real_arrays(path, files, tree, ["FlowObservables_0", obs_key])
180
181    if not np.allclose(corr_data["FlowObservables_0"][0], corr_data["FlowObservables_0"][:]):
182        raise Exception("Not all flow times were equal.")
183
184    t2E_dict = {}
185    for t2, dat in zip(corr_data["FlowObservables_0"][0], corr_data[obs_key].T):
186        t2E_dict[t2] = Obs([dat], [ens_id], idl=[idx]) - 0.3
187
188    return fit_t0(t2E_dict, fit_range, plot_fit=kwargs.get('plot_fit'))
189
190
191def read_DistillationContraction_hd5(path, ens_id, diagrams=["direct"], idl=None):
192    """Read hadrons DistillationContraction hdf5 files in given directory structure
193
194    Parameters
195    -----------------
196    path : str
197        path to the directories to read
198    ens_id : str
199        name of the ensemble, required for internal bookkeeping
200    diagrams : list
201        List of strings of the diagrams to extract, e.g. ["direct", "box", "cross"].
202    idl : range
203        If specified only configurations in the given range are read in.
204
205    Returns
206    -------
207    result : dict
208        extracted DistillationContration data
209    """
210
211    res_dict = {}
212
213    directories, idx = _get_files(path, "data", idl)
214
215    explore_path = Path(path + "/" + directories[0])
216
217    for explore_file in explore_path.iterdir():
218        if explore_file.is_file():
219            stem = explore_file.with_suffix("").with_suffix("").as_posix().split("/")[-1]
220        else:
221            continue
222
223        file_list = []
224        for dir in directories:
225            tmp_path = Path(path + "/" + dir)
226            file_list.append((tmp_path / stem).as_posix() + tmp_path.suffix + ".h5")
227
228        corr_data = {}
229
230        for diagram in diagrams:
231            corr_data[diagram] = []
232
233        try:
234            for n_file, (hd5_file, n_traj) in enumerate(zip(file_list, list(idx))):
235                h5file = h5py.File(hd5_file)
236
237                if n_file == 0:
238                    if h5file["DistillationContraction/Metadata"].attrs.get("TimeSources")[0].decode() != "0...":
239                        raise Exception("Routine is only implemented for files containing inversions on all timeslices.")
240
241                    Nt = h5file["DistillationContraction/Metadata"].attrs.get("Nt")[0]
242
243                    identifier = []
244                    for in_file in range(len(h5file["DistillationContraction/Metadata/DmfInputFiles"].attrs.keys()) - 1):
245                        encoded_info = h5file["DistillationContraction/Metadata/DmfInputFiles"].attrs.get("DmfInputFiles_" + str(in_file))
246                        full_info = encoded_info[0].decode().split("/")[-1].replace(".h5", "").split("_")
247                        my_tuple = (full_info[0], full_info[1][1:], full_info[2], full_info[3])
248                        identifier.append(my_tuple)
249                    identifier = tuple(identifier)
250                    # "DistillationContraction/Metadata/DmfSuffix" contains info about different quarks, irrelevant in the SU(3) case.
251
252                for diagram in diagrams:
253
254                    if diagram == "triangle" and "Identity" not in str(identifier):
255                        part = "im"
256                    else:
257                        part = "re"
258
259                    real_data = np.zeros(Nt)
260                    for x0 in range(Nt):
261                        raw_data = h5file["DistillationContraction/Correlators/" + diagram + "/" + str(x0)][:][part].astype(np.double)
262                        real_data += np.roll(raw_data, -x0)
263                    real_data /= Nt
264
265                    corr_data[diagram].append(real_data)
266                h5file.close()
267
268            res_dict[str(identifier)] = {}
269
270            for diagram in diagrams:
271
272                tmp_data = np.array(corr_data[diagram])
273
274                l_obs = []
275                for c in tmp_data.T:
276                    l_obs.append(Obs([c], [ens_id], idl=[idx]))
277
278                corr = Corr(l_obs)
279                corr.tag = str(identifier)
280
281                res_dict[str(identifier)][diagram] = corr
282        except FileNotFoundError:
283            print("Skip", stem)
284
285    return res_dict
286
287
288class Npr_matrix(np.ndarray):
289
290    def __new__(cls, input_array, mom_in=None, mom_out=None):
291        obj = np.asarray(input_array).view(cls)
292        obj.mom_in = mom_in
293        obj.mom_out = mom_out
294        return obj
295
296    @property
297    def g5H(self):
298        """Gamma_5 hermitean conjugate
299
300        Uses the fact that the propagator is gamma5 hermitean, so just the
301        in and out momenta of the propagator are exchanged.
302        """
303        return Npr_matrix(self,
304                          mom_in=self.mom_out,
305                          mom_out=self.mom_in)
306
307    def _propagate_mom(self, other, name):
308        s_mom = getattr(self, name, None)
309        o_mom = getattr(other, name, None)
310        if s_mom is not None and o_mom is not None:
311            if not np.allclose(s_mom, o_mom):
312                raise Exception(name + ' does not match.')
313        return o_mom if o_mom is not None else s_mom
314
315    def __matmul__(self, other):
316        return self.__new__(Npr_matrix,
317                            super().__matmul__(other),
318                            self._propagate_mom(other, 'mom_in'),
319                            self._propagate_mom(other, 'mom_out'))
320
321    def __array_finalize__(self, obj):
322        if obj is None:
323            return
324        self.mom_in = getattr(obj, 'mom_in', None)
325        self.mom_out = getattr(obj, 'mom_out', None)
326
327
328def read_ExternalLeg_hd5(path, filestem, ens_id, idl=None):
329    """Read hadrons ExternalLeg hdf5 file and output an array of CObs
330
331    Parameters
332    ----------
333    path : str
334        path to the files to read
335    filestem : str
336        namestem of the files to read
337    ens_id : str
338        name of the ensemble, required for internal bookkeeping
339    idl : range
340        If specified only configurations in the given range are read in.
341
342    Returns
343    -------
344    result : Npr_matrix
345        read Cobs-matrix
346    """
347
348    files, idx = _get_files(path, filestem, idl)
349
350    mom = None
351
352    corr_data = []
353    for hd5_file in files:
354        file = h5py.File(path + '/' + hd5_file, "r")
355        raw_data = file['ExternalLeg/corr'][0][0].view('complex')
356        corr_data.append(raw_data)
357        if mom is None:
358            mom = np.array(str(file['ExternalLeg/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float)
359        file.close()
360    corr_data = np.array(corr_data)
361
362    rolled_array = np.rollaxis(corr_data, 0, 5)
363
364    matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
365    for si, sj, ci, cj in np.ndindex(rolled_array.shape[:-1]):
366        real = Obs([rolled_array[si, sj, ci, cj].real], [ens_id], idl=[idx])
367        imag = Obs([rolled_array[si, sj, ci, cj].imag], [ens_id], idl=[idx])
368        matrix[si, sj, ci, cj] = CObs(real, imag)
369
370    return Npr_matrix(matrix, mom_in=mom)
371
372
373def read_Bilinear_hd5(path, filestem, ens_id, idl=None):
374    """Read hadrons Bilinear hdf5 file and output an array of CObs
375
376    Parameters
377    ----------
378    path : str
379        path to the files to read
380    filestem : str
381        namestem of the files to read
382    ens_id : str
383        name of the ensemble, required for internal bookkeeping
384    idl : range
385        If specified only configurations in the given range are read in.
386
387    Returns
388    -------
389    result_dict: dict[Npr_matrix]
390        extracted Bilinears
391    """
392
393    files, idx = _get_files(path, filestem, idl)
394
395    mom_in = None
396    mom_out = None
397
398    corr_data = {}
399    for hd5_file in files:
400        file = h5py.File(path + '/' + hd5_file, "r")
401        for i in range(16):
402            name = file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['gamma'][0].decode('UTF-8')
403            if name not in corr_data:
404                corr_data[name] = []
405            raw_data = file['Bilinear/Bilinear_' + str(i) + '/corr'][0][0].view('complex')
406            corr_data[name].append(raw_data)
407            if mom_in is None:
408                mom_in = np.array(str(file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float)
409            if mom_out is None:
410                mom_out = np.array(str(file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['pOut'])[3:-2].strip().split(), dtype=float)
411
412        file.close()
413
414    result_dict = {}
415
416    for key, data in corr_data.items():
417        local_data = np.array(data)
418
419        rolled_array = np.rollaxis(local_data, 0, 5)
420
421        matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
422        for si, sj, ci, cj in np.ndindex(rolled_array.shape[:-1]):
423            real = Obs([rolled_array[si, sj, ci, cj].real], [ens_id], idl=[idx])
424            imag = Obs([rolled_array[si, sj, ci, cj].imag], [ens_id], idl=[idx])
425            matrix[si, sj, ci, cj] = CObs(real, imag)
426
427        result_dict[key] = Npr_matrix(matrix, mom_in=mom_in, mom_out=mom_out)
428
429    return result_dict
430
431
432def read_Fourquark_hd5(path, filestem, ens_id, idl=None, vertices=["VA", "AV"]):
433    """Read hadrons FourquarkFullyConnected hdf5 file and output an array of CObs
434
435    Parameters
436    ----------
437    path : str
438        path to the files to read
439    filestem : str
440        namestem of the files to read
441    ens_id : str
442        name of the ensemble, required for internal bookkeeping
443    idl : range
444        If specified only configurations in the given range are read in.
445    vertices : list
446        Vertex functions to be extracted.
447
448    Returns
449    -------
450    result_dict : dict
451        extracted fourquark matrizes
452    """
453
454    files, idx = _get_files(path, filestem, idl)
455
456    mom_in = None
457    mom_out = None
458
459    vertex_names = []
460    for vertex in vertices:
461        vertex_names += _get_lorentz_names(vertex)
462
463    corr_data = {}
464
465    tree = 'FourQuarkFullyConnected/FourQuarkFullyConnected_'
466
467    for hd5_file in files:
468        file = h5py.File(path + '/' + hd5_file, "r")
469
470        for i in range(32):
471            name = (file[tree + str(i) + '/info'].attrs['gammaA'][0].decode('UTF-8'), file[tree + str(i) + '/info'].attrs['gammaB'][0].decode('UTF-8'))
472            if name in vertex_names:
473                if name not in corr_data:
474                    corr_data[name] = []
475                raw_data = file[tree + str(i) + '/corr'][0][0].view('complex')
476                corr_data[name].append(raw_data)
477                if mom_in is None:
478                    mom_in = np.array(str(file[tree + str(i) + '/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float)
479                if mom_out is None:
480                    mom_out = np.array(str(file[tree + str(i) + '/info'].attrs['pOut'])[3:-2].strip().split(), dtype=float)
481
482        file.close()
483
484    intermediate_dict = {}
485
486    for vertex in vertices:
487        lorentz_names = _get_lorentz_names(vertex)
488        for v_name in lorentz_names:
489            if v_name in [('SigmaXY', 'SigmaZT'),
490                          ('SigmaXT', 'SigmaYZ'),
491                          ('SigmaYZ', 'SigmaXT'),
492                          ('SigmaZT', 'SigmaXY')]:
493                sign = -1
494            else:
495                sign = 1
496            if vertex not in intermediate_dict:
497                intermediate_dict[vertex] = sign * np.array(corr_data[v_name])
498            else:
499                intermediate_dict[vertex] += sign * np.array(corr_data[v_name])
500
501    result_dict = {}
502
503    for key, data in intermediate_dict.items():
504
505        rolled_array = np.moveaxis(data, 0, 8)
506
507        matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
508        for index in np.ndindex(rolled_array.shape[:-1]):
509            real = Obs([rolled_array[index].real], [ens_id], idl=[idx])
510            imag = Obs([rolled_array[index].imag], [ens_id], idl=[idx])
511            matrix[index] = CObs(real, imag)
512
513        result_dict[key] = Npr_matrix(matrix, mom_in=mom_in, mom_out=mom_out)
514
515    return result_dict
516
517
518def _get_lorentz_names(name):
519    lorentz_index = ['X', 'Y', 'Z', 'T']
520
521    res = []
522
523    if name == "TT":
524        for i in range(4):
525            for j in range(i + 1, 4):
526                res.append(("Sigma" + lorentz_index[i] + lorentz_index[j], "Sigma" + lorentz_index[i] + lorentz_index[j]))
527        return res
528
529    if name == "TTtilde":
530        for i in range(4):
531            for j in range(i + 1, 4):
532                for k in range(4):
533                    for o in range(k + 1, 4):
534                        fac = epsilon_tensor_rank4(i, j, k, o)
535                        if not np.isclose(fac, 0.0):
536                            res.append(("Sigma" + lorentz_index[i] + lorentz_index[j], "Sigma" + lorentz_index[k] + lorentz_index[o]))
537        return res
538
539    assert len(name) == 2
540
541    if 'S' in name or 'P' in name:
542        if not set(name) <= set(['S', 'P']):
543            raise Exception("'" + name + "' is not a Lorentz scalar")
544
545        g_names = {'S': 'Identity',
546                   'P': 'Gamma5'}
547
548        res.append((g_names[name[0]], g_names[name[1]]))
549
550    else:
551        if not set(name) <= set(['V', 'A']):
552            raise Exception("'" + name + "' is not a Lorentz scalar")
553
554        for ind in lorentz_index:
555            res.append(('Gamma' + ind + (name[0] == 'A') * 'Gamma5',
556                        'Gamma' + ind + (name[1] == 'A') * 'Gamma5'))
557
558    return res