# pyerrors.input.hadrons

View Source

      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 source and sink.
     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

def read_meson_hd5(path, filestem, ens_id, meson='meson_0', idl=None, gammas=None): View Source

     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    Returns
     81    -------
     82    corr : Corr
     83        Correlator of the source sink combination in question.
     84    '''
     85
     86    files, idx = _get_files(path, filestem, idl)
     87
     88    tree = meson.rsplit('_')[0]
     89    if gammas is not None:
     90        h5file = h5py.File(path + '/' + files[0], "r")
     91        found_meson = None
     92        for key in h5file[tree].keys():
     93            if gammas[0] == h5file[tree][key].attrs["gamma_snk"][0].decode() and h5file[tree][key].attrs["gamma_src"][0].decode() == gammas[1]:
     94                found_meson = key
     95                break
     96        h5file.close()
     97        if found_meson:
     98            meson = found_meson
     99        else:
    100            raise Exception("Source Sink combination " + str(gammas) + " not found.")
    101
    102    corr_data = []
    103    infos = []
    104    for hd5_file in files:
    105        h5file = h5py.File(path + '/' + hd5_file, "r")
    106        if not tree + '/' + meson in h5file:
    107            raise Exception("Entry '" + meson + "' not contained in the files.")
    108        raw_data = h5file[tree + '/' + meson + '/corr']
    109        real_data = raw_data[:]["re"].astype(np.double)
    110        corr_data.append(real_data)
    111        if not infos:
    112            for k, i in h5file[tree + '/' + meson].attrs.items():
    113                infos.append(k + ': ' + i[0].decode())
    114        h5file.close()
    115    corr_data = np.array(corr_data)
    116
    117    l_obs = []
    118    for c in corr_data.T:
    119        l_obs.append(Obs([c], [ens_id], idl=[idx]))
    120
    121    corr = Corr(l_obs)
    122    corr.tag = r", ".join(infos)
    123    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.

###### Returns

- corr (Corr): Correlator of the source sink combination in question.

def extract_t0_hd5( path, filestem, ens_id, obs='Clover energy density', fit_range=5, idl=None, **kwargs): View Source

    144def extract_t0_hd5(path, filestem, ens_id, obs='Clover energy density', fit_range=5, idl=None, **kwargs):
    145    r'''Read hadrons FlowObservables hdf5 file and extract t0
    146
    147    Parameters
    148    -----------------
    149    path : str
    150        path to the files to read
    151    filestem : str
    152        namestem of the files to read
    153    ens_id : str
    154        name of the ensemble, required for internal bookkeeping
    155    obs : str
    156        label of the observable from which t0 should be extracted.
    157        Options: 'Clover energy density' and 'Plaquette energy density'
    158    fit_range : int
    159        Number of data points left and right of the zero
    160        crossing to be included in the linear fit. (Default: 5)
    161    idl : range
    162        If specified only configurations in the given range are read in.
    163    plot_fit : bool
    164        If true, the fit for the extraction of t0 is shown together with the data.
    165    '''
    166
    167    files, idx = _get_files(path, filestem, idl)
    168    tree = "FlowObservables"
    169
    170    h5file = h5py.File(path + '/' + files[0], "r")
    171    obs_key = None
    172    for key in h5file[tree].keys():
    173        if obs == h5file[tree][key].attrs["description"][0].decode():
    174            obs_key = key
    175            break
    176    h5file.close()
    177    if obs_key is None:
    178        raise Exception(f"Observable {obs} not found.")
    179
    180    corr_data = _extract_real_arrays(path, files, tree, ["FlowObservables_0", obs_key])
    181
    182    if not np.allclose(corr_data["FlowObservables_0"][0], corr_data["FlowObservables_0"][:]):
    183        raise Exception("Not all flow times were equal.")
    184
    185    t2E_dict = {}
    186    for t2, dat in zip(corr_data["FlowObservables_0"][0], corr_data[obs_key].T):
    187        t2E_dict[t2] = Obs([dat], [ens_id], idl=[idx]) - 0.3
    188
    189    return fit_t0(t2E_dict, fit_range, plot_fit=kwargs.get('plot_fit'))

Read hadrons FlowObservables hdf5 file and extract t0

###### 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
- obs (str): label of the observable from which t0 should be extracted. Options: 'Clover energy density' and 'Plaquette energy density'
- fit_range (int): Number of data points left and right of the zero crossing to be included in the linear fit. (Default: 5)
- idl (range): If specified only configurations in the given range are read in.
- plot_fit (bool): If true, the fit for the extraction of t0 is shown together with the data.

def read_DistillationContraction_hd5(path, ens_id, diagrams=['direct'], idl=None): View Source

    192def read_DistillationContraction_hd5(path, ens_id, diagrams=["direct"], idl=None):
    193    """Read hadrons DistillationContraction hdf5 files in given directory structure
    194
    195    Parameters
    196    -----------------
    197    path : str
    198        path to the directories to read
    199    ens_id : str
    200        name of the ensemble, required for internal bookkeeping
    201    diagrams : list
    202        List of strings of the diagrams to extract, e.g. ["direct", "box", "cross"].
    203    idl : range
    204        If specified only configurations in the given range are read in.
    205
    206    Returns
    207    -------
    208    result : dict
    209        extracted DistillationContration data
    210    """
    211
    212    res_dict = {}
    213
    214    directories, idx = _get_files(path, "data", idl)
    215
    216    explore_path = Path(path + "/" + directories[0])
    217
    218    for explore_file in explore_path.iterdir():
    219        if explore_file.is_file():
    220            stem = explore_file.with_suffix("").with_suffix("").as_posix().split("/")[-1]
    221        else:
    222            continue
    223
    224        file_list = []
    225        for dir in directories:
    226            tmp_path = Path(path + "/" + dir)
    227            file_list.append((tmp_path / stem).as_posix() + tmp_path.suffix + ".h5")
    228
    229        corr_data = {}
    230
    231        for diagram in diagrams:
    232            corr_data[diagram] = []
    233
    234        try:
    235            for n_file, (hd5_file, n_traj) in enumerate(zip(file_list, list(idx))):
    236                h5file = h5py.File(hd5_file)
    237
    238                if n_file == 0:
    239                    if h5file["DistillationContraction/Metadata"].attrs.get("TimeSources")[0].decode() != "0...":
    240                        raise Exception("Routine is only implemented for files containing inversions on all timeslices.")
    241
    242                    Nt = h5file["DistillationContraction/Metadata"].attrs.get("Nt")[0]
    243
    244                    identifier = []
    245                    for in_file in range(len(h5file["DistillationContraction/Metadata/DmfInputFiles"].attrs.keys()) - 1):
    246                        encoded_info = h5file["DistillationContraction/Metadata/DmfInputFiles"].attrs.get("DmfInputFiles_" + str(in_file))
    247                        full_info = encoded_info[0].decode().split("/")[-1].replace(".h5", "").split("_")
    248                        my_tuple = (full_info[0], full_info[1][1:], full_info[2], full_info[3])
    249                        identifier.append(my_tuple)
    250                    identifier = tuple(identifier)
    251                    # "DistillationContraction/Metadata/DmfSuffix" contains info about different quarks, irrelevant in the SU(3) case.
    252
    253                for diagram in diagrams:
    254
    255                    if diagram == "triangle" and "Identity" not in str(identifier):
    256                        part = "im"
    257                    else:
    258                        part = "re"
    259
    260                    real_data = np.zeros(Nt)
    261                    for x0 in range(Nt):
    262                        raw_data = h5file["DistillationContraction/Correlators/" + diagram + "/" + str(x0)][:][part].astype(np.double)
    263                        real_data += np.roll(raw_data, -x0)
    264                    real_data /= Nt
    265
    266                    corr_data[diagram].append(real_data)
    267                h5file.close()
    268
    269            res_dict[str(identifier)] = {}
    270
    271            for diagram in diagrams:
    272
    273                tmp_data = np.array(corr_data[diagram])
    274
    275                l_obs = []
    276                for c in tmp_data.T:
    277                    l_obs.append(Obs([c], [ens_id], idl=[idx]))
    278
    279                corr = Corr(l_obs)
    280                corr.tag = str(identifier)
    281
    282                res_dict[str(identifier)][diagram] = corr
    283        except FileNotFoundError:
    284            print("Skip", stem)
    285
    286    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.

###### Returns

- result (dict): extracted DistillationContration data

class Npr_matrix(numpy.ndarray): View Source

    289class Npr_matrix(np.ndarray):
    290
    291    def __new__(cls, input_array, mom_in=None, mom_out=None):
    292        obj = np.asarray(input_array).view(cls)
    293        obj.mom_in = mom_in
    294        obj.mom_out = mom_out
    295        return obj
    296
    297    @property
    298    def g5H(self):
    299        """Gamma_5 hermitean conjugate
    300
    301        Uses the fact that the propagator is gamma5 hermitean, so just the
    302        in and out momenta of the propagator are exchanged.
    303        """
    304        return Npr_matrix(self,
    305                          mom_in=self.mom_out,
    306                          mom_out=self.mom_in)
    307
    308    def _propagate_mom(self, other, name):
    309        s_mom = getattr(self, name, None)
    310        o_mom = getattr(other, name, None)
    311        if s_mom is not None and o_mom is not None:
    312            if not np.allclose(s_mom, o_mom):
    313                raise Exception(name + ' does not match.')
    314        return o_mom if o_mom is not None else s_mom
    315
    316    def __matmul__(self, other):
    317        return self.__new__(Npr_matrix,
    318                            super().__matmul__(other),
    319                            self._propagate_mom(other, 'mom_in'),
    320                            self._propagate_mom(other, 'mom_out'))
    321
    322    def __array_finalize__(self, obj):
    323        if obj is None:
    324            return
    325        self.mom_in = getattr(obj, 'mom_in', None)
    326        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])

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

    329def read_ExternalLeg_hd5(path, filestem, ens_id, idl=None):
    330    """Read hadrons ExternalLeg hdf5 file and output an array of CObs
    331
    332    Parameters
    333    ----------
    334    path : str
    335        path to the files to read
    336    filestem : str
    337        namestem of the files to read
    338    ens_id : str
    339        name of the ensemble, required for internal bookkeeping
    340    idl : range
    341        If specified only configurations in the given range are read in.
    342
    343    Returns
    344    -------
    345    result : Npr_matrix
    346        read Cobs-matrix
    347    """
    348
    349    files, idx = _get_files(path, filestem, idl)
    350
    351    mom = None
    352
    353    corr_data = []
    354    for hd5_file in files:
    355        file = h5py.File(path + '/' + hd5_file, "r")
    356        raw_data = file['ExternalLeg/corr'][0][0].view('complex')
    357        corr_data.append(raw_data)
    358        if mom is None:
    359            mom = np.array(str(file['ExternalLeg/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float)
    360        file.close()
    361    corr_data = np.array(corr_data)
    362
    363    rolled_array = np.rollaxis(corr_data, 0, 5)
    364
    365    matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
    366    for si, sj, ci, cj in np.ndindex(rolled_array.shape[:-1]):
    367        real = Obs([rolled_array[si, sj, ci, cj].real], [ens_id], idl=[idx])
    368        imag = Obs([rolled_array[si, sj, ci, cj].imag], [ens_id], idl=[idx])
    369        matrix[si, sj, ci, cj] = CObs(real, imag)
    370
    371    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.

###### Returns

- result (Npr_matrix): read Cobs-matrix

def read_Bilinear_hd5(path, filestem, ens_id, idl=None): View Source

    374def read_Bilinear_hd5(path, filestem, ens_id, idl=None):
    375    """Read hadrons Bilinear hdf5 file and output an array of CObs
    376
    377    Parameters
    378    ----------
    379    path : str
    380        path to the files to read
    381    filestem : str
    382        namestem of the files to read
    383    ens_id : str
    384        name of the ensemble, required for internal bookkeeping
    385    idl : range
    386        If specified only configurations in the given range are read in.
    387
    388    Returns
    389    -------
    390    result_dict: dict[Npr_matrix]
    391        extracted Bilinears
    392    """
    393
    394    files, idx = _get_files(path, filestem, idl)
    395
    396    mom_in = None
    397    mom_out = None
    398
    399    corr_data = {}
    400    for hd5_file in files:
    401        file = h5py.File(path + '/' + hd5_file, "r")
    402        for i in range(16):
    403            name = file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['gamma'][0].decode('UTF-8')
    404            if name not in corr_data:
    405                corr_data[name] = []
    406            raw_data = file['Bilinear/Bilinear_' + str(i) + '/corr'][0][0].view('complex')
    407            corr_data[name].append(raw_data)
    408            if mom_in is None:
    409                mom_in = np.array(str(file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float)
    410            if mom_out is None:
    411                mom_out = np.array(str(file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['pOut'])[3:-2].strip().split(), dtype=float)
    412
    413        file.close()
    414
    415    result_dict = {}
    416
    417    for key, data in corr_data.items():
    418        local_data = np.array(data)
    419
    420        rolled_array = np.rollaxis(local_data, 0, 5)
    421
    422        matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
    423        for si, sj, ci, cj in np.ndindex(rolled_array.shape[:-1]):
    424            real = Obs([rolled_array[si, sj, ci, cj].real], [ens_id], idl=[idx])
    425            imag = Obs([rolled_array[si, sj, ci, cj].imag], [ens_id], idl=[idx])
    426            matrix[si, sj, ci, cj] = CObs(real, imag)
    427
    428        result_dict[key] = Npr_matrix(matrix, mom_in=mom_in, mom_out=mom_out)
    429
    430    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.

###### Returns

- result_dict (dict[Npr_matrix]): extracted Bilinears

def read_Fourquark_hd5(path, filestem, ens_id, idl=None, vertices=['VA', 'AV']): View Source

    433def read_Fourquark_hd5(path, filestem, ens_id, idl=None, vertices=["VA", "AV"]):
    434    """Read hadrons FourquarkFullyConnected hdf5 file and output an array of CObs
    435
    436    Parameters
    437    ----------
    438    path : str
    439        path to the files to read
    440    filestem : str
    441        namestem of the files to read
    442    ens_id : str
    443        name of the ensemble, required for internal bookkeeping
    444    idl : range
    445        If specified only configurations in the given range are read in.
    446    vertices : list
    447        Vertex functions to be extracted.
    448
    449    Returns
    450    -------
    451    result_dict : dict
    452        extracted fourquark matrizes
    453    """
    454
    455    files, idx = _get_files(path, filestem, idl)
    456
    457    mom_in = None
    458    mom_out = None
    459
    460    vertex_names = []
    461    for vertex in vertices:
    462        vertex_names += _get_lorentz_names(vertex)
    463
    464    corr_data = {}
    465
    466    tree = 'FourQuarkFullyConnected/FourQuarkFullyConnected_'
    467
    468    for hd5_file in files:
    469        file = h5py.File(path + '/' + hd5_file, "r")
    470
    471        for i in range(32):
    472            name = (file[tree + str(i) + '/info'].attrs['gammaA'][0].decode('UTF-8'), file[tree + str(i) + '/info'].attrs['gammaB'][0].decode('UTF-8'))
    473            if name in vertex_names:
    474                if name not in corr_data:
    475                    corr_data[name] = []
    476                raw_data = file[tree + str(i) + '/corr'][0][0].view('complex')
    477                corr_data[name].append(raw_data)
    478                if mom_in is None:
    479                    mom_in = np.array(str(file[tree + str(i) + '/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float)
    480                if mom_out is None:
    481                    mom_out = np.array(str(file[tree + str(i) + '/info'].attrs['pOut'])[3:-2].strip().split(), dtype=float)
    482
    483        file.close()
    484
    485    intermediate_dict = {}
    486
    487    for vertex in vertices:
    488        lorentz_names = _get_lorentz_names(vertex)
    489        for v_name in lorentz_names:
    490            if v_name in [('SigmaXY', 'SigmaZT'),
    491                          ('SigmaXT', 'SigmaYZ'),
    492                          ('SigmaYZ', 'SigmaXT'),
    493                          ('SigmaZT', 'SigmaXY')]:
    494                sign = -1
    495            else:
    496                sign = 1
    497            if vertex not in intermediate_dict:
    498                intermediate_dict[vertex] = sign * np.array(corr_data[v_name])
    499            else:
    500                intermediate_dict[vertex] += sign * np.array(corr_data[v_name])
    501
    502    result_dict = {}
    503
    504    for key, data in intermediate_dict.items():
    505
    506        rolled_array = np.moveaxis(data, 0, 8)
    507
    508        matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
    509        for index in np.ndindex(rolled_array.shape[:-1]):
    510            real = Obs([rolled_array[index].real], [ens_id], idl=[idx])
    511            imag = Obs([rolled_array[index].imag], [ens_id], idl=[idx])
    512            matrix[index] = CObs(real, imag)
    513
    514        result_dict[key] = Npr_matrix(matrix, mom_in=mom_in, mom_out=mom_out)
    515
    516    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.

###### Returns

- result_dict (dict): extracted fourquark matrizes