pyerrors.input.misc

  1import os
  2import fnmatch
  3import re
  4import struct
  5import numpy as np  # Thinly-wrapped numpy
  6from ..obs import Obs
  7
  8
  9def read_pbp(path, prefix, **kwargs):
 10    """Read pbp format from given folder structure. Returns a list of length nrw
 11
 12    Parameters
 13    ----------
 14    r_start : list
 15        list which contains the first config to be read for each replicum
 16    r_stop : list
 17        list which contains the last config to be read for each replicum
 18    """
 19
 20    ls = []
 21    for (dirpath, dirnames, filenames) in os.walk(path):
 22        ls.extend(filenames)
 23        break
 24
 25    if not ls:
 26        raise Exception('Error, directory not found')
 27
 28    # Exclude files with different names
 29    for exc in ls:
 30        if not fnmatch.fnmatch(exc, prefix + '*.dat'):
 31            ls = list(set(ls) - set([exc]))
 32    if len(ls) > 1:
 33        ls.sort(key=lambda x: int(re.findall(r'\d+', x[len(prefix):])[0]))
 34    replica = len(ls)
 35
 36    if 'r_start' in kwargs:
 37        r_start = kwargs.get('r_start')
 38        if len(r_start) != replica:
 39            raise Exception('r_start does not match number of replicas')
 40        # Adjust Configuration numbering to python index
 41        r_start = [o - 1 if o else None for o in r_start]
 42    else:
 43        r_start = [None] * replica
 44
 45    if 'r_stop' in kwargs:
 46        r_stop = kwargs.get('r_stop')
 47        if len(r_stop) != replica:
 48            raise Exception('r_stop does not match number of replicas')
 49    else:
 50        r_stop = [None] * replica
 51
 52    print(r'Read <bar{psi}\psi> from', prefix[:-1], ',', replica, 'replica', end='')
 53
 54    print_err = 0
 55    if 'print_err' in kwargs:
 56        print_err = 1
 57        print()
 58
 59    deltas = []
 60
 61    for rep in range(replica):
 62        tmp_array = []
 63        with open(path + '/' + ls[rep], 'rb') as fp:
 64
 65            t = fp.read(4)  # number of reweighting factors
 66            if rep == 0:
 67                nrw = struct.unpack('i', t)[0]
 68                for k in range(nrw):
 69                    deltas.append([])
 70            else:
 71                if nrw != struct.unpack('i', t)[0]:
 72                    raise Exception('Error: different number of factors for replicum', rep)
 73
 74            for k in range(nrw):
 75                tmp_array.append([])
 76
 77            # This block is necessary for openQCD1.6 ms1 files
 78            nfct = []
 79            for i in range(nrw):
 80                t = fp.read(4)
 81                nfct.append(struct.unpack('i', t)[0])
 82            print('nfct: ', nfct)  # Hasenbusch factor, 1 for rat reweighting
 83
 84            nsrc = []
 85            for i in range(nrw):
 86                t = fp.read(4)
 87                nsrc.append(struct.unpack('i', t)[0])
 88
 89            # body
 90            while True:
 91                t = fp.read(4)
 92                if len(t) < 4:
 93                    break
 94                if print_err:
 95                    config_no = struct.unpack('i', t)
 96                for i in range(nrw):
 97                    tmp_nfct = 1.0
 98                    for j in range(nfct[i]):
 99                        t = fp.read(8 * nsrc[i])
100                        t = fp.read(8 * nsrc[i])
101                        tmp_rw = struct.unpack('d' * nsrc[i], t)
102                        tmp_nfct *= np.mean(np.asarray(tmp_rw))
103                        if print_err:
104                            print(config_no, i, j, np.mean(np.asarray(tmp_rw)), np.std(np.asarray(tmp_rw)))
105                            print('Sources:', np.asarray(tmp_rw))
106                            print('Partial factor:', tmp_nfct)
107                    tmp_array[i].append(tmp_nfct)
108
109            for k in range(nrw):
110                deltas[k].append(tmp_array[k][r_start[rep]:r_stop[rep]])
111
112    rep_names = []
113    for entry in ls:
114        truncated_entry = entry.split('.')[0]
115        idx = truncated_entry.index('r')
116        rep_names.append(truncated_entry[:idx] + '|' + truncated_entry[idx:])
117    print(',', nrw, r'<bar{psi}\psi> with', nsrc, 'sources')
118    result = []
119    for t in range(nrw):
120        result.append(Obs(deltas[t], rep_names))
121
122    return result