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

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

 
 
@@ -1889,57 +1925,57 @@ Extracted t0
- 
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)
+            
591def read_qtop(path, prefix, c, dtr_cnfg=1, version="openQCD", **kwargs):
+592    """Read the topologial charge based on openQCD gradient flow measurements.
+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    version : str
+612        Either openQCD or sfqcd, depending on the data.
+613    L : int
+614        spatial length of the lattice in L/a.
+615        HAS to be set if version != sfqcd, since openQCD does not provide
+616        this in the header
+617    r_start : list
+618        list which contains the first config to be read for each replicum.
+619    r_stop : list
+620        list which contains the last config to be read for each replicum.
+621    files : list
+622        specify the exact files that need to be read
+623        from path, practical if e.g. only one replicum is needed
+624    postfix : str
+625        postfix of the file to read, e.g. '.gfms.dat' for openQCD-files
+626    names : list
+627        Alternative labeling for replicas/ensembles.
+628        Has to have the appropriate length.
+629    Zeuthen_flow : bool
+630        (optional) If True, the Zeuthen flow is used for Qtop. Only possible
+631        for version=='sfqcd' If False, the Wilson flow is used.
+632    integer_charge : bool
+633        If True, the charge is rounded towards the nearest integer on each config.
+634
+635    Returns
+636    -------
+637    result : Obs
+638        Read topological charge
+639    """
+640
+641    return _read_flow_obs(path, prefix, c, dtr_cnfg=dtr_cnfg, version=version, obspos=0, **kwargs)
 

 
 
@@ -2009,76 +2045,76 @@ Read topological charge
- 
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}
+            
644def read_gf_coupling(path, prefix, c, dtr_cnfg=1, Zeuthen_flow=True, **kwargs):
+645    """Read the gradient flow coupling based on sfqcd gradient flow measurements. See 1607.06423 for details.
+646
+647    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.
+648
+649    Parameters
+650    ----------
+651    path : str
+652        path of the measurement files
+653    prefix : str
+654        prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat.
+655        Ignored if file names are passed explicitly via keyword files.
+656    c : double
+657        Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L.
+658    dtr_cnfg : int
+659        (optional) parameter that specifies the number of measurements
+660        between two configs.
+661        If it is not set, the distance between two measurements
+662        in the file is assumed to be the distance between two configurations.
+663    steps : int
+664        (optional) Distance between two configurations in units of trajectories /
+665         cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
+666    r_start : list
+667        list which contains the first config to be read for each replicum.
+668    r_stop : list
+669        list which contains the last config to be read for each replicum.
+670    files : list
+671        specify the exact files that need to be read
+672        from path, practical if e.g. only one replicum is needed
+673    names : list
+674        Alternative labeling for replicas/ensembles.
+675        Has to have the appropriate length.
+676    postfix : str
+677        postfix of the file to read, e.g. '.gfms.dat' for openQCD-files
+678    Zeuthen_flow : bool
+679        (optional) If True, the Zeuthen flow is used for the coupling. If False, the Wilson flow is used.
+680    """
 681
-682    return t * t * (5 / 3 * plaq - 1 / 12 * C2x1) / normdict[L]
+682    if c != 0.3:
+683        raise Exception("The required lattice norm is only implemented for c=0.3 at the moment.")
+684
+685    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)
+686    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)
+687    L = plaq.tag["L"]
+688    T = plaq.tag["T"]
+689
+690    if T != L:
+691        raise Exception("The required lattice norm is only implemented for T=L at the moment.")
+692
+693    if Zeuthen_flow is not True:
+694        raise Exception("The required lattice norm is only implemented for the Zeuthen flow at the moment.")
+695
+696    t = (c * L) ** 2 / 8
+697
+698    normdict = {4: 0.012341170468270,
+699                6: 0.010162691462430,
+700                8: 0.009031614807931,
+701                10: 0.008744966371393,
+702                12: 0.008650917856809,
+703                14: 8.611154391267955E-03,
+704                16: 0.008591758449508,
+705                20: 0.008575359627103,
+706                24: 0.008569387847540,
+707                28: 8.566803713382559E-03,
+708                32: 0.008565541650006,
+709                40: 8.564480684962046E-03,
+710                48: 8.564098025073460E-03,
+711                64: 8.563853943383087E-03}
+712
+713    return t * t * (5 / 3 * plaq - 1 / 12 * C2x1) / normdict[L]
 

 
 
@@ -2134,30 +2170,30 @@ postfix of the file to read, e.g. '.gfms.dat' for openQCD-files
- 
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
+            
 988def qtop_projection(qtop, target=0):
