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