From cb993cb92e3b0e59d4a5620af9197b1a66169588 Mon Sep 17 00:00:00 2001 From: fjosw Date: Wed, 8 Feb 2023 14:53:27 +0000 Subject: [PATCH] Documentation updated --- docs/pyerrors/input/openQCD.html | 3896 +++++++++++++++--------------- 1 file changed, 1954 insertions(+), 1942 deletions(-) diff --git a/docs/pyerrors/input/openQCD.html b/docs/pyerrors/input/openQCD.html index b8974f9d..af63fa04 100644 --- a/docs/pyerrors/input/openQCD.html +++ b/docs/pyerrors/input/openQCD.html @@ -108,1159 +108,1181 @@ 12from ..correlators import Corr 13 14 - 15def read_rwms(path, prefix, version='2.0', names=None, **kwargs): - 16 """Read rwms format from given folder structure. Returns a list of length nrw - 17 - 18 Parameters - 19 ---------- - 20 path : str - 21 path that contains the data files - 22 prefix : str - 23 all files in path that start with prefix are considered as input files. - 24 May be used together postfix to consider only special file endings. - 25 Prefix is ignored, if the keyword 'files' is used. - 26 version : str - 27 version of openQCD, default 2.0 - 28 names : list - 29 list of names that is assigned to the data according according - 30 to the order in the file list. Use careful, if you do not provide file names! - 31 r_start : list - 32 list which contains the first config to be read for each replicum - 33 r_stop : list - 34 list which contains the last config to be read for each replicum - 35 r_step : int - 36 integer that defines a fixed step size between two measurements (in units of configs) - 37 If not given, r_step=1 is assumed. - 38 postfix : str - 39 postfix of the file to read, e.g. '.ms1' for openQCD-files - 40 files : list - 41 list which contains the filenames to be read. No automatic detection of - 42 files performed if given. - 43 print_err : bool - 44 Print additional information that is useful for debugging. - 45 - 46 Returns - 47 ------- - 48 rwms : Obs - 49 Reweighting factors read - 50 """ - 51 known_oqcd_versions = ['1.4', '1.6', '2.0'] - 52 if not (version in known_oqcd_versions): - 53 raise Exception('Unknown openQCD version defined!') - 54 print("Working with openQCD version " + version) - 55 if 'postfix' in kwargs: - 56 postfix = kwargs.get('postfix') - 57 else: - 58 postfix = '' - 59 ls = [] - 60 for (dirpath, dirnames, filenames) in os.walk(path): - 61 ls.extend(filenames) - 62 break - 63 - 64 if not ls: - 65 raise Exception(f"Error, directory '{path}' not found") - 66 if 'files' in kwargs: - 67 ls = kwargs.get('files') - 68 else: - 69 for exc in ls: - 70 if not fnmatch.fnmatch(exc, prefix + '*' + postfix + '.dat'): - 71 ls = list(set(ls) - set([exc])) - 72 if len(ls) > 1: - 73 ls.sort(key=lambda x: int(re.findall(r'\d+', x[len(prefix):])[0])) - 74 replica = len(ls) - 75 - 76 if 'r_start' in kwargs: - 77 r_start = kwargs.get('r_start') - 78 if len(r_start) != replica: - 79 raise Exception('r_start does not match number of replicas') - 80 r_start = [o if o else None for o in r_start] - 81 else: - 82 r_start = [None] * replica - 83 - 84 if 'r_stop' in kwargs: - 85 r_stop = kwargs.get('r_stop') - 86 if len(r_stop) != replica: - 87 raise Exception('r_stop does not match number of replicas') - 88 else: - 89 r_stop = [None] * replica - 90 - 91 if 'r_step' in kwargs: - 92 r_step = kwargs.get('r_step') - 93 else: - 94 r_step = 1 - 95 - 96 print('Read reweighting factors from', prefix[:-1], ',', - 97 replica, 'replica', end='') - 98 - 99 if names is None: - 100 rep_names = [] - 101 for entry in ls: - 102 truncated_entry = entry - 103 suffixes = [".dat", ".rwms", ".ms1"] - 104 for suffix in suffixes: - 105 if truncated_entry.endswith(suffix): - 106 truncated_entry = truncated_entry[0:-len(suffix)] - 107 idx = truncated_entry.index('r') - 108 rep_names.append(truncated_entry[:idx] + '|' + truncated_entry[idx:]) - 109 else: - 110 rep_names = names - 111 - 112 print_err = 0 - 113 if 'print_err' in kwargs: - 114 print_err = 1 - 115 print() + 15def _find_files(path, prefix, postfix, ext, known_files=[]): + 16 found = [] + 17 files = [] + 18 + 19 if postfix != "": + 20 if postfix[-1] != ".": + 21 postfix = postfix + "." + 22 if postfix[0] != ".": + 23 postfix = "." + postfix + 24 + 25 if ext[0] == ".": + 26 ext = ext[1:] + 27 + 28 pattern = prefix + "*" + postfix + ext + 29 + 30 for (dirpath, dirnames, filenames) in os.walk(path + "/"): + 31 found.extend(filenames) + 32 break + 33 + 34 if known_files != []: + 35 for kf in known_files: + 36 if kf not in found: + 37 raise FileNotFoundError("Given file " + kf + " does not exist!") + 38 + 39 return known_files + 40 + 41 if not found: + 42 raise FileNotFoundError(f"Error, directory '{path}' not found") + 43 + 44 for f in found: + 45 if fnmatch.fnmatch(f, pattern): + 46 files.append(f) + 47 + 48 if files == []: + 49 raise Exception("No files found after pattern filter!") + 50 + 51 files.sort(key=lambda x: int(re.findall(r'\d+', x[len(prefix):])[0])) + 52 return files + 53 + 54 + 55def read_rwms(path, prefix, version='2.0', names=None, **kwargs): + 56 """Read rwms format from given folder structure. Returns a list of length nrw + 57 + 58 Parameters + 59 ---------- + 60 path : str + 61 path that contains the data files + 62 prefix : str + 63 all files in path that start with prefix are considered as input files. + 64 May be used together postfix to consider only special file endings. + 65 Prefix is ignored, if the keyword 'files' is used. + 66 version : str + 67 version of openQCD, default 2.0 + 68 names : list + 69 list of names that is assigned to the data according according + 70 to the order in the file list. Use careful, if you do not provide file names! + 71 r_start : list + 72 list which contains the first config to be read for each replicum + 73 r_stop : list + 74 list which contains the last config to be read for each replicum + 75 r_step : int + 76 integer that defines a fixed step size between two measurements (in units of configs) + 77 If not given, r_step=1 is assumed. + 78 postfix : str + 79 postfix of the file to read, e.g. '.ms1' for openQCD-files + 80 files : list + 81 list which contains the filenames to be read. No automatic detection of + 82 files performed if given. + 83 print_err : bool + 84 Print additional information that is useful for debugging. + 85 + 86 Returns + 87 ------- + 88 rwms : Obs + 89 Reweighting factors read + 90 """ + 91 known_oqcd_versions = ['1.4', '1.6', '2.0'] + 92 if not (version in known_oqcd_versions): + 93 raise Exception('Unknown openQCD version defined!') + 94 print("Working with openQCD version " + version) + 95 if 'postfix' in kwargs: + 96 postfix = kwargs.get('postfix') + 97 else: + 98 postfix = '' + 99 + 100 if 'files' in kwargs: + 101 known_files = kwargs.get('files') + 102 else: + 103 known_files = [] + 104 + 105 ls = _find_files(path, prefix, postfix, 'dat', known_files=known_files) + 106 + 107 replica = len(ls) + 108 + 109 if 'r_start' in kwargs: + 110 r_start = kwargs.get('r_start') + 111 if len(r_start) != replica: + 112 raise Exception('r_start does not match number of replicas') + 113 r_start = [o if o else None for o in r_start] + 114 else: + 115 r_start = [None] * replica 116 - 117 deltas = [] - 118 - 119 configlist = [] - 120 r_start_index = [] - 121 r_stop_index = [] - 122 - 123 for rep in range(replica): - 124 tmp_array = [] - 125 with open(path + '/' + ls[rep], 'rb') as fp: - 126 - 127 t = fp.read(4) # number of reweighting factors - 128 if rep == 0: - 129 nrw = struct.unpack('i', t)[0] - 130 if version == '2.0': - 131 nrw = int(nrw / 2) - 132 for k in range(nrw): - 133 deltas.append([]) - 134 else: - 135 if ((nrw != struct.unpack('i', t)[0] and (not version == '2.0')) or (nrw != struct.unpack('i', t)[0] / 2 and version == '2.0')): - 136 raise Exception('Error: different number of reweighting factors for replicum', rep) - 137 - 138 for k in range(nrw): - 139 tmp_array.append([]) - 140 - 141 # This block is necessary for openQCD1.6 and openQCD2.0 ms1 files - 142 nfct = [] - 143 if version in ['1.6', '2.0']: - 144 for i in range(nrw): - 145 t = fp.read(4) - 146 nfct.append(struct.unpack('i', t)[0]) - 147 else: - 148 for i in range(nrw): - 149 nfct.append(1) - 150 - 151 nsrc = [] - 152 for i in range(nrw): - 153 t = fp.read(4) - 154 nsrc.append(struct.unpack('i', t)[0]) - 155 if version == '2.0': - 156 if not struct.unpack('i', fp.read(4))[0] == 0: - 157 print('something is wrong!') - 158 - 159 configlist.append([]) - 160 while True: - 161 t = fp.read(4) - 162 if len(t) < 4: - 163 break - 164 config_no = struct.unpack('i', t)[0] - 165 configlist[-1].append(config_no) - 166 for i in range(nrw): - 167 if (version == '2.0'): - 168 tmpd = _read_array_openQCD2(fp) - 169 tmpd = _read_array_openQCD2(fp) - 170 tmp_rw = tmpd['arr'] - 171 tmp_nfct = 1.0 - 172 for j in range(tmpd['n'][0]): - 173 tmp_nfct *= np.mean(np.exp(-np.asarray(tmp_rw[j]))) - 174 if print_err: - 175 print(config_no, i, j, - 176 np.mean(np.exp(-np.asarray(tmp_rw[j]))), - 177 np.std(np.exp(-np.asarray(tmp_rw[j])))) - 178 print('Sources:', - 179 np.exp(-np.asarray(tmp_rw[j]))) - 180 print('Partial factor:', tmp_nfct) - 181 elif version == '1.6' or version == '1.4': - 182 tmp_nfct = 1.0 - 183 for j in range(nfct[i]): - 184 t = fp.read(8 * nsrc[i]) - 185 t = fp.read(8 * nsrc[i]) - 186 tmp_rw = struct.unpack('d' * nsrc[i], t) - 187 tmp_nfct *= np.mean(np.exp(-np.asarray(tmp_rw))) - 188 if print_err: - 189 print(config_no, i, j, - 190 np.mean(np.exp(-np.asarray(tmp_rw))), - 191 np.std(np.exp(-np.asarray(tmp_rw)))) - 192 print('Sources:', np.exp(-np.asarray(tmp_rw))) - 193 print('Partial factor:', tmp_nfct) - 194 tmp_array[i].append(tmp_nfct) - 195 - 196 diffmeas = configlist[-1][-1] - configlist[-1][-2] - 197 configlist[-1] = [item // diffmeas for item in configlist[-1]] - 198 if configlist[-1][0] > 1 and diffmeas > 1: - 199 warnings.warn('Assume thermalization and that the first measurement belongs to the first config.') - 200 offset = configlist[-1][0] - 1 - 201 configlist[-1] = [item - offset for item in configlist[-1]] - 202 - 203 if r_start[rep] is None: - 204 r_start_index.append(0) - 205 else: - 206 try: - 207 r_start_index.append(configlist[-1].index(r_start[rep])) - 208 except ValueError: - 209 raise Exception('Config %d not in file with range [%d, %d]' % ( - 210 r_start[rep], configlist[-1][0], configlist[-1][-1])) from None - 211 - 212 if r_stop[rep] is None: - 213 r_stop_index.append(len(configlist[-1]) - 1) - 214 else: - 215 try: - 216 r_stop_index.append(configlist[-1].index(r_stop[rep])) - 217 except ValueError: - 218 raise Exception('Config %d not in file with range [%d, %d]' % ( - 219 r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None - 220 - 221 for k in range(nrw): - 222 deltas[k].append(tmp_array[k][r_start_index[rep]:r_stop_index[rep] + 1][::r_step]) - 223 - 224 if np.any([len(np.unique(np.diff(cl))) != 1 for cl in configlist]): - 225 raise Exception('Irregular spaced data in input file!', [len(np.unique(np.diff(cl))) for cl in configlist]) - 226 stepsizes = [list(np.unique(np.diff(cl)))[0] for cl in configlist] - 227 if np.any([step != 1 for step in stepsizes]): - 228 warnings.warn('Stepsize between configurations is greater than one!' + str(stepsizes), RuntimeWarning) - 229 - 230 print(',', nrw, 'reweighting factors with', nsrc, 'sources') - 231 result = [] - 232 idl = [range(configlist[rep][r_start_index[rep]], configlist[rep][r_stop_index[rep]] + 1, r_step) for rep in range(replica)] - 233 - 234 for t in range(nrw): - 235 result.append(Obs(deltas[t], rep_names, idl=idl)) - 236 return result - 237 - 238 - 239def extract_t0(path, prefix, dtr_read, xmin, spatial_extent, fit_range=5, **kwargs): - 240 """Extract t0 from given .ms.dat files. Returns t0 as Obs. - 241 - 242 It is assumed that all boundary effects have - 243 sufficiently decayed at x0=xmin. - 244 The data around the zero crossing of t^2<E> - 0.3 - 245 is fitted with a linear function - 246 from which the exact root is extracted. - 247 - 248 It is assumed that one measurement is performed for each config. - 249 If this is not the case, the resulting idl, as well as the handling - 250 of r_start, r_stop and r_step is wrong and the user has to correct - 251 this in the resulting observable. - 252 - 253 Parameters - 254 ---------- - 255 path : str - 256 Path to .ms.dat files - 257 prefix : str - 258 Ensemble prefix - 259 dtr_read : int - 260 Determines how many trajectories should be skipped - 261 when reading the ms.dat files. - 262 Corresponds to dtr_cnfg / dtr_ms in the openQCD input file. - 263 xmin : int - 264 First timeslice where the boundary - 265 effects have sufficiently decayed. - 266 spatial_extent : int - 267 spatial extent of the lattice, required for normalization. - 268 fit_range : int - 269 Number of data points left and right of the zero - 270 crossing to be included in the linear fit. (Default: 5) - 271 r_start : list - 272 list which contains the first config to be read for each replicum. - 273 r_stop : list - 274 list which contains the last config to be read for each replicum. - 275 r_step : int - 276 integer that defines a fixed step size between two measurements (in units of configs) - 277 If not given, r_step=1 is assumed. - 278 plaquette : bool - 279 If true extract the plaquette estimate of t0 instead. - 280 names : list - 281 list of names that is assigned to the data according according - 282 to the order in the file list. Use careful, if you do not provide file names! - 283 files : list - 284 list which contains the filenames to be read. No automatic detection of - 285 files performed if given. - 286 plot_fit : bool - 287 If true, the fit for the extraction of t0 is shown together with the data. - 288 assume_thermalization : bool - 289 If True: If the first record divided by the distance between two measurements is larger than - 290 1, it is assumed that this is due to thermalization and the first measurement belongs - 291 to the first config (default). - 292 If False: The config numbers are assumed to be traj_number // difference - 293 - 294 Returns - 295 ------- - 296 t0 : Obs - 297 Extracted t0 - 298 """ - 299 - 300 ls = [] - 301 for (dirpath, dirnames, filenames) in os.walk(path): - 302 ls.extend(filenames) - 303 break - 304 - 305 if not ls: - 306 raise Exception('Error, directory not found') - 307 - 308 if 'files' in kwargs: - 309 ls = kwargs.get('files') - 310 else: - 311 for exc in ls: - 312 if not fnmatch.fnmatch(exc, prefix + '*.ms.dat'): - 313 ls = list(set(ls) - set([exc])) - 314 if len(ls) > 1: - 315 ls.sort(key=lambda x: int(re.findall(r'\d+', x[len(prefix):])[0])) - 316 replica = len(ls) - 317 - 318 if 'r_start' in kwargs: - 319 r_start = kwargs.get('r_start') - 320 if len(r_start) != replica: - 321 raise Exception('r_start does not match number of replicas') - 322 r_start = [o if o else None for o in r_start] - 323 else: - 324 r_start = [None] * replica - 325 - 326 if 'r_stop' in kwargs: - 327 r_stop = kwargs.get('r_stop') - 328 if len(r_stop) != replica: - 329 raise Exception('r_stop does not match number of replicas') - 330 else: - 331 r_stop = [None] * replica + 117 if 'r_stop' in kwargs: + 118 r_stop = kwargs.get('r_stop') + 119 if len(r_stop) != replica: + 120 raise Exception('r_stop does not match number of replicas') + 121 else: + 122 r_stop = [None] * replica + 123 + 124 if 'r_step' in kwargs: + 125 r_step = kwargs.get('r_step') + 126 else: + 127 r_step = 1 + 128 + 129 print('Read reweighting factors from', prefix[:-1], ',', + 130 replica, 'replica', end='') + 131 + 132 if names is None: + 133 rep_names = [] + 134 for entry in ls: + 135 truncated_entry = entry + 136 suffixes = [".dat", ".rwms", ".ms1"] + 137 for suffix in suffixes: + 138 if truncated_entry.endswith(suffix): + 139 truncated_entry = truncated_entry[0:-len(suffix)] + 140 idx = truncated_entry.index('r') + 141 rep_names.append(truncated_entry[:idx] + '|' + truncated_entry[idx:]) + 142 else: + 143 rep_names = names + 144 + 145 print_err = 0 + 146 if 'print_err' in kwargs: + 147 print_err = 1 + 148 print() + 149 + 150 deltas = [] + 151 + 152 configlist = [] + 153 r_start_index = [] + 154 r_stop_index = [] + 155 + 156 for rep in range(replica): + 157 tmp_array = [] + 158 with open(path + '/' + ls[rep], 'rb') as fp: + 159 + 160 t = fp.read(4) # number of reweighting factors + 161 if rep == 0: + 162 nrw = struct.unpack('i', t)[0] + 163 if version == '2.0': + 164 nrw = int(nrw / 2) + 165 for k in range(nrw): + 166 deltas.append([]) + 167 else: + 168 if ((nrw != struct.unpack('i', t)[0] and (not version == '2.0')) or (nrw != struct.unpack('i', t)[0] / 2 and version == '2.0')): + 169 raise Exception('Error: different number of reweighting factors for replicum', rep) + 170 + 171 for k in range(nrw): + 172 tmp_array.append([]) + 173 + 174 # This block is necessary for openQCD1.6 and openQCD2.0 ms1 files + 175 nfct = [] + 176 if version in ['1.6', '2.0']: + 177 for i in range(nrw): + 178 t = fp.read(4) + 179 nfct.append(struct.unpack('i', t)[0]) + 180 else: + 181 for i in range(nrw): + 182 nfct.append(1) + 183 + 184 nsrc = [] + 185 for i in range(nrw): + 186 t = fp.read(4) + 187 nsrc.append(struct.unpack('i', t)[0]) + 188 if version == '2.0': + 189 if not struct.unpack('i', fp.read(4))[0] == 0: + 190 raise Exception("You are using the input for openQCD version 2.0, this is not correct.") + 191 + 192 configlist.append([]) + 193 while True: + 194 t = fp.read(4) + 195 if len(t) < 4: + 196 break + 197 config_no = struct.unpack('i', t)[0] + 198 configlist[-1].append(config_no) + 199 for i in range(nrw): + 200 if (version == '2.0'): + 201 tmpd = _read_array_openQCD2(fp) + 202 tmpd = _read_array_openQCD2(fp) + 203 tmp_rw = tmpd['arr'] + 204 tmp_nfct = 1.0 + 205 for j in range(tmpd['n'][0]): + 206 tmp_nfct *= np.mean(np.exp(-np.asarray(tmp_rw[j]))) + 207 if print_err: + 208 print(config_no, i, j, + 209 np.mean(np.exp(-np.asarray(tmp_rw[j]))), + 210 np.std(np.exp(-np.asarray(tmp_rw[j])))) + 211 print('Sources:', + 212 np.exp(-np.asarray(tmp_rw[j]))) + 213 print('Partial factor:', tmp_nfct) + 214 elif version == '1.6' or version == '1.4': + 215 tmp_nfct = 1.0 + 216 for j in range(nfct[i]): + 217 t = fp.read(8 * nsrc[i]) + 218 t = fp.read(8 * nsrc[i]) + 219 tmp_rw = struct.unpack('d' * nsrc[i], t) + 220 tmp_nfct *= np.mean(np.exp(-np.asarray(tmp_rw))) + 221 if print_err: + 222 print(config_no, i, j, + 223 np.mean(np.exp(-np.asarray(tmp_rw))), + 224 np.std(np.exp(-np.asarray(tmp_rw)))) + 225 print('Sources:', np.exp(-np.asarray(tmp_rw))) + 226 print('Partial factor:', tmp_nfct) + 227 tmp_array[i].append(tmp_nfct) + 228 + 229 diffmeas = configlist[-1][-1] - configlist[-1][-2] + 230 configlist[-1] = [item // diffmeas for item in configlist[-1]] + 231 if configlist[-1][0] > 1 and diffmeas > 1: + 232 warnings.warn('Assume thermalization and that the first measurement belongs to the first config.') + 233 offset = configlist[-1][0] - 1 + 234 configlist[-1] = [item - offset for item in configlist[-1]] + 235 + 236 if r_start[rep] is None: + 237 r_start_index.append(0) + 238 else: + 239 try: + 240 r_start_index.append(configlist[-1].index(r_start[rep])) + 241 except ValueError: + 242 raise Exception('Config %d not in file with range [%d, %d]' % ( + 243 r_start[rep], configlist[-1][0], configlist[-1][-1])) from None + 244 + 245 if r_stop[rep] is None: + 246 r_stop_index.append(len(configlist[-1]) - 1) + 247 else: + 248 try: + 249 r_stop_index.append(configlist[-1].index(r_stop[rep])) + 250 except ValueError: + 251 raise Exception('Config %d not in file with range [%d, %d]' % ( + 252 r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None + 253 + 254 for k in range(nrw): + 255 deltas[k].append(tmp_array[k][r_start_index[rep]:r_stop_index[rep] + 1][::r_step]) + 256 + 257 if np.any([len(np.unique(np.diff(cl))) != 1 for cl in configlist]): + 258 raise Exception('Irregular spaced data in input file!', [len(np.unique(np.diff(cl))) for cl in configlist]) + 259 stepsizes = [list(np.unique(np.diff(cl)))[0] for cl in configlist] + 260 if np.any([step != 1 for step in stepsizes]): + 261 warnings.warn('Stepsize between configurations is greater than one!' + str(stepsizes), RuntimeWarning) + 262 + 263 print(',', nrw, 'reweighting factors with', nsrc, 'sources') + 264 result = [] + 265 idl = [range(configlist[rep][r_start_index[rep]], configlist[rep][r_stop_index[rep]] + 1, r_step) for rep in range(replica)] + 266 + 267 for t in range(nrw): + 268 result.append(Obs(deltas[t], rep_names, idl=idl)) + 269 return result + 270 + 271 + 272def extract_t0(path, prefix, dtr_read, xmin, spatial_extent, fit_range=5, **kwargs): + 273 """Extract t0 from given .ms.dat files. Returns t0 as Obs. + 274 + 275 It is assumed that all boundary effects have + 276 sufficiently decayed at x0=xmin. + 277 The data around the zero crossing of t^2<E> - 0.3 + 278 is fitted with a linear function + 279 from which the exact root is extracted. + 280 + 281 It is assumed that one measurement is performed for each config. + 282 If this is not the case, the resulting idl, as well as the handling + 283 of r_start, r_stop and r_step is wrong and the user has to correct + 284 this in the resulting observable. + 285 + 286 Parameters + 287 ---------- + 288 path : str + 289 Path to .ms.dat files + 290 prefix : str + 291 Ensemble prefix + 292 dtr_read : int + 293 Determines how many trajectories should be skipped + 294 when reading the ms.dat files. + 295 Corresponds to dtr_cnfg / dtr_ms in the openQCD input file. + 296 xmin : int + 297 First timeslice where the boundary + 298 effects have sufficiently decayed. + 299 spatial_extent : int + 300 spatial extent of the lattice, required for normalization. + 301 fit_range : int + 302 Number of data points left and right of the zero + 303 crossing to be included in the linear fit. (Default: 5) + 304 r_start : list + 305 list which contains the first config to be read for each replicum. + 306 r_stop : list + 307 list which contains the last config to be read for each replicum. + 308 r_step : int + 309 integer that defines a fixed step size between two measurements (in units of configs) + 310 If not given, r_step=1 is assumed. + 311 plaquette : bool + 312 If true extract the plaquette estimate of t0 instead. + 313 names : list + 314 list of names that is assigned to the data according according + 315 to the order in the file list. Use careful, if you do not provide file names! + 316 files : list + 317 list which contains the filenames to be read. No automatic detection of + 318 files performed if given. + 319 plot_fit : bool + 320 If true, the fit for the extraction of t0 is shown together with the data. + 321 assume_thermalization : bool + 322 If True: If the first record divided by the distance between two measurements is larger than + 323 1, it is assumed that this is due to thermalization and the first measurement belongs + 324 to the first config (default). + 325 If False: The config numbers are assumed to be traj_number // difference + 326 + 327 Returns + 328 ------- + 329 t0 : Obs + 330 Extracted t0 + 331 """ 332 - 333 if 'r_step' in kwargs: - 334 r_step = kwargs.get('r_step') + 333 if 'files' in kwargs: + 334 known_files = kwargs.get('files') 335 else: - 336 r_step = 1 + 336 known_files = [] 337 - 338 print('Extract t0 from', prefix, ',', replica, 'replica') + 338 ls = _find_files(path, prefix, 'ms', 'dat', known_files=known_files) 339 - 340 if 'names' in kwargs: - 341 rep_names = kwargs.get('names') - 342 else: - 343 rep_names = [] - 344 for entry in ls: - 345 truncated_entry = entry.split('.')[0] - 346 idx = truncated_entry.index('r') - 347 rep_names.append(truncated_entry[:idx] + '|' + truncated_entry[idx:]) - 348 - 349 Ysum = [] - 350 - 351 configlist = [] - 352 r_start_index = [] - 353 r_stop_index = [] - 354 - 355 for rep in range(replica): + 340 replica = len(ls) + 341 + 342 if 'r_start' in kwargs: + 343 r_start = kwargs.get('r_start') + 344 if len(r_start) != replica: + 345 raise Exception('r_start does not match number of replicas') + 346 r_start = [o if o else None for o in r_start] + 347 else: + 348 r_start = [None] * replica + 349 + 350 if 'r_stop' in kwargs: + 351 r_stop = kwargs.get('r_stop') + 352 if len(r_stop) != replica: + 353 raise Exception('r_stop does not match number of replicas') + 354 else: + 355 r_stop = [None] * replica 356 - 357 with open(path + '/' + ls[rep], 'rb') as fp: - 358 t = fp.read(12) - 359 header = struct.unpack('iii', t) - 360 if rep == 0: - 361 dn = header[0] - 362 nn = header[1] - 363 tmax = header[2] - 364 elif dn != header[0] or nn != header[1] or tmax != header[2]: - 365 raise Exception('Replica parameters do not match.') - 366 - 367 t = fp.read(8) - 368 if rep == 0: - 369 eps = struct.unpack('d', t)[0] - 370 print('Step size:', eps, ', Maximal t value:', dn * (nn) * eps) - 371 elif eps != struct.unpack('d', t)[0]: - 372 raise Exception('Values for eps do not match among replica.') - 373 - 374 Ysl = [] - 375 - 376 configlist.append([]) - 377 while True: - 378 t = fp.read(4) - 379 if (len(t) < 4): - 380 break - 381 nc = struct.unpack('i', t)[0] - 382 configlist[-1].append(nc) - 383 - 384 t = fp.read(8 * tmax * (nn + 1)) - 385 if kwargs.get('plaquette'): - 386 if nc % dtr_read == 0: - 387 Ysl.append(struct.unpack('d' * tmax * (nn + 1), t)) - 388 t = fp.read(8 * tmax * (nn + 1)) - 389 if not kwargs.get('plaquette'): - 390 if nc % dtr_read == 0: - 391 Ysl.append(struct.unpack('d' * tmax * (nn + 1), t)) - 392 t = fp.read(8 * tmax * (nn + 1)) - 393 - 394 Ysum.append([]) - 395 for i, item in enumerate(Ysl): - 396 Ysum[-1].append([np.mean(item[current + xmin: - 397 current + tmax - xmin]) - 398 for current in range(0, len(item), tmax)]) + 357 if 'r_step' in kwargs: + 358 r_step = kwargs.get('r_step') + 359 else: + 360 r_step = 1 + 361 + 362 print('Extract t0 from', prefix, ',', replica, 'replica') + 363 + 364 if 'names' in kwargs: + 365 rep_names = kwargs.get('names') + 366 else: + 367 rep_names = [] + 368 for entry in ls: + 369 truncated_entry = entry.split('.')[0] + 370 idx = truncated_entry.index('r') + 371 rep_names.append(truncated_entry[:idx] + '|' + truncated_entry[idx:]) + 372 + 373 Ysum = [] + 374 + 375 configlist = [] + 376 r_start_index = [] + 377 r_stop_index = [] + 378 + 379 for rep in range(replica): + 380 + 381 with open(path + '/' + ls[rep], 'rb') as fp: + 382 t = fp.read(12) + 383 header = struct.unpack('iii', t) + 384 if rep == 0: + 385 dn = header[0] + 386 nn = header[1] + 387 tmax = header[2] + 388 elif dn != header[0] or nn != header[1] or tmax != header[2]: + 389 raise Exception('Replica parameters do not match.') + 390 + 391 t = fp.read(8) + 392 if rep == 0: + 393 eps = struct.unpack('d', t)[0] + 394 print('Step size:', eps, ', Maximal t value:', dn * (nn) * eps) + 395 elif eps != struct.unpack('d', t)[0]: + 396 raise Exception('Values for eps do not match among replica.') + 397 + 398 Ysl = [] 399 - 400 diffmeas = configlist[-1][-1] - configlist[-1][-2] - 401 configlist[-1] = [item // diffmeas for item in configlist[-1]] - 402 if kwargs.get('assume_thermalization', True) and configlist[-1][0] > 1: - 403 warnings.warn('Assume thermalization and that the first measurement belongs to the first config.') - 404 offset = configlist[-1][0] - 1 - 405 configlist[-1] = [item - offset for item in configlist[-1]] - 406 - 407 if r_start[rep] is None: - 408 r_start_index.append(0) - 409 else: - 410 try: - 411 r_start_index.append(configlist[-1].index(r_start[rep])) - 412 except ValueError: - 413 raise Exception('Config %d not in file with range [%d, %d]' % ( - 414 r_start[rep], configlist[-1][0], configlist[-1][-1])) from None - 415 - 416 if r_stop[rep] is None: - 417 r_stop_index.append(len(configlist[-1]) - 1) - 418 else: - 419 try: - 420 r_stop_index.append(configlist[-1].index(r_stop[rep])) - 421 except ValueError: - 422 raise Exception('Config %d not in file with range [%d, %d]' % ( - 423 r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None - 424 - 425 if np.any([len(np.unique(np.diff(cl))) != 1 for cl in configlist]): - 426 raise Exception('Irregular spaced data in input file!', [len(np.unique(np.diff(cl))) for cl in configlist]) - 427 stepsizes = [list(np.unique(np.diff(cl)))[0] for cl in configlist] - 428 if np.any([step != 1 for step in stepsizes]): - 429 warnings.warn('Stepsize between configurations is greater than one!' + str(stepsizes), RuntimeWarning) + 400 configlist.append([]) + 401 while True: + 402 t = fp.read(4) + 403 if (len(t) < 4): + 404 break + 405 nc = struct.unpack('i', t)[0] + 406 configlist[-1].append(nc) + 407 + 408 t = fp.read(8 * tmax * (nn + 1)) + 409 if kwargs.get('plaquette'): + 410 if nc % dtr_read == 0: + 411 Ysl.append(struct.unpack('d' * tmax * (nn + 1), t)) + 412 t = fp.read(8 * tmax * (nn + 1)) + 413 if not kwargs.get('plaquette'): + 414 if nc % dtr_read == 0: + 415 Ysl.append(struct.unpack('d' * tmax * (nn + 1), t)) + 416 t = fp.read(8 * tmax * (nn + 1)) + 417 + 418 Ysum.append([]) + 419 for i, item in enumerate(Ysl): + 420 Ysum[-1].append([np.mean(item[current + xmin: + 421 current + tmax - xmin]) + 422 for current in range(0, len(item), tmax)]) + 423 + 424 diffmeas = configlist[-1][-1] - configlist[-1][-2] + 425 configlist[-1] = [item // diffmeas for item in configlist[-1]] + 426 if kwargs.get('assume_thermalization', True) and configlist[-1][0] > 1: + 427 warnings.warn('Assume thermalization and that the first measurement belongs to the first config.') + 428 offset = configlist[-1][0] - 1 + 429 configlist[-1] = [item - offset for item in configlist[-1]] 430 - 431 idl = [range(configlist[rep][r_start_index[rep]], configlist[rep][r_stop_index[rep]] + 1, r_step) for rep in range(replica)] - 432 t2E_dict = {} - 433 for n in range(nn + 1): - 434 samples = [] - 435 for nrep, rep in enumerate(Ysum): - 436 samples.append([]) - 437 for cnfg in rep: - 438 samples[-1].append(cnfg[n]) - 439 samples[-1] = samples[-1][r_start_index[nrep]:r_stop_index[nrep] + 1][::r_step] - 440 new_obs = Obs(samples, rep_names, idl=idl) - 441 t2E_dict[n * dn * eps] = (n * dn * eps) ** 2 * new_obs / (spatial_extent ** 3) - 0.3 - 442 - 443 zero_crossing = np.argmax(np.array( - 444 [o.value for o in t2E_dict.values()]) > 0.0) - 445 - 446 x = list(t2E_dict.keys())[zero_crossing - fit_range: - 447 zero_crossing + fit_range] - 448 y = list(t2E_dict.values())[zero_crossing - fit_range: - 449 zero_crossing + fit_range] - 450 [o.gamma_method() for o in y] - 451 - 452 fit_result = fit_lin(x, y) - 453 - 454 if kwargs.get('plot_fit'): - 455 plt.figure() - 456 gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1], wspace=0.0, hspace=0.0) - 457 ax0 = plt.subplot(gs[0]) - 458 xmore = list(t2E_dict.keys())[zero_crossing - fit_range - 2: zero_crossing + fit_range + 2] - 459 ymore = list(t2E_dict.values())[zero_crossing - fit_range - 2: zero_crossing + fit_range + 2] - 460 [o.gamma_method() for o in ymore] - 461 ax0.errorbar(xmore, [yi.value for yi in ymore], yerr=[yi.dvalue for yi in ymore], fmt='x') - 462 xplot = np.linspace(np.min(x), np.max(x)) - 463 yplot = [fit_result[0] + fit_result[1] * xi for xi in xplot] - 464 [yi.gamma_method() for yi in yplot] - 465 ax0.fill_between(xplot, y1=[yi.value - yi.dvalue for yi in yplot], y2=[yi.value + yi.dvalue for yi in yplot]) - 466 retval = (-fit_result[0] / fit_result[1]) - 467 retval.gamma_method() - 468 ylim = ax0.get_ylim() - 469 ax0.fill_betweenx(ylim, x1=retval.value - retval.dvalue, x2=retval.value + retval.dvalue, color='gray', alpha=0.4) - 470 ax0.set_ylim(ylim) - 471 ax0.set_ylabel(r'$t^2 \langle E(t) \rangle - 0.3 $') - 472 xlim = ax0.get_xlim() - 473 - 474 fit_res = [fit_result[0] + fit_result[1] * xi for xi in x] - 475 residuals = (np.asarray([o.value for o in y]) - [o.value for o in fit_res]) / np.asarray([o.dvalue for o in y]) - 476 ax1 = plt.subplot(gs[1]) - 477 ax1.plot(x, residuals, 'ko', ls='none', markersize=5) - 478 ax1.tick_params(direction='out') - 479 ax1.tick_params(axis="x", bottom=True, top=True, labelbottom=True) - 480 ax1.axhline(y=0.0, ls='--', color='k') - 481 ax1.fill_between(xlim, -1.0, 1.0, alpha=0.1, facecolor='k') - 482 ax1.set_xlim(xlim) - 483 ax1.set_ylabel('Residuals') - 484 ax1.set_xlabel(r'$t/a^2$') - 485 - 486 plt.draw() - 487 return -fit_result[0] / fit_result[1] - 488 - 489 - 490def _parse_array_openQCD2(d, n, size, wa, quadrupel=False): - 491 arr = [] - 492 if d == 2: - 493 for i in range(n[0]): - 494 tmp = wa[i * n[1]:(i + 1) * n[1]] - 495 if quadrupel: - 496 tmp2 = [] - 497 for j in range(0, len(tmp), 2): - 498 tmp2.append(tmp[j]) - 499 arr.append(tmp2) - 500 else: - 501 arr.append(np.asarray(tmp)) - 502 - 503 else: - 504 raise Exception('Only two-dimensional arrays supported!') - 505 - 506 return arr - 507 - 508 - 509def _read_array_openQCD2(fp): - 510 t = fp.read(4) - 511 d = struct.unpack('i', t)[0] - 512 t = fp.read(4 * d) - 513 n = struct.unpack('%di' % (d), t) - 514 t = fp.read(4) - 515 size = struct.unpack('i', t)[0] - 516 if size == 4: - 517 types = 'i' - 518 elif size == 8: - 519 types = 'd' - 520 elif size == 16: - 521 types = 'dd' - 522 else: - 523 raise Exception("Type for size '" + str(size) + "' not known.") - 524 m = n[0] - 525 for i in range(1, d): - 526 m *= n[i] - 527 - 528 t = fp.read(m * size) - 529 tmp = struct.unpack('%d%s' % (m, types), t) - 530 - 531 arr = _parse_array_openQCD2(d, n, size, tmp, quadrupel=True) - 532 return {'d': d, 'n': n, 'size': size, 'arr': arr} - 533 - 534 - 535def read_qtop(path, prefix, c, dtr_cnfg=1, version="openQCD", **kwargs): - 536 """Read the topologial charge based on openQCD gradient flow measurements. - 537 - 538 Parameters - 539 ---------- - 540 path : str - 541 path of the measurement files - 542 prefix : str - 543 prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat. - 544 Ignored if file names are passed explicitly via keyword files. - 545 c : double - 546 Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L. - 547 dtr_cnfg : int - 548 (optional) parameter that specifies the number of measurements - 549 between two configs. - 550 If it is not set, the distance between two measurements - 551 in the file is assumed to be the distance between two configurations. - 552 steps : int - 553 (optional) Distance between two configurations in units of trajectories / - 554 cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given - 555 version : str - 556 Either openQCD or sfqcd, depending on the data. - 557 L : int - 558 spatial length of the lattice in L/a. - 559 HAS to be set if version != sfqcd, since openQCD does not provide - 560 this in the header - 561 r_start : list - 562 list which contains the first config to be read for each replicum. - 563 r_stop : list - 564 list which contains the last config to be read for each replicum. - 565 files : list - 566 specify the exact files that need to be read - 567 from path, practical if e.g. only one replicum is needed - 568 postfix : str - 569 postfix of the file to read, e.g. '.gfms.dat' for openQCD-files - 570 names : list - 571 Alternative labeling for replicas/ensembles. - 572 Has to have the appropriate length. - 573 Zeuthen_flow : bool - 574 (optional) If True, the Zeuthen flow is used for Qtop. Only possible - 575 for version=='sfqcd' If False, the Wilson flow is used. - 576 integer_charge : bool - 577 If True, the charge is rounded towards the nearest integer on each config. - 578 - 579 Returns - 580 ------- - 581 result : Obs - 582 Read topological charge - 583 """ - 584 - 585 return _read_flow_obs(path, prefix, c, dtr_cnfg=dtr_cnfg, version=version, obspos=0, **kwargs) - 586 - 587 - 588def read_gf_coupling(path, prefix, c, dtr_cnfg=1, Zeuthen_flow=True, **kwargs): - 589 """Read the gradient flow coupling based on sfqcd gradient flow measurements. See 1607.06423 for details. - 590 - 591 Note: The current implementation only works for c=0.3 and T=L. The definition of the coupling in 1607.06423 requires projection to topological charge zero which is not done within this function but has to be performed in a separate step. - 592 - 593 Parameters - 594 ---------- - 595 path : str - 596 path of the measurement files - 597 prefix : str - 598 prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat. - 599 Ignored if file names are passed explicitly via keyword files. - 600 c : double - 601 Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L. - 602 dtr_cnfg : int - 603 (optional) parameter that specifies the number of measurements - 604 between two configs. - 605 If it is not set, the distance between two measurements - 606 in the file is assumed to be the distance between two configurations. - 607 steps : int - 608 (optional) Distance between two configurations in units of trajectories / - 609 cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given - 610 r_start : list - 611 list which contains the first config to be read for each replicum. - 612 r_stop : list - 613 list which contains the last config to be read for each replicum. - 614 files : list - 615 specify the exact files that need to be read - 616 from path, practical if e.g. only one replicum is needed - 617 names : list - 618 Alternative labeling for replicas/ensembles. - 619 Has to have the appropriate length. - 620 postfix : str - 621 postfix of the file to read, e.g. '.gfms.dat' for openQCD-files - 622 Zeuthen_flow : bool - 623 (optional) If True, the Zeuthen flow is used for the coupling. If False, the Wilson flow is used. - 624 """ - 625 - 626 if c != 0.3: - 627 raise Exception("The required lattice norm is only implemented for c=0.3 at the moment.") - 628 - 629 plaq = _read_flow_obs(path, prefix, c, dtr_cnfg=dtr_cnfg, version="sfqcd", obspos=6, sum_t=False, Zeuthen_flow=Zeuthen_flow, integer_charge=False, **kwargs) - 630 C2x1 = _read_flow_obs(path, prefix, c, dtr_cnfg=dtr_cnfg, version="sfqcd", obspos=7, sum_t=False, Zeuthen_flow=Zeuthen_flow, integer_charge=False, **kwargs) - 631 L = plaq.tag["L"] - 632 T = plaq.tag["T"] - 633 - 634 if T != L: - 635 raise Exception("The required lattice norm is only implemented for T=L at the moment.") - 636 - 637 if Zeuthen_flow is not True: - 638 raise Exception("The required lattice norm is only implemented for the Zeuthen flow at the moment.") - 639 - 640 t = (c * L) ** 2 / 8 - 641 - 642 normdict = {4: 0.012341170468270, - 643 6: 0.010162691462430, - 644 8: 0.009031614807931, - 645 10: 0.008744966371393, - 646 12: 0.008650917856809, - 647 14: 8.611154391267955E-03, - 648 16: 0.008591758449508, - 649 20: 0.008575359627103, - 650 24: 0.008569387847540, - 651 28: 8.566803713382559E-03, - 652 32: 0.008565541650006, - 653 40: 8.564480684962046E-03, - 654 48: 8.564098025073460E-03, - 655 64: 8.563853943383087E-03} - 656 - 657 return t * t * (5 / 3 * plaq - 1 / 12 * C2x1) / normdict[L] - 658 - 659 - 660def _read_flow_obs(path, prefix, c, dtr_cnfg=1, version="openQCD", obspos=0, sum_t=True, **kwargs): - 661 """Read a flow observable based on openQCD gradient flow measurements. - 662 - 663 Parameters - 664 ---------- - 665 path : str - 666 path of the measurement files - 667 prefix : str - 668 prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat. - 669 Ignored if file names are passed explicitly via keyword files. - 670 c : double - 671 Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L. - 672 dtr_cnfg : int - 673 (optional) parameter that specifies the number of measurements - 674 between two configs. - 675 If it is not set, the distance between two measurements - 676 in the file is assumed to be the distance between two configurations. - 677 steps : int - 678 (optional) Distance between two configurations in units of trajectories / - 679 cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given - 680 version : str - 681 Either openQCD or sfqcd, depending on the data. - 682 obspos : int - 683 position of the obeservable in the measurement file. Only relevant for sfqcd files. - 684 sum_t : bool - 685 If true sum over all timeslices, if false only take the value at T/2. - 686 L : int - 687 spatial length of the lattice in L/a. - 688 HAS to be set if version != sfqcd, since openQCD does not provide - 689 this in the header - 690 r_start : list - 691 list which contains the first config to be read for each replicum. - 692 r_stop : list - 693 list which contains the last config to be read for each replicum. - 694 files : list - 695 specify the exact files that need to be read - 696 from path, practical if e.g. only one replicum is needed - 697 names : list - 698 Alternative labeling for replicas/ensembles. - 699 Has to have the appropriate length. - 700 postfix : str - 701 postfix of the file to read, e.g. '.gfms.dat' for openQCD-files - 702 Zeuthen_flow : bool - 703 (optional) If True, the Zeuthen flow is used for Qtop. Only possible - 704 for version=='sfqcd' If False, the Wilson flow is used. - 705 integer_charge : bool - 706 If True, the charge is rounded towards the nearest integer on each config. - 707 - 708 Returns - 709 ------- - 710 result : Obs - 711 flow observable specified - 712 """ - 713 known_versions = ["openQCD", "sfqcd"] - 714 - 715 if version not in known_versions: - 716 raise Exception("Unknown openQCD version.") - 717 if "steps" in kwargs: - 718 steps = kwargs.get("steps") - 719 if version == "sfqcd": - 720 if "L" in kwargs: - 721 supposed_L = kwargs.get("L") - 722 else: - 723 supposed_L = None - 724 postfix = ".gfms.dat" - 725 else: - 726 if "L" not in kwargs: - 727 raise Exception("This version of openQCD needs you to provide the spatial length of the lattice as parameter 'L'.") - 728 else: - 729 L = kwargs.get("L") - 730 postfix = ".ms.dat" + 431 if r_start[rep] is None: + 432 r_start_index.append(0) + 433 else: + 434 try: + 435 r_start_index.append(configlist[-1].index(r_start[rep])) + 436 except ValueError: + 437 raise Exception('Config %d not in file with range [%d, %d]' % ( + 438 r_start[rep], configlist[-1][0], configlist[-1][-1])) from None + 439 + 440 if r_stop[rep] is None: + 441 r_stop_index.append(len(configlist[-1]) - 1) + 442 else: + 443 try: + 444 r_stop_index.append(configlist[-1].index(r_stop[rep])) + 445 except ValueError: + 446 raise Exception('Config %d not in file with range [%d, %d]' % ( + 447 r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None + 448 + 449 if np.any([len(np.unique(np.diff(cl))) != 1 for cl in configlist]): + 450 raise Exception('Irregular spaced data in input file!', [len(np.unique(np.diff(cl))) for cl in configlist]) + 451 stepsizes = [list(np.unique(np.diff(cl)))[0] for cl in configlist] + 452 if np.any([step != 1 for step in stepsizes]): + 453 warnings.warn('Stepsize between configurations is greater than one!' + str(stepsizes), RuntimeWarning) + 454 + 455 idl = [range(configlist[rep][r_start_index[rep]], configlist[rep][r_stop_index[rep]] + 1, r_step) for rep in range(replica)] + 456 t2E_dict = {} + 457 for n in range(nn + 1): + 458 samples = [] + 459 for nrep, rep in enumerate(Ysum): + 460 samples.append([]) + 461 for cnfg in rep: + 462 samples[-1].append(cnfg[n]) + 463 samples[-1] = samples[-1][r_start_index[nrep]:r_stop_index[nrep] + 1][::r_step] + 464 new_obs = Obs(samples, rep_names, idl=idl) + 465 t2E_dict[n * dn * eps] = (n * dn * eps) ** 2 * new_obs / (spatial_extent ** 3) - 0.3 + 466 + 467 zero_crossing = np.argmax(np.array( + 468 [o.value for o in t2E_dict.values()]) > 0.0) + 469 + 470 x = list(t2E_dict.keys())[zero_crossing - fit_range: + 471 zero_crossing + fit_range] + 472 y = list(t2E_dict.values())[zero_crossing - fit_range: + 473 zero_crossing + fit_range] + 474 [o.gamma_method() for o in y] + 475 + 476 fit_result = fit_lin(x, y) + 477 + 478 if kwargs.get('plot_fit'): + 479 plt.figure() + 480 gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1], wspace=0.0, hspace=0.0) + 481 ax0 = plt.subplot(gs[0]) + 482 xmore = list(t2E_dict.keys())[zero_crossing - fit_range - 2: zero_crossing + fit_range + 2] + 483 ymore = list(t2E_dict.values())[zero_crossing - fit_range - 2: zero_crossing + fit_range + 2] + 484 [o.gamma_method() for o in ymore] + 485 ax0.errorbar(xmore, [yi.value for yi in ymore], yerr=[yi.dvalue for yi in ymore], fmt='x') + 486 xplot = np.linspace(np.min(x), np.max(x)) + 487 yplot = [fit_result[0] + fit_result[1] * xi for xi in xplot] + 488 [yi.gamma_method() for yi in yplot] + 489 ax0.fill_between(xplot, y1=[yi.value - yi.dvalue for yi in yplot], y2=[yi.value + yi.dvalue for yi in yplot]) + 490 retval = (-fit_result[0] / fit_result[1]) + 491 retval.gamma_method() + 492 ylim = ax0.get_ylim() + 493 ax0.fill_betweenx(ylim, x1=retval.value - retval.dvalue, x2=retval.value + retval.dvalue, color='gray', alpha=0.4) + 494 ax0.set_ylim(ylim) + 495 ax0.set_ylabel(r'$t^2 \langle E(t) \rangle - 0.3 $') + 496 xlim = ax0.get_xlim() + 497 + 498 fit_res = [fit_result[0] + fit_result[1] * xi for xi in x] + 499 residuals = (np.asarray([o.value for o in y]) - [o.value for o in fit_res]) / np.asarray([o.dvalue for o in y]) + 500 ax1 = plt.subplot(gs[1]) + 501 ax1.plot(x, residuals, 'ko', ls='none', markersize=5) + 502 ax1.tick_params(direction='out') + 503 ax1.tick_params(axis="x", bottom=True, top=True, labelbottom=True) + 504 ax1.axhline(y=0.0, ls='--', color='k') + 505 ax1.fill_between(xlim, -1.0, 1.0, alpha=0.1, facecolor='k') + 506 ax1.set_xlim(xlim) + 507 ax1.set_ylabel('Residuals') + 508 ax1.set_xlabel(r'$t/a^2$') + 509 + 510 plt.draw() + 511 return -fit_result[0] / fit_result[1] + 512 + 513 + 514def _parse_array_openQCD2(d, n, size, wa, quadrupel=False): + 515 arr = [] + 516 if d == 2: + 517 for i in range(n[0]): + 518 tmp = wa[i * n[1]:(i + 1) * n[1]] + 519 if quadrupel: + 520 tmp2 = [] + 521 for j in range(0, len(tmp), 2): + 522 tmp2.append(tmp[j]) + 523 arr.append(tmp2) + 524 else: + 525 arr.append(np.asarray(tmp)) + 526 + 527 else: + 528 raise Exception('Only two-dimensional arrays supported!') + 529 + 530 return arr + 531 + 532 + 533def _read_array_openQCD2(fp): + 534 t = fp.read(4) + 535 d = struct.unpack('i', t)[0] + 536 t = fp.read(4 * d) + 537 n = struct.unpack('%di' % (d), t) + 538 t = fp.read(4) + 539 size = struct.unpack('i', t)[0] + 540 if size == 4: + 541 types = 'i' + 542 elif size == 8: + 543 types = 'd' + 544 elif size == 16: + 545 types = 'dd' + 546 else: + 547 raise Exception("Type for size '" + str(size) + "' not known.") + 548 m = n[0] + 549 for i in range(1, d): + 550 m *= n[i] + 551 + 552 t = fp.read(m * size) + 553 tmp = struct.unpack('%d%s' % (m, types), t) + 554 + 555 arr = _parse_array_openQCD2(d, n, size, tmp, quadrupel=True) + 556 return {'d': d, 'n': n, 'size': size, 'arr': arr} + 557 + 558 + 559def read_qtop(path, prefix, c, dtr_cnfg=1, version="openQCD", **kwargs): + 560 """Read the topologial charge based on openQCD gradient flow measurements. + 561 + 562 Parameters + 563 ---------- + 564 path : str + 565 path of the measurement files + 566 prefix : str + 567 prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat. + 568 Ignored if file names are passed explicitly via keyword files. + 569 c : double + 570 Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L. + 571 dtr_cnfg : int + 572 (optional) parameter that specifies the number of measurements + 573 between two configs. + 574 If it is not set, the distance between two measurements + 575 in the file is assumed to be the distance between two configurations. + 576 steps : int + 577 (optional) Distance between two configurations in units of trajectories / + 578 cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given + 579 version : str + 580 Either openQCD or sfqcd, depending on the data. + 581 L : int + 582 spatial length of the lattice in L/a. + 583 HAS to be set if version != sfqcd, since openQCD does not provide + 584 this in the header + 585 r_start : list + 586 list which contains the first config to be read for each replicum. + 587 r_stop : list + 588 list which contains the last config to be read for each replicum. + 589 files : list + 590 specify the exact files that need to be read + 591 from path, practical if e.g. only one replicum is needed + 592 postfix : str + 593 postfix of the file to read, e.g. '.gfms.dat' for openQCD-files + 594 names : list + 595 Alternative labeling for replicas/ensembles. + 596 Has to have the appropriate length. + 597 Zeuthen_flow : bool + 598 (optional) If True, the Zeuthen flow is used for Qtop. Only possible + 599 for version=='sfqcd' If False, the Wilson flow is used. + 600 integer_charge : bool + 601 If True, the charge is rounded towards the nearest integer on each config. + 602 + 603 Returns + 604 ------- + 605 result : Obs + 606 Read topological charge + 607 """ + 608 + 609 return _read_flow_obs(path, prefix, c, dtr_cnfg=dtr_cnfg, version=version, obspos=0, **kwargs) + 610 + 611 + 612def read_gf_coupling(path, prefix, c, dtr_cnfg=1, Zeuthen_flow=True, **kwargs): + 613 """Read the gradient flow coupling based on sfqcd gradient flow measurements. See 1607.06423 for details. + 614 + 615 Note: The current implementation only works for c=0.3 and T=L. The definition of the coupling in 1607.06423 requires projection to topological charge zero which is not done within this function but has to be performed in a separate step. + 616 + 617 Parameters + 618 ---------- + 619 path : str + 620 path of the measurement files + 621 prefix : str + 622 prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat. + 623 Ignored if file names are passed explicitly via keyword files. + 624 c : double + 625 Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L. + 626 dtr_cnfg : int + 627 (optional) parameter that specifies the number of measurements + 628 between two configs. + 629 If it is not set, the distance between two measurements + 630 in the file is assumed to be the distance between two configurations. + 631 steps : int + 632 (optional) Distance between two configurations in units of trajectories / + 633 cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given + 634 r_start : list + 635 list which contains the first config to be read for each replicum. + 636 r_stop : list + 637 list which contains the last config to be read for each replicum. + 638 files : list + 639 specify the exact files that need to be read + 640 from path, practical if e.g. only one replicum is needed + 641 names : list + 642 Alternative labeling for replicas/ensembles. + 643 Has to have the appropriate length. + 644 postfix : str + 645 postfix of the file to read, e.g. '.gfms.dat' for openQCD-files + 646 Zeuthen_flow : bool + 647 (optional) If True, the Zeuthen flow is used for the coupling. If False, the Wilson flow is used. + 648 """ + 649 + 650 if c != 0.3: + 651 raise Exception("The required lattice norm is only implemented for c=0.3 at the moment.") + 652 + 653 plaq = _read_flow_obs(path, prefix, c, dtr_cnfg=dtr_cnfg, version="sfqcd", obspos=6, sum_t=False, Zeuthen_flow=Zeuthen_flow, integer_charge=False, **kwargs) + 654 C2x1 = _read_flow_obs(path, prefix, c, dtr_cnfg=dtr_cnfg, version="sfqcd", obspos=7, sum_t=False, Zeuthen_flow=Zeuthen_flow, integer_charge=False, **kwargs) + 655 L = plaq.tag["L"] + 656 T = plaq.tag["T"] + 657 + 658 if T != L: + 659 raise Exception("The required lattice norm is only implemented for T=L at the moment.") + 660 + 661 if Zeuthen_flow is not True: + 662 raise Exception("The required lattice norm is only implemented for the Zeuthen flow at the moment.") + 663 + 664 t = (c * L) ** 2 / 8 + 665 + 666 normdict = {4: 0.012341170468270, + 667 6: 0.010162691462430, + 668 8: 0.009031614807931, + 669 10: 0.008744966371393, + 670 12: 0.008650917856809, + 671 14: 8.611154391267955E-03, + 672 16: 0.008591758449508, + 673 20: 0.008575359627103, + 674 24: 0.008569387847540, + 675 28: 8.566803713382559E-03, + 676 32: 0.008565541650006, + 677 40: 8.564480684962046E-03, + 678 48: 8.564098025073460E-03, + 679 64: 8.563853943383087E-03} + 680 + 681 return t * t * (5 / 3 * plaq - 1 / 12 * C2x1) / normdict[L] + 682 + 683 + 684def _read_flow_obs(path, prefix, c, dtr_cnfg=1, version="openQCD", obspos=0, sum_t=True, **kwargs): + 685 """Read a flow observable based on openQCD gradient flow measurements. + 686 + 687 Parameters + 688 ---------- + 689 path : str + 690 path of the measurement files + 691 prefix : str + 692 prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat. + 693 Ignored if file names are passed explicitly via keyword files. + 694 c : double + 695 Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L. + 696 dtr_cnfg : int + 697 (optional) parameter that specifies the number of measurements + 698 between two configs. + 699 If it is not set, the distance between two measurements + 700 in the file is assumed to be the distance between two configurations. + 701 steps : int + 702 (optional) Distance between two configurations in units of trajectories / + 703 cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given + 704 version : str + 705 Either openQCD or sfqcd, depending on the data. + 706 obspos : int + 707 position of the obeservable in the measurement file. Only relevant for sfqcd files. + 708 sum_t : bool + 709 If true sum over all timeslices, if false only take the value at T/2. + 710 L : int + 711 spatial length of the lattice in L/a. + 712 HAS to be set if version != sfqcd, since openQCD does not provide + 713 this in the header + 714 r_start : list + 715 list which contains the first config to be read for each replicum. + 716 r_stop : list + 717 list which contains the last config to be read for each replicum. + 718 files : list + 719 specify the exact files that need to be read + 720 from path, practical if e.g. only one replicum is needed + 721 names : list + 722 Alternative labeling for replicas/ensembles. + 723 Has to have the appropriate length. + 724 postfix : str + 725 postfix of the file to read, e.g. '.gfms.dat' for openQCD-files + 726 Zeuthen_flow : bool + 727 (optional) If True, the Zeuthen flow is used for Qtop. Only possible + 728 for version=='sfqcd' If False, the Wilson flow is used. + 729 integer_charge : bool + 730 If True, the charge is rounded towards the nearest integer on each config. 731 - 732 if "postfix" in kwargs: - 733 postfix = kwargs.get("postfix") - 734 - 735 if "files" in kwargs: - 736 files = kwargs.get("files") - 737 postfix = '' - 738 else: - 739 found = [] - 740 files = [] - 741 for (dirpath, dirnames, filenames) in os.walk(path + "/"): - 742 found.extend(filenames) - 743 break - 744 for f in found: - 745 if fnmatch.fnmatch(f, prefix + "*" + postfix): - 746 files.append(f) - 747 - 748 files = sorted(files) - 749 - 750 if 'r_start' in kwargs: - 751 r_start = kwargs.get('r_start') - 752 if len(r_start) != len(files): - 753 raise Exception('r_start does not match number of replicas') - 754 r_start = [o if o else None for o in r_start] - 755 else: - 756 r_start = [None] * len(files) - 757 - 758 if 'r_stop' in kwargs: - 759 r_stop = kwargs.get('r_stop') - 760 if len(r_stop) != len(files): - 761 raise Exception('r_stop does not match number of replicas') - 762 else: - 763 r_stop = [None] * len(files) - 764 rep_names = [] + 732 Returns + 733 ------- + 734 result : Obs + 735 flow observable specified + 736 """ + 737 known_versions = ["openQCD", "sfqcd"] + 738 + 739 if version not in known_versions: + 740 raise Exception("Unknown openQCD version.") + 741 if "steps" in kwargs: + 742 steps = kwargs.get("steps") + 743 if version == "sfqcd": + 744 if "L" in kwargs: + 745 supposed_L = kwargs.get("L") + 746 else: + 747 supposed_L = None + 748 postfix = "gfms" + 749 else: + 750 if "L" not in kwargs: + 751 raise Exception("This version of openQCD needs you to provide the spatial length of the lattice as parameter 'L'.") + 752 else: + 753 L = kwargs.get("L") + 754 postfix = "ms" + 755 + 756 if "postfix" in kwargs: + 757 postfix = kwargs.get("postfix") + 758 + 759 if "files" in kwargs: + 760 known_files = kwargs.get("files") + 761 else: + 762 known_files = [] + 763 + 764 files = _find_files(path, prefix, postfix, "dat", known_files=known_files) 765 - 766 zeuthen = kwargs.get('Zeuthen_flow', False) - 767 if zeuthen and version not in ['sfqcd']: - 768 raise Exception('Zeuthen flow can only be used for version==sfqcd') - 769 - 770 r_start_index = [] - 771 r_stop_index = [] - 772 deltas = [] - 773 configlist = [] - 774 if not zeuthen: - 775 obspos += 8 - 776 for rep, file in enumerate(files): - 777 with open(path + "/" + file, "rb") as fp: - 778 - 779 Q = [] - 780 traj_list = [] - 781 if version in ['sfqcd']: - 782 t = fp.read(12) - 783 header = struct.unpack('<iii', t) - 784 zthfl = header[0] # Zeuthen flow -> if it's equal to 2 it means that the Zeuthen flow is also 'measured' (apart from the Wilson flow) - 785 ncs = header[1] # number of different values for c in t_flow=1/8 c² L² -> measurements done for ncs c's - 786 tmax = header[2] # lattice T/a - 787 - 788 t = fp.read(12) - 789 Ls = struct.unpack('<iii', t) - 790 if (Ls[0] == Ls[1] and Ls[1] == Ls[2]): - 791 L = Ls[0] - 792 if not (supposed_L == L) and supposed_L: - 793 raise Exception("It seems the length given in the header and by you contradict each other") - 794 else: - 795 raise Exception("Found more than one spatial length in header!") - 796 - 797 t = fp.read(16) - 798 header2 = struct.unpack('<dd', t) - 799 tol = header2[0] - 800 cmax = header2[1] # highest value of c used - 801 - 802 if c > cmax: - 803 raise Exception('Flow has been determined between c=0 and c=%lf with tolerance %lf' % (cmax, tol)) - 804 - 805 if (zthfl == 2): - 806 nfl = 2 # number of flows - 807 else: - 808 nfl = 1 - 809 iobs = 8 * nfl # number of flow observables calculated - 810 - 811 while True: - 812 t = fp.read(4) - 813 if (len(t) < 4): - 814 break - 815 traj_list.append(struct.unpack('i', t)[0]) # trajectory number when measurement was done - 816 - 817 for j in range(ncs + 1): - 818 for i in range(iobs): - 819 t = fp.read(8 * tmax) - 820 if (i == obspos): # determines the flow observable -> i=0 <-> Zeuthen flow - 821 Q.append(struct.unpack('d' * tmax, t)) - 822 - 823 else: - 824 t = fp.read(12) - 825 header = struct.unpack('<iii', t) - 826 # step size in integration steps "dnms" - 827 dn = header[0] - 828 # number of measurements, so "ntot"/dn - 829 nn = header[1] - 830 # lattice T/a - 831 tmax = header[2] + 766 if 'r_start' in kwargs: + 767 r_start = kwargs.get('r_start') + 768 if len(r_start) != len(files): + 769 raise Exception('r_start does not match number of replicas') + 770 r_start = [o if o else None for o in r_start] + 771 else: + 772 r_start = [None] * len(files) + 773 + 774 if 'r_stop' in kwargs: + 775 r_stop = kwargs.get('r_stop') + 776 if len(r_stop) != len(files): + 777 raise Exception('r_stop does not match number of replicas') + 778 else: + 779 r_stop = [None] * len(files) + 780 rep_names = [] + 781 + 782 zeuthen = kwargs.get('Zeuthen_flow', False) + 783 if zeuthen and version not in ['sfqcd']: + 784 raise Exception('Zeuthen flow can only be used for version==sfqcd') + 785 + 786 r_start_index = [] + 787 r_stop_index = [] + 788 deltas = [] + 789 configlist = [] + 790 if not zeuthen: + 791 obspos += 8 + 792 for rep, file in enumerate(files): + 793 with open(path + "/" + file, "rb") as fp: + 794 + 795 Q = [] + 796 traj_list = [] + 797 if version in ['sfqcd']: + 798 t = fp.read(12) + 799 header = struct.unpack('<iii', t) + 800 zthfl = header[0] # Zeuthen flow -> if it's equal to 2 it means that the Zeuthen flow is also 'measured' (apart from the Wilson flow) + 801 ncs = header[1] # number of different values for c in t_flow=1/8 c² L² -> measurements done for ncs c's + 802 tmax = header[2] # lattice T/a + 803 + 804 t = fp.read(12) + 805 Ls = struct.unpack('<iii', t) + 806 if (Ls[0] == Ls[1] and Ls[1] == Ls[2]): + 807 L = Ls[0] + 808 if not (supposed_L == L) and supposed_L: + 809 raise Exception("It seems the length given in the header and by you contradict each other") + 810 else: + 811 raise Exception("Found more than one spatial length in header!") + 812 + 813 t = fp.read(16) + 814 header2 = struct.unpack('<dd', t) + 815 tol = header2[0] + 816 cmax = header2[1] # highest value of c used + 817 + 818 if c > cmax: + 819 raise Exception('Flow has been determined between c=0 and c=%lf with tolerance %lf' % (cmax, tol)) + 820 + 821 if (zthfl == 2): + 822 nfl = 2 # number of flows + 823 else: + 824 nfl = 1 + 825 iobs = 8 * nfl # number of flow observables calculated + 826 + 827 while True: + 828 t = fp.read(4) + 829 if (len(t) < 4): + 830 break + 831 traj_list.append(struct.unpack('i', t)[0]) # trajectory number when measurement was done 832 - 833 t = fp.read(8) - 834 eps = struct.unpack('d', t)[0] - 835 - 836 while True: - 837 t = fp.read(4) - 838 if (len(t) < 4): - 839 break - 840 traj_list.append(struct.unpack('i', t)[0]) - 841 # Wsl - 842 t = fp.read(8 * tmax * (nn + 1)) - 843 # Ysl - 844 t = fp.read(8 * tmax * (nn + 1)) - 845 # Qsl, which is asked for in this method - 846 t = fp.read(8 * tmax * (nn + 1)) - 847 # unpack the array of Qtops, - 848 # on each timeslice t=0,...,tmax-1 and the - 849 # measurement number in = 0...nn (see README.qcd1) - 850 tmpd = struct.unpack('d' * tmax * (nn + 1), t) - 851 Q.append(tmpd) - 852 - 853 if len(np.unique(np.diff(traj_list))) != 1: - 854 raise Exception("Irregularities in stepsize found") - 855 else: - 856 if 'steps' in kwargs: - 857 if steps != traj_list[1] - traj_list[0]: - 858 raise Exception("steps and the found stepsize are not the same") - 859 else: - 860 steps = traj_list[1] - traj_list[0] - 861 - 862 configlist.append([tr // steps // dtr_cnfg for tr in traj_list]) - 863 if configlist[-1][0] > 1: - 864 offset = configlist[-1][0] - 1 - 865 warnings.warn('Assume thermalization and that the first measurement belongs to the first config. Offset = %d configs (%d trajectories / cycles)' % ( - 866 offset, offset * steps)) - 867 configlist[-1] = [item - offset for item in configlist[-1]] + 833 for j in range(ncs + 1): + 834 for i in range(iobs): + 835 t = fp.read(8 * tmax) + 836 if (i == obspos): # determines the flow observable -> i=0 <-> Zeuthen flow + 837 Q.append(struct.unpack('d' * tmax, t)) + 838 + 839 else: + 840 t = fp.read(12) + 841 header = struct.unpack('<iii', t) + 842 # step size in integration steps "dnms" + 843 dn = header[0] + 844 # number of measurements, so "ntot"/dn + 845 nn = header[1] + 846 # lattice T/a + 847 tmax = header[2] + 848 + 849 t = fp.read(8) + 850 eps = struct.unpack('d', t)[0] + 851 + 852 while True: + 853 t = fp.read(4) + 854 if (len(t) < 4): + 855 break + 856 traj_list.append(struct.unpack('i', t)[0]) + 857 # Wsl + 858 t = fp.read(8 * tmax * (nn + 1)) + 859 # Ysl + 860 t = fp.read(8 * tmax * (nn + 1)) + 861 # Qsl, which is asked for in this method + 862 t = fp.read(8 * tmax * (nn + 1)) + 863 # unpack the array of Qtops, + 864 # on each timeslice t=0,...,tmax-1 and the + 865 # measurement number in = 0...nn (see README.qcd1) + 866 tmpd = struct.unpack('d' * tmax * (nn + 1), t) + 867 Q.append(tmpd) 868 - 869 if r_start[rep] is None: - 870 r_start_index.append(0) + 869 if len(np.unique(np.diff(traj_list))) != 1: + 870 raise Exception("Irregularities in stepsize found") 871 else: - 872 try: - 873 r_start_index.append(configlist[-1].index(r_start[rep])) - 874 except ValueError: - 875 raise Exception('Config %d not in file with range [%d, %d]' % ( - 876 r_start[rep], configlist[-1][0], configlist[-1][-1])) from None + 872 if 'steps' in kwargs: + 873 if steps != traj_list[1] - traj_list[0]: + 874 raise Exception("steps and the found stepsize are not the same") + 875 else: + 876 steps = traj_list[1] - traj_list[0] 877 - 878 if r_stop[rep] is None: - 879 r_stop_index.append(len(configlist[-1]) - 1) - 880 else: - 881 try: - 882 r_stop_index.append(configlist[-1].index(r_stop[rep])) - 883 except ValueError: - 884 raise Exception('Config %d not in file with range [%d, %d]' % ( - 885 r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None - 886 - 887 if version in ['sfqcd']: - 888 cstepsize = cmax / ncs - 889 index_aim = round(c / cstepsize) - 890 else: - 891 t_aim = (c * L) ** 2 / 8 - 892 index_aim = round(t_aim / eps / dn) + 878 configlist.append([tr // steps // dtr_cnfg for tr in traj_list]) + 879 if configlist[-1][0] > 1: + 880 offset = configlist[-1][0] - 1 + 881 warnings.warn('Assume thermalization and that the first measurement belongs to the first config. Offset = %d configs (%d trajectories / cycles)' % ( + 882 offset, offset * steps)) + 883 configlist[-1] = [item - offset for item in configlist[-1]] + 884 + 885 if r_start[rep] is None: + 886 r_start_index.append(0) + 887 else: + 888 try: + 889 r_start_index.append(configlist[-1].index(r_start[rep])) + 890 except ValueError: + 891 raise Exception('Config %d not in file with range [%d, %d]' % ( + 892 r_start[rep], configlist[-1][0], configlist[-1][-1])) from None 893 - 894 Q_sum = [] - 895 for i, item in enumerate(Q): - 896 if sum_t is True: - 897 Q_sum.append([sum(item[current:current + tmax]) - 898 for current in range(0, len(item), tmax)]) - 899 else: - 900 Q_sum.append([item[int(tmax / 2)]]) - 901 Q_top = [] - 902 if version in ['sfqcd']: - 903 for i in range(len(Q_sum) // (ncs + 1)): - 904 Q_top.append(Q_sum[i * (ncs + 1) + index_aim][0]) - 905 else: - 906 for i in range(len(Q) // dtr_cnfg): - 907 Q_top.append(Q_sum[dtr_cnfg * i][index_aim]) - 908 if len(Q_top) != len(traj_list) // dtr_cnfg: - 909 raise Exception("qtops and traj_list dont have the same length") - 910 - 911 if kwargs.get('integer_charge', False): - 912 Q_top = [round(q) for q in Q_top] - 913 - 914 truncated_file = file[:-len(postfix)] - 915 - 916 if "names" not in kwargs: - 917 try: - 918 idx = truncated_file.index('r') - 919 except Exception: - 920 if "names" not in kwargs: - 921 raise Exception("Automatic recognition of replicum failed, please enter the key word 'names'.") - 922 ens_name = truncated_file[:idx] - 923 rep_names.append(ens_name + '|' + truncated_file[idx:]) - 924 else: - 925 names = kwargs.get("names") - 926 rep_names = names - 927 deltas.append(Q_top) - 928 - 929 idl = [range(int(configlist[rep][r_start_index[rep]]), int(configlist[rep][r_stop_index[rep]]) + 1, 1) for rep in range(len(deltas))] - 930 deltas = [deltas[nrep][r_start_index[nrep]:r_stop_index[nrep] + 1] for nrep in range(len(deltas))] - 931 result = Obs(deltas, rep_names, idl=idl) - 932 result.tag = {"T": tmax - 1, - 933 "L": L} - 934 return result - 935 - 936 - 937def qtop_projection(qtop, target=0): - 938 """Returns the projection to the topological charge sector defined by target. - 939 - 940 Parameters - 941 ---------- - 942 path : Obs - 943 Topological charge. - 944 target : int - 945 Specifies the topological sector to be reweighted to (default 0) - 946 - 947 Returns - 948 ------- - 949 reto : Obs - 950 projection to the topological charge sector defined by target - 951 """ - 952 if qtop.reweighted: - 953 raise Exception('You can not use a reweighted observable for reweighting!') - 954 - 955 proj_qtop = [] - 956 for n in qtop.deltas: - 957 proj_qtop.append(np.array([1 if round(qtop.r_values[n] + q) == target else 0 for q in qtop.deltas[n]])) - 958 - 959 reto = Obs(proj_qtop, qtop.names, idl=[qtop.idl[name] for name in qtop.names]) - 960 return reto - 961 + 894 if r_stop[rep] is None: + 895 r_stop_index.append(len(configlist[-1]) - 1) + 896 else: + 897 try: + 898 r_stop_index.append(configlist[-1].index(r_stop[rep])) + 899 except ValueError: + 900 raise Exception('Config %d not in file with range [%d, %d]' % ( + 901 r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None + 902 + 903 if version in ['sfqcd']: + 904 cstepsize = cmax / ncs + 905 index_aim = round(c / cstepsize) + 906 else: + 907 t_aim = (c * L) ** 2 / 8 + 908 index_aim = round(t_aim / eps / dn) + 909 + 910 Q_sum = [] + 911 for i, item in enumerate(Q): + 912 if sum_t is True: + 913 Q_sum.append([sum(item[current:current + tmax]) + 914 for current in range(0, len(item), tmax)]) + 915 else: + 916 Q_sum.append([item[int(tmax / 2)]]) + 917 Q_top = [] + 918 if version in ['sfqcd']: + 919 for i in range(len(Q_sum) // (ncs + 1)): + 920 Q_top.append(Q_sum[i * (ncs + 1) + index_aim][0]) + 921 else: + 922 for i in range(len(Q) // dtr_cnfg): + 923 Q_top.append(Q_sum[dtr_cnfg * i][index_aim]) + 924 if len(Q_top) != len(traj_list) // dtr_cnfg: + 925 raise Exception("qtops and traj_list dont have the same length") + 926 + 927 if kwargs.get('integer_charge', False): + 928 Q_top = [round(q) for q in Q_top] + 929 + 930 truncated_file = file[:-len(postfix)] + 931 + 932 if "names" not in kwargs: + 933 try: + 934 idx = truncated_file.index('r') + 935 except Exception: + 936 if "names" not in kwargs: + 937 raise Exception("Automatic recognition of replicum failed, please enter the key word 'names'.") + 938 ens_name = truncated_file[:idx] + 939 rep_names.append(ens_name + '|' + truncated_file[idx:]) + 940 else: + 941 names = kwargs.get("names") + 942 rep_names = names + 943 deltas.append(Q_top) + 944 + 945 idl = [range(int(configlist[rep][r_start_index[rep]]), int(configlist[rep][r_stop_index[rep]]) + 1, 1) for rep in range(len(deltas))] + 946 deltas = [deltas[nrep][r_start_index[nrep]:r_stop_index[nrep] + 1] for nrep in range(len(deltas))] + 947 result = Obs(deltas, rep_names, idl=idl) + 948 result.tag = {"T": tmax - 1, + 949 "L": L} + 950 return result + 951 + 952 + 953def qtop_projection(qtop, target=0): + 954 """Returns the projection to the topological charge sector defined by target. + 955 + 956 Parameters + 957 ---------- + 958 path : Obs + 959 Topological charge. + 960 target : int + 961 Specifies the topological sector to be reweighted to (default 0) 962 - 963def read_qtop_sector(path, prefix, c, target=0, **kwargs): - 964 """Constructs reweighting factors to a specified topological sector. - 965 - 966 Parameters - 967 ---------- - 968 path : str - 969 path of the measurement files - 970 prefix : str - 971 prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat - 972 c : double - 973 Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L - 974 target : int - 975 Specifies the topological sector to be reweighted to (default 0) - 976 dtr_cnfg : int - 977 (optional) parameter that specifies the number of trajectories - 978 between two configs. - 979 if it is not set, the distance between two measurements - 980 in the file is assumed to be the distance between two configurations. - 981 steps : int - 982 (optional) Distance between two configurations in units of trajectories / - 983 cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given - 984 version : str - 985 version string of the openQCD (sfqcd) version used to create - 986 the ensemble. Default is 2.0. May also be set to sfqcd. - 987 L : int - 988 spatial length of the lattice in L/a. - 989 HAS to be set if version != sfqcd, since openQCD does not provide - 990 this in the header - 991 r_start : list - 992 offset of the first ensemble, making it easier to match - 993 later on with other Obs - 994 r_stop : list - 995 last configurations that need to be read (per replicum) - 996 files : list - 997 specify the exact files that need to be read - 998 from path, practical if e.g. only one replicum is needed - 999 names : list -1000 Alternative labeling for replicas/ensembles. -1001 Has to have the appropriate length -1002 Zeuthen_flow : bool -1003 (optional) If True, the Zeuthen flow is used for Qtop. Only possible -1004 for version=='sfqcd' If False, the Wilson flow is used. -1005 -1006 Returns -1007 ------- -1008 reto : Obs -1009 projection to the topological charge sector defined by target -1010 """ -1011 -1012 if not isinstance(target, int): -1013 raise Exception("'target' has to be an integer.") -1014 -1015 kwargs['integer_charge'] = True -1016 qtop = read_qtop(path, prefix, c, **kwargs) -1017 -1018 return qtop_projection(qtop, target=target) -1019 -1020 -1021def read_ms5_xsf(path, prefix, qc, corr, sep="r", **kwargs): -1022 """ -1023 Read data from files in the specified directory with the specified prefix and quark combination extension, and return a `Corr` object containing the data. -1024 -1025 Parameters -1026 ---------- -1027 path : str -1028 The directory to search for the files in. -1029 prefix : str -1030 The prefix to match the files against. -1031 qc : str -1032 The quark combination extension to match the files against. -1033 corr : str -1034 The correlator to extract data for. -1035 sep : str, optional -1036 The separator to use when parsing the replika names. -1037 **kwargs -1038 Additional keyword arguments. The following keyword arguments are recognized: -1039 -1040 - names (List[str]): A list of names to use for the replicas. -1041 -1042 Returns -1043 ------- -1044 Corr -1045 A complex valued `Corr` object containing the data read from the files. In case of boudary to bulk correlators. -1046 or -1047 CObs -1048 A complex valued `CObs` object containing the data read from the files. In case of boudary to boundary correlators. -1049 -1050 -1051 Raises -1052 ------ -1053 FileNotFoundError -1054 If no files matching the specified prefix and quark combination extension are found in the specified directory. -1055 IOError -1056 If there is an error reading a file. -1057 struct.error -1058 If there is an error unpacking binary data. -1059 """ -1060 -1061 found = [] -1062 files = [] -1063 names = [] -1064 -1065 if "names" in kwargs: -1066 names = kwargs.get("names") -1067 -1068 for (dirpath, dirnames, filenames) in os.walk(path + "/"): -1069 found.extend(filenames) -1070 break -1071 -1072 for f in found: -1073 if fnmatch.fnmatch(f, prefix + "*.ms5_xsf_" + qc + ".dat"): -1074 files.append(f) -1075 if "names" not in kwargs: -1076 if not sep == "": -1077 se = f.split(".")[0] -1078 for s in f.split(".")[1:-1]: -1079 se += "." + s -1080 names.append(se.split(sep)[0] + "|r" + se.split(sep)[1]) -1081 else: -1082 names.append(prefix) -1083 -1084 names = sorted(names) -1085 files = sorted(files) -1086 -1087 cnfgs = [] -1088 realsamples = [] -1089 imagsamples = [] -1090 repnum = 0 -1091 for file in files: -1092 with open(path + "/" + file, "rb") as fp: + 963 Returns + 964 ------- + 965 reto : Obs + 966 projection to the topological charge sector defined by target + 967 """ + 968 if qtop.reweighted: + 969 raise Exception('You can not use a reweighted observable for reweighting!') + 970 + 971 proj_qtop = [] + 972 for n in qtop.deltas: + 973 proj_qtop.append(np.array([1 if round(qtop.r_values[n] + q) == target else 0 for q in qtop.deltas[n]])) + 974 + 975 reto = Obs(proj_qtop, qtop.names, idl=[qtop.idl[name] for name in qtop.names]) + 976 return reto + 977 + 978 + 979def read_qtop_sector(path, prefix, c, target=0, **kwargs): + 980 """Constructs reweighting factors to a specified topological sector. + 981 + 982 Parameters + 983 ---------- + 984 path : str + 985 path of the measurement files + 986 prefix : str + 987 prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat + 988 c : double + 989 Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L + 990 target : int + 991 Specifies the topological sector to be reweighted to (default 0) + 992 dtr_cnfg : int + 993 (optional) parameter that specifies the number of trajectories + 994 between two configs. + 995 if it is not set, the distance between two measurements + 996 in the file is assumed to be the distance between two configurations. + 997 steps : int + 998 (optional) Distance between two configurations in units of trajectories / + 999 cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given +1000 version : str +1001 version string of the openQCD (sfqcd) version used to create +1002 the ensemble. Default is 2.0. May also be set to sfqcd. +1003 L : int +1004 spatial length of the lattice in L/a. +1005 HAS to be set if version != sfqcd, since openQCD does not provide +1006 this in the header +1007 r_start : list +1008 offset of the first ensemble, making it easier to match +1009 later on with other Obs +1010 r_stop : list +1011 last configurations that need to be read (per replicum) +1012 files : list +1013 specify the exact files that need to be read +1014 from path, practical if e.g. only one replicum is needed +1015 names : list +1016 Alternative labeling for replicas/ensembles. +1017 Has to have the appropriate length +1018 Zeuthen_flow : bool +1019 (optional) If True, the Zeuthen flow is used for Qtop. Only possible +1020 for version=='sfqcd' If False, the Wilson flow is used. +1021 +1022 Returns +1023 ------- +1024 reto : Obs +1025 projection to the topological charge sector defined by target +1026 """ +1027 +1028 if not isinstance(target, int): +1029 raise Exception("'target' has to be an integer.") +1030 +1031 kwargs['integer_charge'] = True +1032 qtop = read_qtop(path, prefix, c, **kwargs) +1033 +1034 return qtop_projection(qtop, target=target) +1035 +1036 +1037def read_ms5_xsf(path, prefix, qc, corr, sep="r", **kwargs): +1038 """ +1039 Read data from files in the specified directory with the specified prefix and quark combination extension, and return a `Corr` object containing the data. +1040 +1041 Parameters +1042 ---------- +1043 path : str +1044 The directory to search for the files in. +1045 prefix : str +1046 The prefix to match the files against. +1047 qc : str +1048 The quark combination extension to match the files against. +1049 corr : str +1050 The correlator to extract data for. +1051 sep : str, optional +1052 The separator to use when parsing the replika names. +1053 **kwargs +1054 Additional keyword arguments. The following keyword arguments are recognized: +1055 +1056 - names (List[str]): A list of names to use for the replicas. +1057 +1058 Returns +1059 ------- +1060 Corr +1061 A complex valued `Corr` object containing the data read from the files. In case of boudary to bulk correlators. +1062 or +1063 CObs +1064 A complex valued `CObs` object containing the data read from the files. In case of boudary to boundary correlators. +1065 +1066 +1067 Raises +1068 ------ +1069 FileNotFoundError +1070 If no files matching the specified prefix and quark combination extension are found in the specified directory. +1071 IOError +1072 If there is an error reading a file. +1073 struct.error +1074 If there is an error unpacking binary data. +1075 """ +1076 +1077 # found = [] +1078 files = [] +1079 names = [] +1080 +1081 # test if the input is correct +1082 if qc not in ['dd', 'ud', 'du', 'uu']: +1083 raise Exception("Unknown quark conbination!") +1084 +1085 if corr not in ["gS", "gP", "gA", "gV", "gVt", "lA", "lV", "lVt", "lT", "lTt", "g1", "l1"]: +1086 raise Exception("Unknown correlator!") +1087 +1088 if "files" in kwargs: +1089 known_files = kwargs.get("files") +1090 else: +1091 known_files = [] +1092 files = _find_files(path, prefix, "ms5_xsf_" + qc, "dat", known_files=known_files) 1093 -1094 t = fp.read(8) -1095 kappa = struct.unpack('d', t)[0] -1096 t = fp.read(8) -1097 csw = struct.unpack('d', t)[0] -1098 t = fp.read(8) -1099 dF = struct.unpack('d', t)[0] -1100 t = fp.read(8) -1101 zF = struct.unpack('d', t)[0] -1102 -1103 t = fp.read(4) -1104 tmax = struct.unpack('i', t)[0] -1105 t = fp.read(4) -1106 bnd = struct.unpack('i', t)[0] -1107 -1108 placesBI = ["gS", "gP", -1109 "gA", "gV", -1110 "gVt", "lA", -1111 "lV", "lVt", -1112 "lT", "lTt"] -1113 placesBB = ["g1", "l1"] -1114 -1115 # the chunks have the following structure: -1116 # confignumber, 10x timedependent complex correlators as doubles, 2x timeindependent complex correlators as doubles -1117 -1118 chunksize = 4 + (8 * 2 * tmax * 10) + (8 * 2 * 2) -1119 packstr = '=i' + ('d' * 2 * tmax * 10) + ('d' * 2 * 2) -1120 cnfgs.append([]) -1121 realsamples.append([]) -1122 imagsamples.append([]) -1123 for t in range(tmax): -1124 realsamples[repnum].append([]) -1125 imagsamples[repnum].append([]) -1126 -1127 while True: -1128 cnfgt = fp.read(chunksize) -1129 if not cnfgt: -1130 break -1131 asascii = struct.unpack(packstr, cnfgt) -1132 cnfg = asascii[0] -1133 cnfgs[repnum].append(cnfg) -1134 -1135 if corr not in placesBB: -1136 tmpcorr = asascii[1 + 2 * tmax * placesBI.index(corr):1 + 2 * tmax * placesBI.index(corr) + 2 * tmax] -1137 else: -1138 tmpcorr = asascii[1 + 2 * tmax * len(placesBI) + 2 * placesBB.index(corr):1 + 2 * tmax * len(placesBI) + 2 * placesBB.index(corr) + 2] +1094 if "names" in kwargs: +1095 names = kwargs.get("names") +1096 else: +1097 for f in files: +1098 if not sep == "": +1099 se = f.split(".")[0] +1100 for s in f.split(".")[1:-2]: +1101 se += "." + s +1102 names.append(se.split(sep)[0] + "|r" + se.split(sep)[1]) +1103 else: +1104 names.append(prefix) +1105 +1106 names = sorted(names) +1107 files = sorted(files) +1108 +1109 cnfgs = [] +1110 realsamples = [] +1111 imagsamples = [] +1112 repnum = 0 +1113 for file in files: +1114 with open(path + "/" + file, "rb") as fp: +1115 +1116 t = fp.read(8) +1117 kappa = struct.unpack('d', t)[0] +1118 t = fp.read(8) +1119 csw = struct.unpack('d', t)[0] +1120 t = fp.read(8) +1121 dF = struct.unpack('d', t)[0] +1122 t = fp.read(8) +1123 zF = struct.unpack('d', t)[0] +1124 +1125 t = fp.read(4) +1126 tmax = struct.unpack('i', t)[0] +1127 t = fp.read(4) +1128 bnd = struct.unpack('i', t)[0] +1129 +1130 placesBI = ["gS", "gP", +1131 "gA", "gV", +1132 "gVt", "lA", +1133 "lV", "lVt", +1134 "lT", "lTt"] +1135 placesBB = ["g1", "l1"] +1136 +1137 # the chunks have the following structure: +1138 # confignumber, 10x timedependent complex correlators as doubles, 2x timeindependent complex correlators as doubles 1139 -1140 corrres = [[], []] -1141 for i in range(len(tmpcorr)): -1142 corrres[i % 2].append(tmpcorr[i]) -1143 for t in range(int(len(tmpcorr) / 2)): -1144 realsamples[repnum][t].append(corrres[0][t]) -1145 for t in range(int(len(tmpcorr) / 2)): -1146 imagsamples[repnum][t].append(corrres[1][t]) -1147 repnum += 1 +1140 chunksize = 4 + (8 * 2 * tmax * 10) + (8 * 2 * 2) +1141 packstr = '=i' + ('d' * 2 * tmax * 10) + ('d' * 2 * 2) +1142 cnfgs.append([]) +1143 realsamples.append([]) +1144 imagsamples.append([]) +1145 for t in range(tmax): +1146 realsamples[repnum].append([]) +1147 imagsamples[repnum].append([]) 1148 -1149 s = "Read correlator " + corr + " from " + str(repnum) + " replika with " + str(len(realsamples[0][t])) -1150 for rep in range(1, repnum): -1151 s += ", " + str(len(realsamples[rep][t])) -1152 s += " samples" -1153 print(s) -1154 print("Asserted run parameters:\n T:", tmax, "kappa:", kappa, "csw:", csw, "dF:", dF, "zF:", zF, "bnd:", bnd) -1155 -1156 # we have the data now... but we need to re format the whole thing and put it into Corr objects. -1157 -1158 compObs = [] -1159 -1160 for t in range(int(len(tmpcorr) / 2)): -1161 compObs.append(CObs(Obs([realsamples[rep][t] for rep in range(repnum)], names=names, idl=cnfgs), -1162 Obs([imagsamples[rep][t] for rep in range(repnum)], names=names, idl=cnfgs))) -1163 -1164 if len(compObs) == 1: -1165 return compObs[0] -1166 else: -1167 return Corr(compObs) +1149 while True: +1150 cnfgt = fp.read(chunksize) +1151 if not cnfgt: +1152 break +1153 asascii = struct.unpack(packstr, cnfgt) +1154 cnfg = asascii[0] +1155 cnfgs[repnum].append(cnfg) +1156 +1157 if corr not in placesBB: +1158 tmpcorr = asascii[1 + 2 * tmax * placesBI.index(corr):1 + 2 * tmax * placesBI.index(corr) + 2 * tmax] +1159 else: +1160 tmpcorr = asascii[1 + 2 * tmax * len(placesBI) + 2 * placesBB.index(corr):1 + 2 * tmax * len(placesBI) + 2 * placesBB.index(corr) + 2] +1161 +1162 corrres = [[], []] +1163 for i in range(len(tmpcorr)): +1164 corrres[i % 2].append(tmpcorr[i]) +1165 for t in range(int(len(tmpcorr) / 2)): +1166 realsamples[repnum][t].append(corrres[0][t]) +1167 for t in range(int(len(tmpcorr) / 2)): +1168 imagsamples[repnum][t].append(corrres[1][t]) +1169 repnum += 1 +1170 +1171 s = "Read correlator " + corr + " from " + str(repnum) + " replika with " + str(len(realsamples[0][t])) +1172 for rep in range(1, repnum): +1173 s += ", " + str(len(realsamples[rep][t])) +1174 s += " samples" +1175 print(s) +1176 print("Asserted run parameters:\n T:", tmax, "kappa:", kappa, "csw:", csw, "dF:", dF, "zF:", zF, "bnd:", bnd) +1177 +1178 # we have the data now... but we need to re format the whole thing and put it into Corr objects. +1179 +1180 compObs = [] +1181 +1182 for t in range(int(len(tmpcorr) / 2)): +1183 compObs.append(CObs(Obs([realsamples[rep][t] for rep in range(repnum)], names=names, idl=cnfgs), +1184 Obs([imagsamples[rep][t] for rep in range(repnum)], names=names, idl=cnfgs))) +1185 +1186 if len(compObs) == 1: +1187 return compObs[0] +1188 else: +1189 return Corr(compObs) @@ -1276,228 +1298,221 @@ - 
 16def read_rwms(path, prefix, version='2.0', names=None, **kwargs):
