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