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