- 17    """Read rwms format from given folder structure. Returns a list of length nrw
- 18
- 19    Parameters
- 20    ----------
- 21    path : str
- 22        path that contains the data files
- 23    prefix : str
- 24        all files in path that start with prefix are considered as input files.
- 25        May be used together postfix to consider only special file endings.
- 26        Prefix is ignored, if the keyword 'files' is used.
- 27    version : str
- 28        version of openQCD, default 2.0
- 29    names : list
- 30        list of names that is assigned to the data according according
- 31        to the order in the file list. Use careful, if you do not provide file names!
- 32    r_start : list
- 33        list which contains the first config to be read for each replicum
- 34    r_stop : list
- 35        list which contains the last config to be read for each replicum
- 36    r_step : int
- 37        integer that defines a fixed step size between two measurements (in units of configs)
- 38        If not given, r_step=1 is assumed.
- 39    postfix : str
- 40        postfix of the file to read, e.g. '.ms1' for openQCD-files
- 41    files : list
- 42        list which contains the filenames to be read. No automatic detection of
- 43        files performed if given.
- 44    print_err : bool
- 45        Print additional information that is useful for debugging.
- 46
- 47    Returns
- 48    -------
- 49    rwms : Obs
- 50        Reweighting factors read
- 51    """
- 52    known_oqcd_versions = ['1.4', '1.6', '2.0']
- 53    if not (version in known_oqcd_versions):
- 54        raise Exception('Unknown openQCD version defined!')
- 55    print("Working with openQCD version " + version)
- 56    if 'postfix' in kwargs:
- 57        postfix = kwargs.get('postfix')
- 58    else:
- 59        postfix = ''
- 60    ls = []
- 61    for (dirpath, dirnames, filenames) in os.walk(path):
- 62        ls.extend(filenames)
- 63        break
- 64
- 65    if not ls:
- 66        raise Exception(f"Error, directory '{path}' not found")
- 67    if 'files' in kwargs:
- 68        ls = kwargs.get('files')
- 69    else:
- 70        for exc in ls:
- 71            if not fnmatch.fnmatch(exc, prefix + '*' + postfix + '.dat'):
- 72                ls = list(set(ls) - set([exc]))
- 73        if len(ls) > 1:
- 74            ls.sort(key=lambda x: int(re.findall(r'\d+', x[len(prefix):])[0]))
- 75    replica = len(ls)
- 76
- 77    if 'r_start' in kwargs:
- 78        r_start = kwargs.get('r_start')
- 79        if len(r_start) != replica:
- 80            raise Exception('r_start does not match number of replicas')
- 81        r_start = [o if o else None for o in r_start]
- 82    else:
- 83        r_start = [None] * replica
- 84
- 85    if 'r_stop' in kwargs:
- 86        r_stop = kwargs.get('r_stop')
- 87        if len(r_stop) != replica:
- 88            raise Exception('r_stop does not match number of replicas')
- 89    else:
- 90        r_stop = [None] * replica
- 91
- 92    if 'r_step' in kwargs:
- 93        r_step = kwargs.get('r_step')
- 94    else:
- 95        r_step = 1
- 96
- 97    print('Read reweighting factors from', prefix[:-1], ',',
- 98          replica, 'replica', end='')
- 99
-100    if names is None:
-101        rep_names = []
-102        for entry in ls:
-103            truncated_entry = entry
-104            suffixes = [".dat", ".rwms", ".ms1"]
-105            for suffix in suffixes:
-106                if truncated_entry.endswith(suffix):
-107                    truncated_entry = truncated_entry[0:-len(suffix)]
-108            idx = truncated_entry.index('r')
-109            rep_names.append(truncated_entry[:idx] + '|' + truncated_entry[idx:])
-110    else:
-111        rep_names = names
-112
-113    print_err = 0
-114    if 'print_err' in kwargs:
-115        print_err = 1
-116        print()
+            
 56def read_rwms(path, prefix, version='2.0', names=None, **kwargs):