+ 989    """Returns the projection to the topological charge sector defined by target.
+ 990
+ 991    Parameters
+ 992    ----------
+ 993    path : Obs
+ 994        Topological charge.
+ 995    target : int
+ 996        Specifies the topological sector to be reweighted to (default 0)
+ 997
+ 998    Returns
+ 999    -------
+1000    reto : Obs
+1001        projection to the topological charge sector defined by target
+1002    """
+1003    if qtop.reweighted:
+1004        raise Exception('You can not use a reweighted observable for reweighting!')
+1005
+1006    proj_qtop = []
+1007    for n in qtop.deltas:
+1008        proj_qtop.append(np.array([1 if round(qtop.r_values[n] + q) == target else 0 for q in qtop.deltas[n]]))
+1009
+1010    reto = Obs(proj_qtop, qtop.names, idl=[qtop.idl[name] for name in qtop.names])
+1011    return reto
 

 
 
@@ -2193,62 +2229,62 @@ projection to the topological charge sector defined by 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)
+            
1014def read_qtop_sector(path, prefix, c, target=0, **kwargs):
+1015    """Constructs reweighting factors to a specified topological sector.
+1016
+1017    Parameters
+1018    ----------
+1019    path : str
+1020        path of the measurement files
+1021    prefix : str
+1022        prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat
+1023    c : double
+1024        Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L
+1025    target : int
+1026        Specifies the topological sector to be reweighted to (default 0)
+1027    dtr_cnfg : int
+1028        (optional) parameter that specifies the number of trajectories
+1029        between two configs.
+1030        if it is not set, the distance between two measurements
+1031        in the file is assumed to be the distance between two configurations.
+1032    steps : int
+1033        (optional) Distance between two configurations in units of trajectories /
+1034         cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
+1035    version : str
+1036        version string of the openQCD (sfqcd) version used to create
+1037        the ensemble. Default is 2.0. May also be set to sfqcd.
+1038    L : int
+1039        spatial length of the lattice in L/a.
+1040        HAS to be set if version != sfqcd, since openQCD does not provide
+1041        this in the header
+1042    r_start : list
+1043        offset of the first ensemble, making it easier to match
+1044        later on with other Obs
+1045    r_stop : list
+1046        last configurations that need to be read (per replicum)
+1047    files : list
+1048        specify the exact files that need to be read
+1049        from path, practical if e.g. only one replicum is needed
+1050    names : list
+1051        Alternative labeling for replicas/ensembles.
+1052        Has to have the appropriate length
+1053    Zeuthen_flow : bool
+1054        (optional) If True, the Zeuthen flow is used for Qtop. Only possible
+1055        for version=='sfqcd' If False, the Wilson flow is used.
+1056
+1057    Returns
+1058    -------
+1059    reto : Obs
+1060        projection to the topological charge sector defined by target
+1061    """
+1062
+1063    if not isinstance(target, int):
+1064        raise Exception("'target' has to be an integer.")
+1065
+1066    kwargs['integer_charge'] = True
+1067    qtop = read_qtop(path, prefix, c, **kwargs)
+1068
+1069    return qtop_projection(qtop, target=target)
 

 
 
@@ -2317,159 +2353,159 @@ projection to the topological charge sector defined by target
- 
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    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
+            
1072def read_ms5_xsf(path, prefix, qc, corr, sep="r", **kwargs):
+1073    """
+1074    Read data from files in the specified directory with the specified prefix and quark combination extension, and return a `Corr` object containing the data.
+1075
+1076    Parameters
+1077    ----------
+1078    path : str
+1079        The directory to search for the files in.
+1080    prefix : str
+1081        The prefix to match the files against.
+1082    qc : str
+1083        The quark combination extension to match the files against.
+1084    corr : str
+1085        The correlator to extract data for.
+1086    sep : str, optional
+1087        The separator to use when parsing the replika names.
+1088    **kwargs
+1089        Additional keyword arguments. The following keyword arguments are recognized:
+1090
+1091        - names (List[str]): A list of names to use for the replicas.
+1092
+1093    Returns
+1094    -------
+1095    Corr
+1096        A complex valued `Corr` object containing the data read from the files. In case of boudary to bulk correlators.
+1097    or
+1098    CObs
+1099        A complex valued `CObs` object containing the data read from the files. In case of boudary to boundary correlators.
+1100
+1101
+1102    Raises
+1103    ------
+1104    FileNotFoundError
+1105        If no files matching the specified prefix and quark combination extension are found in the specified directory.
+1106    IOError
+1107        If there is an error reading a file.
+1108    struct.error
+1109        If there is an error unpacking binary data.
+1110    """
+1111
+1112    # found = []
+1113    files = []
+1114    names = []
+1115
+1116    # test if the input is correct
+1117    if qc not in ['dd', 'ud', 'du', 'uu']:
+1118        raise Exception("Unknown quark conbination!")
+1119
+1120    if corr not in ["gS", "gP", "gA", "gV", "gVt", "lA", "lV", "lVt", "lT", "lTt", "g1", "l1"]:
+1121        raise Exception("Unknown correlator!")
+1122
+1123    if "files" in kwargs:
+1124        known_files = kwargs.get("files")
+1125    else:
+1126        known_files = []
+1127    files = _find_files(path, prefix, "ms5_xsf_" + qc, "dat", known_files=known_files)
+1128
+1129    if "names" in kwargs:
+1130        names = kwargs.get("names")
+1131    else:
+1132        for f in files:
+1133            if not sep == "":
+1134                se = f.split(".")[0]
+1135                for s in f.split(".")[1:-2]:
+1136                    se += "." + s
+1137                names.append(se.split(sep)[0] + "|r" + se.split(sep)[1])
+1138            else:
+1139                names.append(prefix)
 1140
