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
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)
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]
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
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)
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)