+ 57    """Read rwms format from given folder structure. Returns a list of length nrw
+ 58
+ 59    Parameters
+ 60    ----------
+ 61    path : str
+ 62        path that contains the data files
+ 63    prefix : str
+ 64        all files in path that start with prefix are considered as input files.
+ 65        May be used together postfix to consider only special file endings.
+ 66        Prefix is ignored, if the keyword 'files' is used.
+ 67    version : str
+ 68        version of openQCD, default 2.0
+ 69    names : list
+ 70        list of names that is assigned to the data according according
+ 71        to the order in the file list. Use careful, if you do not provide file names!
+ 72    r_start : list
+ 73        list which contains the first config to be read for each replicum
+ 74    r_stop : list
+ 75        list which contains the last config to be read for each replicum
+ 76    r_step : int
+ 77        integer that defines a fixed step size between two measurements (in units of configs)
+ 78        If not given, r_step=1 is assumed.
+ 79    postfix : str
+ 80        postfix of the file to read, e.g. '.ms1' for openQCD-files
+ 81    files : list
+ 82        list which contains the filenames to be read. No automatic detection of
+ 83        files performed if given.
+ 84    print_err : bool
+ 85        Print additional information that is useful for debugging.
+ 86
+ 87    Returns
+ 88    -------
+ 89    rwms : Obs
+ 90        Reweighting factors read
+ 91    """
+ 92    known_oqcd_versions = ['1.4', '1.6', '2.0']
+ 93    if not (version in known_oqcd_versions):
+ 94        raise Exception('Unknown openQCD version defined!')
+ 95    print("Working with openQCD version " + version)
+ 96    if 'postfix' in kwargs:
+ 97        postfix = kwargs.get('postfix')
+ 98    else:
+ 99        postfix = ''
+100
+101    if 'files' in kwargs:
+102        known_files = kwargs.get('files')
+103    else:
+104        known_files = []
+105
+106    ls = _find_files(path, prefix, postfix, 'dat', known_files=known_files)
+107
+108    replica = len(ls)
+109
+110    if 'r_start' in kwargs:
+111        r_start = kwargs.get('r_start')
+112        if len(r_start) != replica:
+113            raise Exception('r_start does not match number of replicas')
+114        r_start = [o if o else None for o in r_start]
+115    else:
+116        r_start = [None] * replica
 117