-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            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
+1141    names = sorted(names)
+1142    files = sorted(files)
+1143
+1144    cnfgs = []
+1145    realsamples = []
+1146    imagsamples = []
+1147    repnum = 0
+1148    for file in files:
+1149        with open(path + "/" + file, "rb") as fp:
+1150
+1151            t = fp.read(8)
+1152            kappa = struct.unpack('d', t)[0]
+1153            t = fp.read(8)
+1154            csw = struct.unpack('d', t)[0]
+1155            t = fp.read(8)
+1156            dF = struct.unpack('d', t)[0]
+1157            t = fp.read(8)
+1158            zF = struct.unpack('d', t)[0]
+1159
+1160            t = fp.read(4)
+1161            tmax = struct.unpack('i', t)[0]
+1162            t = fp.read(4)
+1163            bnd = struct.unpack('i', t)[0]
+1164
+1165            placesBI = ["gS", "gP",
+1166                        "gA", "gV",
+1167                        "gVt", "lA",
+1168                        "lV", "lVt",
+1169                        "lT", "lTt"]
+1170            placesBB = ["g1", "l1"]
 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)
+1172            # the chunks have the following structure:
+1173            # confignumber, 10x timedependent complex correlators as doubles, 2x timeindependent complex correlators as doubles
+1174
+1175            chunksize = 4 + (8 * 2 * tmax * 10) + (8 * 2 * 2)
+1176            packstr = '=i' + ('d' * 2 * tmax * 10) + ('d' * 2 * 2)
+1177            cnfgs.append([])
+1178            realsamples.append([])
+1179            imagsamples.append([])
+1180            for t in range(tmax):
+1181                realsamples[repnum].append([])
+1182                imagsamples[repnum].append([])
+1183
+1184            while True:
+1185                cnfgt = fp.read(chunksize)
+1186                if not cnfgt:
+1187                    break
+1188                asascii = struct.unpack(packstr, cnfgt)
+1189                cnfg = asascii[0]
+1190                cnfgs[repnum].append(cnfg)
+1191
+1192                if corr not in placesBB:
+1193                    tmpcorr = asascii[1 + 2 * tmax * placesBI.index(corr):1 + 2 * tmax * placesBI.index(corr) + 2 * tmax]
+1194                else:
+1195                    tmpcorr = asascii[1 + 2 * tmax * len(placesBI) + 2 * placesBB.index(corr):1 + 2 * tmax * len(placesBI) + 2 * placesBB.index(corr) + 2]
+1196
+1197                corrres = [[], []]
+1198                for i in range(len(tmpcorr)):
+1199                    corrres[i % 2].append(tmpcorr[i])
+1200                for t in range(int(len(tmpcorr) / 2)):
+1201                    realsamples[repnum][t].append(corrres[0][t])
+1202                for t in range(int(len(tmpcorr) / 2)):
+1203                    imagsamples[repnum][t].append(corrres[1][t])
+1204        repnum += 1
+1205
+1206    s = "Read correlator " + corr + " from " + str(repnum) + " replika with " + str(len(realsamples[0][t]))
+1207    for rep in range(1, repnum):
+1208        s += ", " + str(len(realsamples[rep][t]))
+1209    s += " samples"
+1210    print(s)
+1211    print("Asserted run parameters:\n T:", tmax, "kappa:", kappa, "csw:", csw, "dF:", dF, "zF:", zF, "bnd:", bnd)
+1212
+1213    # we have the data now... but we need to re format the whole thing and put it into Corr objects.
+1214
+1215    compObs = []
+1216
+1217    for t in range(int(len(tmpcorr) / 2)):
+1218        compObs.append(CObs(Obs([realsamples[rep][t] for rep in range(repnum)], names=names, idl=cnfgs),
+1219                            Obs([imagsamples[rep][t] for rep in range(repnum)], names=names, idl=cnfgs)))
+1220
+1221    if len(compObs) == 1:
+1222        return compObs[0]
+1223    else:
+1224        return Corr(compObs)