-118    deltas = []
-119
-120    configlist = []
-121    r_start_index = []
-122    r_stop_index = []
-123
-124    for rep in range(replica):
-125        tmp_array = []
-126        with open(path + '/' + ls[rep], 'rb') as fp:
-127
-128            t = fp.read(4)  # number of reweighting factors
-129            if rep == 0:
-130                nrw = struct.unpack('i', t)[0]
-131                if version == '2.0':
-132                    nrw = int(nrw / 2)
-133                for k in range(nrw):
-134                    deltas.append([])
-135            else:
-136                if ((nrw != struct.unpack('i', t)[0] and (not version == '2.0')) or (nrw != struct.unpack('i', t)[0] / 2 and version == '2.0')):
-137                    raise Exception('Error: different number of reweighting factors for replicum', rep)
-138
-139            for k in range(nrw):
-140                tmp_array.append([])
-141
-142            # This block is necessary for openQCD1.6 and openQCD2.0 ms1 files
-143            nfct = []
-144            if version in ['1.6', '2.0']:
-145                for i in range(nrw):
-146                    t = fp.read(4)
-147                    nfct.append(struct.unpack('i', t)[0])
-148            else:
-149                for i in range(nrw):
-150                    nfct.append(1)
-151
-152            nsrc = []
-153            for i in range(nrw):
-154                t = fp.read(4)
-155                nsrc.append(struct.unpack('i', t)[0])
-156            if version == '2.0':
-157                if not struct.unpack('i', fp.read(4))[0] == 0:
-158                    print('something is wrong!')
-159
-160            configlist.append([])
-161            while True:
-162                t = fp.read(4)
-163                if len(t) < 4:
-164                    break
-165                config_no = struct.unpack('i', t)[0]
-166                configlist[-1].append(config_no)
-167                for i in range(nrw):
-168                    if (version == '2.0'):
-169                        tmpd = _read_array_openQCD2(fp)
-170                        tmpd = _read_array_openQCD2(fp)
-171                        tmp_rw = tmpd['arr']
-172                        tmp_nfct = 1.0
-173                        for j in range(tmpd['n'][0]):
-174                            tmp_nfct *= np.mean(np.exp(-np.asarray(tmp_rw[j])))
-175                            if print_err:
-176                                print(config_no, i, j,
-177                                      np.mean(np.exp(-np.asarray(tmp_rw[j]))),
-178                                      np.std(np.exp(-np.asarray(tmp_rw[j]))))
-179                                print('Sources:',
-180                                      np.exp(-np.asarray(tmp_rw[j])))
-181                                print('Partial factor:', tmp_nfct)
-182                    elif version == '1.6' or version == '1.4':
-183                        tmp_nfct = 1.0
-184                        for j in range(nfct[i]):
-185                            t = fp.read(8 * nsrc[i])
-186                            t = fp.read(8 * nsrc[i])
-187                            tmp_rw = struct.unpack('d' * nsrc[i], t)
-188                            tmp_nfct *= np.mean(np.exp(-np.asarray(tmp_rw)))
-189                            if print_err:
-190                                print(config_no, i, j,
-191                                      np.mean(np.exp(-np.asarray(tmp_rw))),
-192                                      np.std(np.exp(-np.asarray(tmp_rw))))
-193                                print('Sources:', np.exp(-np.asarray(tmp_rw)))
-194                                print('Partial factor:', tmp_nfct)
-195                    tmp_array[i].append(tmp_nfct)
-196
-197            diffmeas = configlist[-1][-1] - configlist[-1][-2]
-198            configlist[-1] = [item // diffmeas for item in configlist[-1]]
-199            if configlist[-1][0] > 1 and diffmeas > 1:
-200                warnings.warn('Assume thermalization and that the first measurement belongs to the first config.')
-201                offset = configlist[-1][0] - 1
-202                configlist[-1] = [item - offset for item in configlist[-1]]
-203
-204            if r_start[rep] is None:
-205                r_start_index.append(0)
-206            else:
-207                try:
-208                    r_start_index.append(configlist[-1].index(r_start[rep]))
-209                except ValueError:
-210                    raise Exception('Config %d not in file with range [%d, %d]' % (
-211                        r_start[rep], configlist[-1][0], configlist[-1][-1])) from None
-212
-213            if r_stop[rep] is None:
-214                r_stop_index.append(len(configlist[-1]) - 1)
-215            else:
-216                try:
-217                    r_stop_index.append(configlist[-1].index(r_stop[rep]))
-218                except ValueError:
-219                    raise Exception('Config %d not in file with range [%d, %d]' % (
-220                        r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None
-221
-222            for k in range(nrw):
-223                deltas[k].append(tmp_array[k][r_start_index[rep]:r_stop_index[rep] + 1][::r_step])
-224
-225    if np.any([len(np.unique(np.diff(cl))) != 1 for cl in configlist]):
-226        raise Exception('Irregular spaced data in input file!', [len(np.unique(np.diff(cl))) for cl in configlist])
-227    stepsizes = [list(np.unique(np.diff(cl)))[0] for cl in configlist]
-228    if np.any([step != 1 for step in stepsizes]):
-229        warnings.warn('Stepsize between configurations is greater than one!' + str(stepsizes), RuntimeWarning)
-230
-231    print(',', nrw, 'reweighting factors with', nsrc, 'sources')
-232    result = []
-233    idl = [range(configlist[rep][r_start_index[rep]], configlist[rep][r_stop_index[rep]] + 1, r_step) for rep in range(replica)]
-234
-235    for t in range(nrw):
-236        result.append(Obs(deltas[t], rep_names, idl=idl))
-237    return result
+118    if 'r_stop' in kwargs:
+119        r_stop = kwargs.get('r_stop')
+120        if len(r_stop) != replica:
+121            raise Exception('r_stop does not match number of replicas')
+122    else:
+123        r_stop = [None] * replica
+124
+125    if 'r_step' in kwargs:
+126        r_step = kwargs.get('r_step')
+127    else:
+128        r_step = 1
+129
+130    print('Read reweighting factors from', prefix[:-1], ',',
+131          replica, 'replica', end='')
+132
+133    if names is None:
+134        rep_names = []
+135        for entry in ls:
+136            truncated_entry = entry
+137            suffixes = [".dat", ".rwms", ".ms1"]
+138            for suffix in suffixes:
+139                if truncated_entry.endswith(suffix):
+140                    truncated_entry = truncated_entry[0:-len(suffix)]
+141            idx = truncated_entry.index('r')
+142            rep_names.append(truncated_entry[:idx] + '|' + truncated_entry[idx:])
+143    else:
+144        rep_names = names
+145
+146    print_err = 0
+147    if 'print_err' in kwargs:
+148        print_err = 1
+149        print()
+150
+151    deltas = []
+152
+153    configlist = []
+154    r_start_index = []
+155    r_stop_index = []
+156
+157    for rep in range(replica):
+158        tmp_array = []
+159        with open(path + '/' + ls[rep], 'rb') as fp:
+160
+161            t = fp.read(4)  # number of reweighting factors
+162            if rep == 0:
+163                nrw = struct.unpack('i', t)[0]
+164                if version == '2.0':
+165                    nrw = int(nrw / 2)
+166                for k in range(nrw):
+167                    deltas.append([])
+168            else:
+169                if ((nrw != struct.unpack('i', t)[0] and (not version == '2.0')) or (nrw != struct.unpack('i', t)[0] / 2 and version == '2.0')):
+170                    raise Exception('Error: different number of reweighting factors for replicum', rep)
+171
+172            for k in range(nrw):
+173                tmp_array.append([])
+174
+175            # This block is necessary for openQCD1.6 and openQCD2.0 ms1 files
+176            nfct = []
+177            if version in ['1.6', '2.0']:
+178                for i in range(nrw):
+179                    t = fp.read(4)
+180                    nfct.append(struct.unpack('i', t)[0])
+181            else:
+182                for i in range(nrw):
+183                    nfct.append(1)
+184
+185            nsrc = []
+186            for i in range(nrw):
+187                t = fp.read(4)
+188                nsrc.append(struct.unpack('i', t)[0])
+189            if version == '2.0':
+190                if not struct.unpack('i', fp.read(4))[0] == 0:
+191                    raise Exception("You are using the input for openQCD version 2.0, this is not correct.")
+192
+193            configlist.append([])
+194            while True:
+195                t = fp.read(4)
+196                if len(t) < 4:
+197                    break
+198                config_no = struct.unpack('i', t)[0]
+199                configlist[-1].append(config_no)
+200                for i in range(nrw):
+201                    if (version == '2.0'):
+202                        tmpd = _read_array_openQCD2(fp)
+203                        tmpd = _read_array_openQCD2(fp)
+204                        tmp_rw = tmpd['arr']
+205                        tmp_nfct = 1.0
+206                        for j in range(tmpd['n'][0]):
+207                            tmp_nfct *= np.mean(np.exp(-np.asarray(tmp_rw[j])))
+208                            if print_err:
+209                                print(config_no, i, j,
+210                                      np.mean(np.exp(-np.asarray(tmp_rw[j]))),
+211                                      np.std(np.exp(-np.asarray(tmp_rw[j]))))
+212                                print('Sources:',
+213                                      np.exp(-np.asarray(tmp_rw[j])))
+214                                print('Partial factor:', tmp_nfct)
+215                    elif version == '1.6' or version == '1.4':
+216                        tmp_nfct = 1.0
+217                        for j in range(nfct[i]):
+218                            t = fp.read(8 * nsrc[i])
+219                            t = fp.read(8 * nsrc[i])
+220                            tmp_rw = struct.unpack('d' * nsrc[i], t)
+221                            tmp_nfct *= np.mean(np.exp(-np.asarray(tmp_rw)))
+222                            if print_err:
+223                                print(config_no, i, j,
+224                                      np.mean(np.exp(-np.asarray(tmp_rw))),
+225                                      np.std(np.exp(-np.asarray(tmp_rw))))
+226                                print('Sources:', np.exp(-np.asarray(tmp_rw)))
+227                                print('Partial factor:', tmp_nfct)
+228                    tmp_array[i].append(tmp_nfct)
+229
+230            diffmeas = configlist[-1][-1] - configlist[-1][-2]
+231            configlist[-1] = [item // diffmeas for item in configlist[-1]]
+232            if configlist[-1][0] > 1 and diffmeas > 1:
+233                warnings.warn('Assume thermalization and that the first measurement belongs to the first config.')
+234                offset = configlist[-1][0] - 1
+235                configlist[-1] = [item - offset for item in configlist[-1]]
+236
+237            if r_start[rep] is None:
+238                r_start_index.append(0)
+239            else:
+240                try:
+241                    r_start_index.append(configlist[-1].index(r_start[rep]))
+242                except ValueError:
+243                    raise Exception('Config %d not in file with range [%d, %d]' % (
+244                        r_start[rep], configlist[-1][0], configlist[-1][-1])) from None
+245
+246            if r_stop[rep] is None:
+247                r_stop_index.append(len(configlist[-1]) - 1)
+248            else:
+249                try:
+250                    r_stop_index.append(configlist[-1].index(r_stop[rep]))
+251                except ValueError:
+252                    raise Exception('Config %d not in file with range [%d, %d]' % (
+253                        r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None
+254
+255            for k in range(nrw):
+256                deltas[k].append(tmp_array[k][r_start_index[rep]:r_stop_index[rep] + 1][::r_step])
+257
+258    if np.any([len(np.unique(np.diff(cl))) != 1 for cl in configlist]):
+259        raise Exception('Irregular spaced data in input file!', [len(np.unique(np.diff(cl))) for cl in configlist])
+260    stepsizes = [list(np.unique(np.diff(cl)))[0] for cl in configlist]
+261    if np.any([step != 1 for step in stepsizes]):
+262        warnings.warn('Stepsize between configurations is greater than one!' + str(stepsizes), RuntimeWarning)
+263
+264    print(',', nrw, 'reweighting factors with', nsrc, 'sources')
+265    result = []
+266    idl = [range(configlist[rep][r_start_index[rep]], configlist[rep][r_stop_index[rep]] + 1, r_step) for rep in range(replica)]
+267
+268    for t in range(nrw):
+269        result.append(Obs(deltas[t], rep_names, idl=idl))
+270    return result
 

 
 
@@ -1554,255 +1569,246 @@ Reweighting factors read
- 
240def extract_t0(path, prefix, dtr_read, xmin, spatial_extent, fit_range=5, **kwargs):
-241    """Extract t0 from given .ms.dat files. Returns t0 as Obs.
-242
-243    It is assumed that all boundary effects have
-244    sufficiently decayed at x0=xmin.
-245    The data around the zero crossing of t^2<E> - 0.3
-246    is fitted with a linear function
-247    from which the exact root is extracted.
-248
-249    It is assumed that one measurement is performed for each config.
-250    If this is not the case, the resulting idl, as well as the handling
-251    of r_start, r_stop and r_step is wrong and the user has to correct
-252    this in the resulting observable.
-253
-254    Parameters
-255    ----------
-256    path : str
-257        Path to .ms.dat files
-258    prefix : str
-259        Ensemble prefix
-260    dtr_read : int
-261        Determines how many trajectories should be skipped
-262        when reading the ms.dat files.
-263        Corresponds to dtr_cnfg / dtr_ms in the openQCD input file.
-264    xmin : int
-265        First timeslice where the boundary
-266        effects have sufficiently decayed.
-267    spatial_extent : int
-268        spatial extent of the lattice, required for normalization.
-269    fit_range : int
-270        Number of data points left and right of the zero
-271        crossing to be included in the linear fit. (Default: 5)
-272    r_start : list
-273        list which contains the first config to be read for each replicum.
-274    r_stop : list
-275        list which contains the last config to be read for each replicum.
-276    r_step : int
-277        integer that defines a fixed step size between two measurements (in units of configs)
-278        If not given, r_step=1 is assumed.
-279    plaquette : bool
-280        If true extract the plaquette estimate of t0 instead.
-281    names : list
-282        list of names that is assigned to the data according according
-283        to the order in the file list. Use careful, if you do not provide file names!
-284    files : list
-285        list which contains the filenames to be read. No automatic detection of
-286        files performed if given.
-287    plot_fit : bool
-288        If true, the fit for the extraction of t0 is shown together with the data.
-289    assume_thermalization : bool
-290        If True: If the first record divided by the distance between two measurements is larger than
-291        1, it is assumed that this is due to thermalization and the first measurement belongs
-292        to the first config (default).
-293        If False: The config numbers are assumed to be traj_number // difference
-294
-295    Returns
-296    -------
-297    t0 : Obs
-298        Extracted t0
-299    """
-300
-301    ls = []
-302    for (dirpath, dirnames, filenames) in os.walk(path):
-303        ls.extend(filenames)
-304        break
-305
-306    if not ls:
-307        raise Exception('Error, directory not found')
-308
-309    if 'files' in kwargs:
-310        ls = kwargs.get('files')
-311    else:
-312        for exc in ls:
-313            if not fnmatch.fnmatch(exc, prefix + '*.ms.dat'):
-314                ls = list(set(ls) - set([exc]))
-315        if len(ls) > 1:
-316            ls.sort(key=lambda x: int(re.findall(r'\d+', x[len(prefix):])[0]))
-317    replica = len(ls)
-318
-319    if 'r_start' in kwargs:
-320        r_start = kwargs.get('r_start')
-321        if len(r_start) != replica:
-322            raise Exception('r_start does not match number of replicas')
-323        r_start = [o if o else None for o in r_start]
-324    else:
-325        r_start = [None] * replica
-326
-327    if 'r_stop' in kwargs:
-328        r_stop = kwargs.get('r_stop')
-329        if len(r_stop) != replica:
-330            raise Exception('r_stop does not match number of replicas')
-331    else:
-332        r_stop = [None] * replica
+            
273def extract_t0(path, prefix, dtr_read, xmin, spatial_extent, fit_range=5, **kwargs):
+274    """Extract t0 from given .ms.dat files. Returns t0 as Obs.
+275
+276    It is assumed that all boundary effects have
+277    sufficiently decayed at x0=xmin.
+278    The data around the zero crossing of t^2<E> - 0.3
+279    is fitted with a linear function
+280    from which the exact root is extracted.
+281
+282    It is assumed that one measurement is performed for each config.
+283    If this is not the case, the resulting idl, as well as the handling
+284    of r_start, r_stop and r_step is wrong and the user has to correct
+285    this in the resulting observable.
+286
+287    Parameters
+288    ----------
+289    path : str
+290        Path to .ms.dat files
+291    prefix : str
+292        Ensemble prefix
+293    dtr_read : int
+294        Determines how many trajectories should be skipped
+295        when reading the ms.dat files.
+296        Corresponds to dtr_cnfg / dtr_ms in the openQCD input file.
+297    xmin : int
+298        First timeslice where the boundary
+299        effects have sufficiently decayed.
+300    spatial_extent : int
+301        spatial extent of the lattice, required for normalization.
+302    fit_range : int
+303        Number of data points left and right of the zero
+304        crossing to be included in the linear fit. (Default: 5)
+305    r_start : list
+306        list which contains the first config to be read for each replicum.
+307    r_stop : list
+308        list which contains the last config to be read for each replicum.
+309    r_step : int
+310        integer that defines a fixed step size between two measurements (in units of configs)
+311        If not given, r_step=1 is assumed.
+312    plaquette : bool
+313        If true extract the plaquette estimate of t0 instead.
+314    names : list
+315        list of names that is assigned to the data according according
+316        to the order in the file list. Use careful, if you do not provide file names!
+317    files : list
+318        list which contains the filenames to be read. No automatic detection of
+319        files performed if given.
+320    plot_fit : bool
+321        If true, the fit for the extraction of t0 is shown together with the data.
+322    assume_thermalization : bool
+323        If True: If the first record divided by the distance between two measurements is larger than
+324        1, it is assumed that this is due to thermalization and the first measurement belongs
+325        to the first config (default).
+326        If False: The config numbers are assumed to be traj_number // difference
+327
+328    Returns
+329    -------
+330    t0 : Obs
+331        Extracted t0
+332    """
 333
-334    if 'r_step' in kwargs:
-335        r_step = kwargs.get('r_step')
+334    if 'files' in kwargs:
+335        known_files = kwargs.get('files')
 336    else:
-337        r_step = 1
+337        known_files = []
 338
-339    print('Extract t0 from', prefix, ',', replica, 'replica')
+339    ls = _find_files(path, prefix, 'ms', 'dat', known_files=known_files)
 340
-341    if 'names' in kwargs:
-342        rep_names = kwargs.get('names')
-343    else:
-344        rep_names = []
-345        for entry in ls:
-346            truncated_entry = entry.split('.')[0]
-347            idx = truncated_entry.index('r')
-348            rep_names.append(truncated_entry[:idx] + '|' + truncated_entry[idx:])
-349
-350    Ysum = []
-351
-352    configlist = []
-353    r_start_index = []
-354    r_stop_index = []
-355
-356    for rep in range(replica):
+341    replica = len(ls)
+342
+343    if 'r_start' in kwargs:
+344        r_start = kwargs.get('r_start')
+345        if len(r_start) != replica:
+346            raise Exception('r_start does not match number of replicas')
+347        r_start = [o if o else None for o in r_start]
+348    else:
+349        r_start = [None] * replica
+350
+351    if 'r_stop' in kwargs:
+352        r_stop = kwargs.get('r_stop')
+353        if len(r_stop) != replica:
+354            raise Exception('r_stop does not match number of replicas')
+355    else:
+356        r_stop = [None] * replica
 357
-358        with open(path + '/' + ls[rep], 'rb') as fp:
-359            t = fp.read(12)
-360            header = struct.unpack('iii', t)
-361            if rep == 0:
-362                dn = header[0]
-363                nn = header[1]
-364                tmax = header[2]
-365            elif dn != header[0] or nn != header[1] or tmax != header[2]:
-366                raise Exception('Replica parameters do not match.')
-367
-368            t = fp.read(8)
-369            if rep == 0:
-370                eps = struct.unpack('d', t)[0]
-371                print('Step size:', eps, ', Maximal t value:', dn * (nn) * eps)
-372            elif eps != struct.unpack('d', t)[0]:
-373                raise Exception('Values for eps do not match among replica.')
-374
-375            Ysl = []
-376
-377            configlist.append([])
-378            while True:
-379                t = fp.read(4)
-380                if (len(t) < 4):
-381                    break
-382                nc = struct.unpack('i', t)[0]
-383                configlist[-1].append(nc)
-384
-385                t = fp.read(8 * tmax * (nn + 1))
-386                if kwargs.get('plaquette'):
-387                    if nc % dtr_read == 0:
-388                        Ysl.append(struct.unpack('d' * tmax * (nn + 1), t))
-389                t = fp.read(8 * tmax * (nn + 1))
-390                if not kwargs.get('plaquette'):
-391                    if nc % dtr_read == 0:
-392                        Ysl.append(struct.unpack('d' * tmax * (nn + 1), t))
-393                t = fp.read(8 * tmax * (nn + 1))
-394
-395        Ysum.append([])
-396        for i, item in enumerate(Ysl):
-397            Ysum[-1].append([np.mean(item[current + xmin:
-398                             current + tmax - xmin])
-399                            for current in range(0, len(item), tmax)])
+358    if 'r_step' in kwargs:
+359        r_step = kwargs.get('r_step')
+360    else:
+361        r_step = 1
+362
+363    print('Extract t0 from', prefix, ',', replica, 'replica')
+364
+365    if 'names' in kwargs:
+366        rep_names = kwargs.get('names')
+367    else:
+368        rep_names = []
+369        for entry in ls:
+370            truncated_entry = entry.split('.')[0]
+371            idx = truncated_entry.index('r')
+372            rep_names.append(truncated_entry[:idx] + '|' + truncated_entry[idx:])
+373
+374    Ysum = []
+375
+376    configlist = []
+377    r_start_index = []
+378    r_stop_index = []
+379
+380    for rep in range(replica):
+381
+382        with open(path + '/' + ls[rep], 'rb') as fp:
+383            t = fp.read(12)
+384            header = struct.unpack('iii', t)
+385            if rep == 0:
+386                dn = header[0]
+387                nn = header[1]
+388                tmax = header[2]
+389            elif dn != header[0] or nn != header[1] or tmax != header[2]:
+390                raise Exception('Replica parameters do not match.')
+391
+392            t = fp.read(8)
+393            if rep == 0:
+394                eps = struct.unpack('d', t)[0]
+395                print('Step size:', eps, ', Maximal t value:', dn * (nn) * eps)
+396            elif eps != struct.unpack('d', t)[0]:
+397                raise Exception('Values for eps do not match among replica.')
+398
+399            Ysl = []
 400
-401        diffmeas = configlist[-1][-1] - configlist[-1][-2]
-402        configlist[-1] = [item // diffmeas for item in configlist[-1]]
-403        if kwargs.get('assume_thermalization', True) and configlist[-1][0] > 1:
-404            warnings.warn('Assume thermalization and that the first measurement belongs to the first config.')
-405            offset = configlist[-1][0] - 1
-406            configlist[-1] = [item - offset for item in configlist[-1]]
-407
-408        if r_start[rep] is None:
-409            r_start_index.append(0)
-410        else:
-411            try:
-412                r_start_index.append(configlist[-1].index(r_start[rep]))
-413            except ValueError:
-414                raise Exception('Config %d not in file with range [%d, %d]' % (
-415                    r_start[rep], configlist[-1][0], configlist[-1][-1])) from None
-416
-417        if r_stop[rep] is None:
-418            r_stop_index.append(len(configlist[-1]) - 1)
-419        else:
-420            try:
-421                r_stop_index.append(configlist[-1].index(r_stop[rep]))
-422            except ValueError:
-423                raise Exception('Config %d not in file with range [%d, %d]' % (
-424                    r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None
-425
-426    if np.any([len(np.unique(np.diff(cl))) != 1 for cl in configlist]):
-427        raise Exception('Irregular spaced data in input file!', [len(np.unique(np.diff(cl))) for cl in configlist])
-428    stepsizes = [list(np.unique(np.diff(cl)))[0] for cl in configlist]
-429    if np.any([step != 1 for step in stepsizes]):
-430        warnings.warn('Stepsize between configurations is greater than one!' + str(stepsizes), RuntimeWarning)
+401            configlist.append([])
+402            while True:
+403                t = fp.read(4)
+404                if (len(t) < 4):
+405                    break
+406                nc = struct.unpack('i', t)[0]
+407                configlist[-1].append(nc)
+408
+409                t = fp.read(8 * tmax * (nn + 1))
+410                if kwargs.get('plaquette'):
+411                    if nc % dtr_read == 0:
+412                        Ysl.append(struct.unpack('d' * tmax * (nn + 1), t))
+413                t = fp.read(8 * tmax * (nn + 1))
+414                if not kwargs.get('plaquette'):
+415                    if nc % dtr_read == 0:
+416                        Ysl.append(struct.unpack('d' * tmax * (nn + 1), t))
+417                t = fp.read(8 * tmax * (nn + 1))
+418
+419        Ysum.append([])
+420        for i, item in enumerate(Ysl):
+421            Ysum[-1].append([np.mean(item[current + xmin:
+422                             current + tmax - xmin])
+423                            for current in range(0, len(item), tmax)])
+424
+425        diffmeas = configlist[-1][-1] - configlist[-1][-2]
+426        configlist[-1] = [item // diffmeas for item in configlist[-1]]
+427        if kwargs.get('assume_thermalization', True) and configlist[-1][0] > 1:
+428            warnings.warn('Assume thermalization and that the first measurement belongs to the first config.')
+429            offset = configlist[-1][0] - 1
+430            configlist[-1] = [item - offset for item in configlist[-1]]
 431
-432    idl = [range(configlist[rep][r_start_index[rep]], configlist[rep][r_stop_index[rep]] + 1, r_step) for rep in range(replica)]
-433    t2E_dict = {}
-434    for n in range(nn + 1):
-435        samples = []
-436        for nrep, rep in enumerate(Ysum):
-437            samples.append([])
-438            for cnfg in rep:
-439                samples[-1].append(cnfg[n])
-440            samples[-1] = samples[-1][r_start_index[nrep]:r_stop_index[nrep] + 1][::r_step]
-441        new_obs = Obs(samples, rep_names, idl=idl)
-442        t2E_dict[n * dn * eps] = (n * dn * eps) ** 2 * new_obs / (spatial_extent ** 3) - 0.3
-443
-444    zero_crossing = np.argmax(np.array(
-445        [o.value for o in t2E_dict.values()]) > 0.0)
-446
-447    x = list(t2E_dict.keys())[zero_crossing - fit_range:
-448                              zero_crossing + fit_range]
-449    y = list(t2E_dict.values())[zero_crossing - fit_range:
-450                                zero_crossing + fit_range]
-451    [o.gamma_method() for o in y]
-452
-453    fit_result = fit_lin(x, y)
-454
-455    if kwargs.get('plot_fit'):
-456        plt.figure()
-457        gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1], wspace=0.0, hspace=0.0)
-458        ax0 = plt.subplot(gs[0])
-459        xmore = list(t2E_dict.keys())[zero_crossing - fit_range - 2: zero_crossing + fit_range + 2]
-460        ymore = list(t2E_dict.values())[zero_crossing - fit_range - 2: zero_crossing + fit_range + 2]
-461        [o.gamma_method() for o in ymore]
-462        ax0.errorbar(xmore, [yi.value for yi in ymore], yerr=[yi.dvalue for yi in ymore], fmt='x')
-463        xplot = np.linspace(np.min(x), np.max(x))
-464        yplot = [fit_result[0] + fit_result[1] * xi for xi in xplot]
-465        [yi.gamma_method() for yi in yplot]
-466        ax0.fill_between(xplot, y1=[yi.value - yi.dvalue for yi in yplot], y2=[yi.value + yi.dvalue for yi in yplot])
-467        retval = (-fit_result[0] / fit_result[1])
-468        retval.gamma_method()
-469        ylim = ax0.get_ylim()
-470        ax0.fill_betweenx(ylim, x1=retval.value - retval.dvalue, x2=retval.value + retval.dvalue, color='gray', alpha=0.4)
-471        ax0.set_ylim(ylim)
-472        ax0.set_ylabel(r'$t^2 \langle E(t) \rangle - 0.3 $')
-473        xlim = ax0.get_xlim()
-474
-475        fit_res = [fit_result[0] + fit_result[1] * xi for xi in x]
-476        residuals = (np.asarray([o.value for o in y]) - [o.value for o in fit_res]) / np.asarray([o.dvalue for o in y])
-477        ax1 = plt.subplot(gs[1])
-478        ax1.plot(x, residuals, 'ko', ls='none', markersize=5)
-479        ax1.tick_params(direction='out')
-480        ax1.tick_params(axis="x", bottom=True, top=True, labelbottom=True)
-481        ax1.axhline(y=0.0, ls='--', color='k')
-482        ax1.fill_between(xlim, -1.0, 1.0, alpha=0.1, facecolor='k')
-483        ax1.set_xlim(xlim)
-484        ax1.set_ylabel('Residuals')
-485        ax1.set_xlabel(r'$t/a^2$')
-486
-487        plt.draw()
-488    return -fit_result[0] / fit_result[1]
+432        if r_start[rep] is None:
+433            r_start_index.append(0)
+434        else:
+435            try:
+436                r_start_index.append(configlist[-1].index(r_start[rep]))
+437            except ValueError:
+438                raise Exception('Config %d not in file with range [%d, %d]' % (
+439                    r_start[rep], configlist[-1][0], configlist[-1][-1])) from None
+440
+441        if r_stop[rep] is None:
+442            r_stop_index.append(len(configlist[-1]) - 1)
+443        else:
+444            try:
+445                r_stop_index.append(configlist[-1].index(r_stop[rep]))
+446            except ValueError:
+447                raise Exception('Config %d not in file with range [%d, %d]' % (
+448                    r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None
+449
+450    if np.any([len(np.unique(np.diff(cl))) != 1 for cl in configlist]):
+451        raise Exception('Irregular spaced data in input file!', [len(np.unique(np.diff(cl))) for cl in configlist])
+452    stepsizes = [list(np.unique(np.diff(cl)))[0] for cl in configlist]
+453    if np.any([step != 1 for step in stepsizes]):
+454        warnings.warn('Stepsize between configurations is greater than one!' + str(stepsizes), RuntimeWarning)
+455
+456    idl = [range(configlist[rep][r_start_index[rep]], configlist[rep][r_stop_index[rep]] + 1, r_step) for rep in range(replica)]
+457    t2E_dict = {}
+458    for n in range(nn + 1):
+459        samples = []
+460        for nrep, rep in enumerate(Ysum):
+461            samples.append([])
+462            for cnfg in rep:
+463                samples[-1].append(cnfg[n])
+464            samples[-1] = samples[-1][r_start_index[nrep]:r_stop_index[nrep] + 1][::r_step]
+465        new_obs = Obs(samples, rep_names, idl=idl)
+466        t2E_dict[n * dn * eps] = (n * dn * eps) ** 2 * new_obs / (spatial_extent ** 3) - 0.3
+467
+468    zero_crossing = np.argmax(np.array(
+469        [o.value for o in t2E_dict.values()]) > 0.0)
+470
+471    x = list(t2E_dict.keys())[zero_crossing - fit_range:
+472                              zero_crossing + fit_range]
+473    y = list(t2E_dict.values())[zero_crossing - fit_range:
+474                                zero_crossing + fit_range]
+475    [o.gamma_method() for o in y]
+476
+477    fit_result = fit_lin(x, y)
+478
+479    if kwargs.get('plot_fit'):
+480        plt.figure()
+481        gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1], wspace=0.0, hspace=0.0)
+482        ax0 = plt.subplot(gs[0])
+483        xmore = list(t2E_dict.keys())[zero_crossing - fit_range - 2: zero_crossing + fit_range + 2]
+484        ymore = list(t2E_dict.values())[zero_crossing - fit_range - 2: zero_crossing + fit_range + 2]
+485        [o.gamma_method() for o in ymore]
+486        ax0.errorbar(xmore, [yi.value for yi in ymore], yerr=[yi.dvalue for yi in ymore], fmt='x')
+487        xplot = np.linspace(np.min(x), np.max(x))
+488        yplot = [fit_result[0] + fit_result[1] * xi for xi in xplot]
+489        [yi.gamma_method() for yi in yplot]
+490        ax0.fill_between(xplot, y1=[yi.value - yi.dvalue for yi in yplot], y2=[yi.value + yi.dvalue for yi in yplot])
+491        retval = (-fit_result[0] / fit_result[1])
+492        retval.gamma_method()
+493        ylim = ax0.get_ylim()
+494        ax0.fill_betweenx(ylim, x1=retval.value - retval.dvalue, x2=retval.value + retval.dvalue, color='gray', alpha=0.4)
+495        ax0.set_ylim(ylim)
+496        ax0.set_ylabel(r'$t^2 \langle E(t) \rangle - 0.3 $')
+497        xlim = ax0.get_xlim()
+498
+499        fit_res = [fit_result[0] + fit_result[1] * xi for xi in x]
+500        residuals = (np.asarray([o.value for o in y]) - [o.value for o in fit_res]) / np.asarray([o.dvalue for o in y])
+501        ax1 = plt.subplot(gs[1])
+502        ax1.plot(x, residuals, 'ko', ls='none', markersize=5)
+503        ax1.tick_params(direction='out')
+504        ax1.tick_params(axis="x", bottom=True, top=True, labelbottom=True)
+505        ax1.axhline(y=0.0, ls='--', color='k')
+506        ax1.fill_between(xlim, -1.0, 1.0, alpha=0.1, facecolor='k')
+507        ax1.set_xlim(xlim)
+508        ax1.set_ylabel('Residuals')
+509        ax1.set_xlabel(r'$t/a^2$')
+510
+511        plt.draw()
+512    return -fit_result[0] / fit_result[1]
 

 
 
@@ -1883,57 +1889,57 @@ Extracted t0
- 
536def read_qtop(path, prefix, c, dtr_cnfg=1, version="openQCD", **kwargs):
-537    """Read the topologial charge based on openQCD gradient flow measurements.
-538
-539    Parameters
-540    ----------
-541    path : str
-542        path of the measurement files
-543    prefix : str
-544        prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat.
-545        Ignored if file names are passed explicitly via keyword files.
-546    c : double
-547        Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L.
-548    dtr_cnfg : int
-549        (optional) parameter that specifies the number of measurements
-550        between two configs.
-551        If it is not set, the distance between two measurements
-552        in the file is assumed to be the distance between two configurations.
-553    steps : int
-554        (optional) Distance between two configurations in units of trajectories /
-555         cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
-556    version : str
-557        Either openQCD or sfqcd, depending on the data.
-558    L : int
-559        spatial length of the lattice in L/a.
-560        HAS to be set if version != sfqcd, since openQCD does not provide
-561        this in the header
-562    r_start : list
-563        list which contains the first config to be read for each replicum.
-564    r_stop : list
-565        list which contains the last config to be read for each replicum.
-566    files : list
-567        specify the exact files that need to be read
-568        from path, practical if e.g. only one replicum is needed
-569    postfix : str
-570        postfix of the file to read, e.g. '.gfms.dat' for openQCD-files
-571    names : list
-572        Alternative labeling for replicas/ensembles.
-573        Has to have the appropriate length.
-574    Zeuthen_flow : bool
-575        (optional) If True, the Zeuthen flow is used for Qtop. Only possible
-576        for version=='sfqcd' If False, the Wilson flow is used.
-577    integer_charge : bool
-578        If True, the charge is rounded towards the nearest integer on each config.
-579
-580    Returns
-581    -------
-582    result : Obs
-583        Read topological charge
-584    """
-585
-586    return _read_flow_obs(path, prefix, c, dtr_cnfg=dtr_cnfg, version=version, obspos=0, **kwargs)
+            
560def read_qtop(path, prefix, c, dtr_cnfg=1, version="openQCD", **kwargs):
+561    """Read the topologial charge based on openQCD gradient flow measurements.
+562
+563    Parameters
+564    ----------
+565    path : str
+566        path of the measurement files
+567    prefix : str
+568        prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat.
+569        Ignored if file names are passed explicitly via keyword files.
+570    c : double
+571        Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L.
+572    dtr_cnfg : int
+573        (optional) parameter that specifies the number of measurements
+574        between two configs.
+575        If it is not set, the distance between two measurements
+576        in the file is assumed to be the distance between two configurations.
+577    steps : int
+578        (optional) Distance between two configurations in units of trajectories /
+579         cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
+580    version : str
+581        Either openQCD or sfqcd, depending on the data.
+582    L : int
+583        spatial length of the lattice in L/a.
+584        HAS to be set if version != sfqcd, since openQCD does not provide
+585        this in the header
+586    r_start : list
+587        list which contains the first config to be read for each replicum.
+588    r_stop : list
+589        list which contains the last config to be read for each replicum.
+590    files : list
+591        specify the exact files that need to be read
+592        from path, practical if e.g. only one replicum is needed
+593    postfix : str
+594        postfix of the file to read, e.g. '.gfms.dat' for openQCD-files
+595    names : list
+596        Alternative labeling for replicas/ensembles.
+597        Has to have the appropriate length.
+598    Zeuthen_flow : bool
+599        (optional) If True, the Zeuthen flow is used for Qtop. Only possible
+600        for version=='sfqcd' If False, the Wilson flow is used.
+601    integer_charge : bool
+602        If True, the charge is rounded towards the nearest integer on each config.
+603
+604    Returns
+605    -------
+606    result : Obs
+607        Read topological charge
+608    """
+609
+610    return _read_flow_obs(path, prefix, c, dtr_cnfg=dtr_cnfg, version=version, obspos=0, **kwargs)
 

 
 
@@ -2003,76 +2009,76 @@ Read topological charge
- 
589def read_gf_coupling(path, prefix, c, dtr_cnfg=1, Zeuthen_flow=True, **kwargs):
-590    """Read the gradient flow coupling based on sfqcd gradient flow measurements. See 1607.06423 for details.
-591
-592    Note: The current implementation only works for c=0.3 and T=L. The definition of the coupling in 1607.06423 requires projection to topological charge zero which is not done within this function but has to be performed in a separate step.
-593
-594    Parameters
-595    ----------
-596    path : str
-597        path of the measurement files
-598    prefix : str
-599        prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat.
-600        Ignored if file names are passed explicitly via keyword files.
-601    c : double
-602        Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L.
-603    dtr_cnfg : int
-604        (optional) parameter that specifies the number of measurements
-605        between two configs.
-606        If it is not set, the distance between two measurements
-607        in the file is assumed to be the distance between two configurations.
-608    steps : int
-609        (optional) Distance between two configurations in units of trajectories /
-610         cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
-611    r_start : list
-612        list which contains the first config to be read for each replicum.
-613    r_stop : list
-614        list which contains the last config to be read for each replicum.
-615    files : list
-616        specify the exact files that need to be read
-617        from path, practical if e.g. only one replicum is needed
-618    names : list
-619        Alternative labeling for replicas/ensembles.
-620        Has to have the appropriate length.
-621    postfix : str
-622        postfix of the file to read, e.g. '.gfms.dat' for openQCD-files
-623    Zeuthen_flow : bool
-624        (optional) If True, the Zeuthen flow is used for the coupling. If False, the Wilson flow is used.
-625    """
-626
-627    if c != 0.3:
-628        raise Exception("The required lattice norm is only implemented for c=0.3 at the moment.")
-629
-630    plaq = _read_flow_obs(path, prefix, c, dtr_cnfg=dtr_cnfg, version="sfqcd", obspos=6, sum_t=False, Zeuthen_flow=Zeuthen_flow, integer_charge=False, **kwargs)
-631    C2x1 = _read_flow_obs(path, prefix, c, dtr_cnfg=dtr_cnfg, version="sfqcd", obspos=7, sum_t=False, Zeuthen_flow=Zeuthen_flow, integer_charge=False, **kwargs)
-632    L = plaq.tag["L"]
-633    T = plaq.tag["T"]
-634
-635    if T != L:
-636        raise Exception("The required lattice norm is only implemented for T=L at the moment.")
-637
-638    if Zeuthen_flow is not True:
-639        raise Exception("The required lattice norm is only implemented for the Zeuthen flow at the moment.")
-640
-641    t = (c * L) ** 2 / 8
-642
-643    normdict = {4: 0.012341170468270,
-644                6: 0.010162691462430,
-645                8: 0.009031614807931,
-646                10: 0.008744966371393,
-647                12: 0.008650917856809,
-648                14: 8.611154391267955E-03,
-649                16: 0.008591758449508,
-650                20: 0.008575359627103,
-651                24: 0.008569387847540,
-652                28: 8.566803713382559E-03,
-653                32: 0.008565541650006,
-654                40: 8.564480684962046E-03,
-655                48: 8.564098025073460E-03,
-656                64: 8.563853943383087E-03}
-657
-658    return t * t * (5 / 3 * plaq - 1 / 12 * C2x1) / normdict[L]
+            
613def read_gf_coupling(path, prefix, c, dtr_cnfg=1, Zeuthen_flow=True, **kwargs):
+614    """Read the gradient flow coupling based on sfqcd gradient flow measurements. See 1607.06423 for details.
+615
+616    Note: The current implementation only works for c=0.3 and T=L. The definition of the coupling in 1607.06423 requires projection to topological charge zero which is not done within this function but has to be performed in a separate step.
+617
+618    Parameters
+619    ----------
+620    path : str
+621        path of the measurement files
+622    prefix : str
+623        prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat.
+624        Ignored if file names are passed explicitly via keyword files.
+625    c : double
+626        Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L.
+627    dtr_cnfg : int
+628        (optional) parameter that specifies the number of measurements
+629        between two configs.
+630        If it is not set, the distance between two measurements
+631        in the file is assumed to be the distance between two configurations.
+632    steps : int
+633        (optional) Distance between two configurations in units of trajectories /
+634         cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
+635    r_start : list
+636        list which contains the first config to be read for each replicum.
+637    r_stop : list
+638        list which contains the last config to be read for each replicum.
+639    files : list
+640        specify the exact files that need to be read
+641        from path, practical if e.g. only one replicum is needed
+642    names : list
+643        Alternative labeling for replicas/ensembles.
+644        Has to have the appropriate length.
+645    postfix : str
+646        postfix of the file to read, e.g. '.gfms.dat' for openQCD-files
+647    Zeuthen_flow : bool
+648        (optional) If True, the Zeuthen flow is used for the coupling. If False, the Wilson flow is used.
+649    """
+650
+651    if c != 0.3:
+652        raise Exception("The required lattice norm is only implemented for c=0.3 at the moment.")
+653
+654    plaq = _read_flow_obs(path, prefix, c, dtr_cnfg=dtr_cnfg, version="sfqcd", obspos=6, sum_t=False, Zeuthen_flow=Zeuthen_flow, integer_charge=False, **kwargs)
+655    C2x1 = _read_flow_obs(path, prefix, c, dtr_cnfg=dtr_cnfg, version="sfqcd", obspos=7, sum_t=False, Zeuthen_flow=Zeuthen_flow, integer_charge=False, **kwargs)
+656    L = plaq.tag["L"]
+657    T = plaq.tag["T"]
+658
+659    if T != L:
+660        raise Exception("The required lattice norm is only implemented for T=L at the moment.")
+661
+662    if Zeuthen_flow is not True:
+663        raise Exception("The required lattice norm is only implemented for the Zeuthen flow at the moment.")
+664
+665    t = (c * L) ** 2 / 8
+666
+667    normdict = {4: 0.012341170468270,
+668                6: 0.010162691462430,
+669                8: 0.009031614807931,
+670                10: 0.008744966371393,
+671                12: 0.008650917856809,
+672                14: 8.611154391267955E-03,
+673                16: 0.008591758449508,
+674                20: 0.008575359627103,
+675                24: 0.008569387847540,
+676                28: 8.566803713382559E-03,
+677                32: 0.008565541650006,
+678                40: 8.564480684962046E-03,
+679                48: 8.564098025073460E-03,
+680                64: 8.563853943383087E-03}
+681
+682    return t * t * (5 / 3 * plaq - 1 / 12 * C2x1) / normdict[L]
 

 
 
@@ -2128,30 +2134,30 @@ postfix of the file to read, e.g. '.gfms.dat' for openQCD-files
- 
938def qtop_projection(qtop, target=0):
-939    """Returns the projection to the topological charge sector defined by target.
-940
-941    Parameters
-942    ----------
-943    path : Obs
-944        Topological charge.
-945    target : int
-946        Specifies the topological sector to be reweighted to (default 0)
-947
-948    Returns
-949    -------
-950    reto : Obs
-951        projection to the topological charge sector defined by target
-952    """
-953    if qtop.reweighted:
-954        raise Exception('You can not use a reweighted observable for reweighting!')
-955
-956    proj_qtop = []
-957    for n in qtop.deltas:
-958        proj_qtop.append(np.array([1 if round(qtop.r_values[n] + q) == target else 0 for q in qtop.deltas[n]]))
-959
-960    reto = Obs(proj_qtop, qtop.names, idl=[qtop.idl[name] for name in qtop.names])
-961    return reto
+            
954def qtop_projection(qtop, target=0):
+955    """Returns the projection to the topological charge sector defined by target.
+956
+957    Parameters
+958    ----------
+959    path : Obs
+960        Topological charge.
+961    target : int
+962        Specifies the topological sector to be reweighted to (default 0)
+963
+964    Returns
+965    -------
+966    reto : Obs
+967        projection to the topological charge sector defined by target
+968    """
+969    if qtop.reweighted:
+970        raise Exception('You can not use a reweighted observable for reweighting!')
+971
+972    proj_qtop = []
+973    for n in qtop.deltas:
+974        proj_qtop.append(np.array([1 if round(qtop.r_values[n] + q) == target else 0 for q in qtop.deltas[n]]))
+975
+976    reto = Obs(proj_qtop, qtop.names, idl=[qtop.idl[name] for name in qtop.names])
+977    return reto
 

 
 
@@ -2187,62 +2193,62 @@ projection to the topological charge sector defined by target
- 
 964def read_qtop_sector(path, prefix, c, target=0, **kwargs):
- 965    """Constructs reweighting factors to a specified topological sector.
- 966
- 967    Parameters
- 968    ----------
- 969    path : str
- 970        path of the measurement files
- 971    prefix : str
- 972        prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat
- 973    c : double
- 974        Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L
- 975    target : int
- 976        Specifies the topological sector to be reweighted to (default 0)
- 977    dtr_cnfg : int
- 978        (optional) parameter that specifies the number of trajectories
- 979        between two configs.
- 980        if it is not set, the distance between two measurements
- 981        in the file is assumed to be the distance between two configurations.
- 982    steps : int
- 983        (optional) Distance between two configurations in units of trajectories /
- 984         cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
- 985    version : str
- 986        version string of the openQCD (sfqcd) version used to create
- 987        the ensemble. Default is 2.0. May also be set to sfqcd.
- 988    L : int
- 989        spatial length of the lattice in L/a.
- 990        HAS to be set if version != sfqcd, since openQCD does not provide
- 991        this in the header
- 992    r_start : list
- 993        offset of the first ensemble, making it easier to match
- 994        later on with other Obs
- 995    r_stop : list
- 996        last configurations that need to be read (per replicum)
- 997    files : list
- 998        specify the exact files that need to be read
- 999        from path, practical if e.g. only one replicum is needed
-1000    names : list
-1001        Alternative labeling for replicas/ensembles.
-1002        Has to have the appropriate length
-1003    Zeuthen_flow : bool
-1004        (optional) If True, the Zeuthen flow is used for Qtop. Only possible
-1005        for version=='sfqcd' If False, the Wilson flow is used.
-1006
-1007    Returns
-1008    -------
-1009    reto : Obs
-1010        projection to the topological charge sector defined by target
-1011    """
-1012
-1013    if not isinstance(target, int):
-1014        raise Exception("'target' has to be an integer.")
-1015
-1016    kwargs['integer_charge'] = True
-1017    qtop = read_qtop(path, prefix, c, **kwargs)
-1018
-1019    return qtop_projection(qtop, target=target)
+            
 980def read_qtop_sector(path, prefix, c, target=0, **kwargs):
+ 981    """Constructs reweighting factors to a specified topological sector.
+ 982
+ 983    Parameters
+ 984    ----------
+ 985    path : str
+ 986        path of the measurement files
+ 987    prefix : str
+ 988        prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat
+ 989    c : double
+ 990        Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L
+ 991    target : int
+ 992        Specifies the topological sector to be reweighted to (default 0)
+ 993    dtr_cnfg : int
+ 994        (optional) parameter that specifies the number of trajectories
+ 995        between two configs.
+ 996        if it is not set, the distance between two measurements
+ 997        in the file is assumed to be the distance between two configurations.
+ 998    steps : int
+ 999        (optional) Distance between two configurations in units of trajectories /
+1000         cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
+1001    version : str
+1002        version string of the openQCD (sfqcd) version used to create
+1003        the ensemble. Default is 2.0. May also be set to sfqcd.
+1004    L : int
+1005        spatial length of the lattice in L/a.
+1006        HAS to be set if version != sfqcd, since openQCD does not provide
+1007        this in the header
+1008    r_start : list
+1009        offset of the first ensemble, making it easier to match
+1010        later on with other Obs
+1011    r_stop : list
+1012        last configurations that need to be read (per replicum)
+1013    files : list
+1014        specify the exact files that need to be read
+1015        from path, practical if e.g. only one replicum is needed
+1016    names : list
+1017        Alternative labeling for replicas/ensembles.
+1018        Has to have the appropriate length
+1019    Zeuthen_flow : bool
+1020        (optional) If True, the Zeuthen flow is used for Qtop. Only possible
+1021        for version=='sfqcd' If False, the Wilson flow is used.
+1022
+1023    Returns
+1024    -------
+1025    reto : Obs
+1026        projection to the topological charge sector defined by target
+1027    """
+1028
+1029    if not isinstance(target, int):
+1030        raise Exception("'target' has to be an integer.")
+1031
+1032    kwargs['integer_charge'] = True
+1033    qtop = read_qtop(path, prefix, c, **kwargs)
+1034
+1035    return qtop_projection(qtop, target=target)
 

 
 
@@ -2311,153 +2317,159 @@ projection to the topological charge sector defined by target
- 
1022def read_ms5_xsf(path, prefix, qc, corr, sep="r", **kwargs):
-1023    """
-1024    Read data from files in the specified directory with the specified prefix and quark combination extension, and return a `Corr` object containing the data.
-1025
-1026    Parameters
-1027    ----------
-1028    path : str
-1029        The directory to search for the files in.
-1030    prefix : str
-1031        The prefix to match the files against.
-1032    qc : str
-1033        The quark combination extension to match the files against.
-1034    corr : str
-1035        The correlator to extract data for.
-1036    sep : str, optional
-1037        The separator to use when parsing the replika names.
-1038    **kwargs
-1039        Additional keyword arguments. The following keyword arguments are recognized:
-1040
-1041        - names (List[str]): A list of names to use for the replicas.
-1042
-1043    Returns
-1044    -------
-1045    Corr
-1046        A complex valued `Corr` object containing the data read from the files. In case of boudary to bulk correlators.
-1047    or
-1048    CObs
-1049        A complex valued `CObs` object containing the data read from the files. In case of boudary to boundary correlators.
-1050
-1051
-1052    Raises
-1053    ------
-1054    FileNotFoundError
-1055        If no files matching the specified prefix and quark combination extension are found in the specified directory.
-1056    IOError
-1057        If there is an error reading a file.
-1058    struct.error
-1059        If there is an error unpacking binary data.
-1060    """
-1061
-1062    found = []
-1063    files = []
-1064    names = []
-1065
-1066    if "names" in kwargs:
-1067        names = kwargs.get("names")
-1068
-1069    for (dirpath, dirnames, filenames) in os.walk(path + "/"):
-1070        found.extend(filenames)
-1071        break
-1072
-1073    for f in found:
-1074        if fnmatch.fnmatch(f, prefix + "*.ms5_xsf_" + qc + ".dat"):
-1075            files.append(f)
-1076            if "names" not in kwargs:
-1077                if not sep == "":
-1078                    se = f.split(".")[0]
-1079                    for s in f.split(".")[1:-1]:
-1080                        se += "." + s
-1081                    names.append(se.split(sep)[0] + "|r" + se.split(sep)[1])
-1082                else:
-1083                    names.append(prefix)
-1084
-1085    names = sorted(names)
-1086    files = sorted(files)
-1087
-1088    cnfgs = []
-1089    realsamples = []
-1090    imagsamples = []
-1091    repnum = 0
-1092    for file in files:
-1093        with open(path + "/" + file, "rb") as fp:
+            
1038def read_ms5_xsf(path, prefix, qc, corr, sep="r", **kwargs):
+1039    """
+1040    Read data from files in the specified directory with the specified prefix and quark combination extension, and return a `Corr` object containing the data.
+1041
+1042    Parameters
+1043    ----------
+1044    path : str
+1045        The directory to search for the files in.
+1046    prefix : str
+1047        The prefix to match the files against.
+1048    qc : str
+1049        The quark combination extension to match the files against.
+1050    corr : str
+1051        The correlator to extract data for.
+1052    sep : str, optional
+1053        The separator to use when parsing the replika names.
+1054    **kwargs
+1055        Additional keyword arguments. The following keyword arguments are recognized:
+1056
+1057        - names (List[str]): A list of names to use for the replicas.
+1058
+1059    Returns
+1060    -------
+1061    Corr
+1062        A complex valued `Corr` object containing the data read from the files. In case of boudary to bulk correlators.
+1063    or
+1064    CObs
+1065        A complex valued `CObs` object containing the data read from the files. In case of boudary to boundary correlators.
+1066
+1067
+1068    Raises
+1069    ------
+1070    FileNotFoundError
+1071        If no files matching the specified prefix and quark combination extension are found in the specified directory.
+1072    IOError
+1073        If there is an error reading a file.
+1074    struct.error
+1075        If there is an error unpacking binary data.
+1076    """
+1077
+1078    # found = []
+1079    files = []
+1080    names = []
+1081
+1082    # test if the input is correct
+1083    if qc not in ['dd', 'ud', 'du', 'uu']:
+1084        raise Exception("Unknown quark conbination!")
+1085
+1086    if corr not in ["gS", "gP", "gA", "gV", "gVt", "lA", "lV", "lVt", "lT", "lTt", "g1", "l1"]:
+1087        raise Exception("Unknown correlator!")
+1088
+1089    if "files" in kwargs:
+1090        known_files = kwargs.get("files")
+1091    else:
+1092        known_files = []
+1093    files = _find_files(path, prefix, "ms5_xsf_" + qc, "dat", known_files=known_files)
 1094
-1095            t = fp.read(8)
-1096            kappa = struct.unpack('d', t)[0]
-1097            t = fp.read(8)
-1098            csw = struct.unpack('d', t)[0]
-1099            t = fp.read(8)
-1100            dF = struct.unpack('d', t)[0]
-1101            t = fp.read(8)
-1102            zF = struct.unpack('d', t)[0]
-1103
-1104            t = fp.read(4)
-1105            tmax = struct.unpack('i', t)[0]
-1106            t = fp.read(4)
-1107            bnd = struct.unpack('i', t)[0]
-1108
-1109            placesBI = ["gS", "gP",
-1110                        "gA", "gV",
-1111                        "gVt", "lA",
-1112                        "lV", "lVt",
-1113                        "lT", "lTt"]
-1114            placesBB = ["g1", "l1"]
-1115
-1116            # the chunks have the following structure:
-1117            # confignumber, 10x timedependent complex correlators as doubles, 2x timeindependent complex correlators as doubles
-1118
-1119            chunksize = 4 + (8 * 2 * tmax * 10) + (8 * 2 * 2)
-1120            packstr = '=i' + ('d' * 2 * tmax * 10) + ('d' * 2 * 2)
-1121            cnfgs.append([])
-1122            realsamples.append([])
-1123            imagsamples.append([])
-1124            for t in range(tmax):
-1125                realsamples[repnum].append([])
-1126                imagsamples[repnum].append([])
-1127
-1128            while True:
-1129                cnfgt = fp.read(chunksize)
-1130                if not cnfgt:
-1131                    break
-1132                asascii = struct.unpack(packstr, cnfgt)
-1133                cnfg = asascii[0]
-1134                cnfgs[repnum].append(cnfg)
-1135
-1136                if corr not in placesBB:
-1137                    tmpcorr = asascii[1 + 2 * tmax * placesBI.index(corr):1 + 2 * tmax * placesBI.index(corr) + 2 * tmax]
-1138                else:
-1139                    tmpcorr = asascii[1 + 2 * tmax * len(placesBI) + 2 * placesBB.index(corr):1 + 2 * tmax * len(placesBI) + 2 * placesBB.index(corr) + 2]
+1095    if "names" in kwargs:
+1096        names = kwargs.get("names")
+1097    else:
+1098        for f in files:
+1099            if not sep == "":
+1100                se = f.split(".")[0]
+1101                for s in f.split(".")[1:-2]:
+1102                    se += "." + s
+1103                names.append(se.split(sep)[0] + "|r" + se.split(sep)[1])
+1104            else:
+1105                names.append(prefix)
+1106
+1107    names = sorted(names)
+1108    files = sorted(files)
+1109
+1110    cnfgs = []
+1111    realsamples = []
+1112    imagsamples = []
+1113    repnum = 0
+1114    for file in files:
+1115        with open(path + "/" + file, "rb") as fp:
+1116
+1117            t = fp.read(8)
+1118            kappa = struct.unpack('d', t)[0]
+1119            t = fp.read(8)
+1120            csw = struct.unpack('d', t)[0]
+1121            t = fp.read(8)
+1122            dF = struct.unpack('d', t)[0]
+1123            t = fp.read(8)
+1124            zF = struct.unpack('d', t)[0]
+1125
+1126            t = fp.read(4)
+1127            tmax = struct.unpack('i', t)[0]
+1128            t = fp.read(4)
+1129            bnd = struct.unpack('i', t)[0]
+1130
+1131            placesBI = ["gS", "gP",
+1132                        "gA", "gV",
+1133                        "gVt", "lA",
+1134                        "lV", "lVt",
+1135                        "lT", "lTt"]
+1136            placesBB = ["g1", "l1"]
+1137
+1138            # the chunks have the following structure:
+1139            # confignumber, 10x timedependent complex correlators as doubles, 2x timeindependent complex correlators as doubles
 1140
-1141                corrres = [[], []]
-1142                for i in range(len(tmpcorr)):
-1143                    corrres[i % 2].append(tmpcorr[i])
-1144                for t in range(int(len(tmpcorr) / 2)):
-1145                    realsamples[repnum][t].append(corrres[0][t])
-1146                for t in range(int(len(tmpcorr) / 2)):
-1147                    imagsamples[repnum][t].append(corrres[1][t])
-1148        repnum += 1
+1141            chunksize = 4 + (8 * 2 * tmax * 10) + (8 * 2 * 2)
+1142            packstr = '=i' + ('d' * 2 * tmax * 10) + ('d' * 2 * 2)
+1143            cnfgs.append([])
+1144            realsamples.append([])
+1145            imagsamples.append([])
+1146            for t in range(tmax):
+1147                realsamples[repnum].append([])
+1148                imagsamples[repnum].append([])
 1149
-1150    s = "Read correlator " + corr + " from " + str(repnum) + " replika with " + str(len(realsamples[0][t]))
-1151    for rep in range(1, repnum):
-1152        s += ", " + str(len(realsamples[rep][t]))
-1153    s += " samples"
-1154    print(s)
-1155    print("Asserted run parameters:\n T:", tmax, "kappa:", kappa, "csw:", csw, "dF:", dF, "zF:", zF, "bnd:", bnd)
-1156
-1157    # we have the data now... but we need to re format the whole thing and put it into Corr objects.
-1158
-1159    compObs = []
-1160
-1161    for t in range(int(len(tmpcorr) / 2)):
-1162        compObs.append(CObs(Obs([realsamples[rep][t] for rep in range(repnum)], names=names, idl=cnfgs),
-1163                            Obs([imagsamples[rep][t] for rep in range(repnum)], names=names, idl=cnfgs)))
-1164
-1165    if len(compObs) == 1:
-1166        return compObs[0]
-1167    else:
-1168        return Corr(compObs)
+1150            while True:
+1151                cnfgt = fp.read(chunksize)
+1152                if not cnfgt:
+1153                    break
+1154                asascii = struct.unpack(packstr, cnfgt)
+1155                cnfg = asascii[0]
+1156                cnfgs[repnum].append(cnfg)
+1157
+1158                if corr not in placesBB:
+1159                    tmpcorr = asascii[1 + 2 * tmax * placesBI.index(corr):1 + 2 * tmax * placesBI.index(corr) + 2 * tmax]
+1160                else:
+1161                    tmpcorr = asascii[1 + 2 * tmax * len(placesBI) + 2 * placesBB.index(corr):1 + 2 * tmax * len(placesBI) + 2 * placesBB.index(corr) + 2]
+1162
+1163                corrres = [[], []]
+1164                for i in range(len(tmpcorr)):
+1165                    corrres[i % 2].append(tmpcorr[i])
+1166                for t in range(int(len(tmpcorr) / 2)):
+1167                    realsamples[repnum][t].append(corrres[0][t])
+1168                for t in range(int(len(tmpcorr) / 2)):
+1169                    imagsamples[repnum][t].append(corrres[1][t])
+1170        repnum += 1
+1171
+1172    s = "Read correlator " + corr + " from " + str(repnum) + " replika with " + str(len(realsamples[0][t]))
+1173    for rep in range(1, repnum):
+1174        s += ", " + str(len(realsamples[rep][t]))
+1175    s += " samples"
+1176    print(s)
+1177    print("Asserted run parameters:\n T:", tmax, "kappa:", kappa, "csw:", csw, "dF:", dF, "zF:", zF, "bnd:", bnd)
+1178
+1179    # we have the data now... but we need to re format the whole thing and put it into Corr objects.
+1180
+1181    compObs = []
+1182
+1183    for t in range(int(len(tmpcorr) / 2)):
+1184        compObs.append(CObs(Obs([realsamples[rep][t] for rep in range(repnum)], names=names, idl=cnfgs),
+1185                            Obs([imagsamples[rep][t] for rep in range(repnum)], names=names, idl=cnfgs)))
+1186
+1187    if len(compObs) == 1:
+1188        return compObs[0]
+1189    else:
+1190        return Corr(compObs)