pyerrors.input.openQCD
View Source
0import os 1import fnmatch 2import re 3import struct 4import warnings 5import numpy as np # Thinly-wrapped numpy 6import matplotlib.pyplot as plt 7from matplotlib import gridspec 8from ..obs import Obs 9from ..fits import fit_lin 10 11 12def read_rwms(path, prefix, version='2.0', names=None, **kwargs): 13 """Read rwms format from given folder structure. Returns a list of length nrw 14 15 Parameters 16 ---------- 17 path : str 18 path that contains the data files 19 prefix : str 20 all files in path that start with prefix are considered as input files. 21 May be used together postfix to consider only special file endings. 22 Prefix is ignored, if the keyword 'files' is used. 23 version : str 24 version of openQCD, default 2.0 25 names : list 26 list of names that is assigned to the data according according 27 to the order in the file list. Use careful, if you do not provide file names! 28 r_start : list 29 list which contains the first config to be read for each replicum 30 r_stop : list 31 list which contains the last config to be read for each replicum 32 r_step : int 33 integer that defines a fixed step size between two measurements (in units of configs) 34 If not given, r_step=1 is assumed. 35 postfix : str 36 postfix of the file to read, e.g. '.ms1' for openQCD-files 37 files : list 38 list which contains the filenames to be read. No automatic detection of 39 files performed if given. 40 print_err : bool 41 Print additional information that is useful for debugging. 42 """ 43 known_oqcd_versions = ['1.4', '1.6', '2.0'] 44 if not (version in known_oqcd_versions): 45 raise Exception('Unknown openQCD version defined!') 46 print("Working with openQCD version " + version) 47 if 'postfix' in kwargs: 48 postfix = kwargs.get('postfix') 49 else: 50 postfix = '' 51 ls = [] 52 for (dirpath, dirnames, filenames) in os.walk(path): 53 ls.extend(filenames) 54 break 55 56 if not ls: 57 raise Exception('Error, directory not found') 58 if 'files' in kwargs: 59 ls = kwargs.get('files') 60 else: 61 for exc in ls: 62 if not fnmatch.fnmatch(exc, prefix + '*' + postfix + '.dat'): 63 ls = list(set(ls) - set([exc])) 64 if len(ls) > 1: 65 ls.sort(key=lambda x: int(re.findall(r'\d+', x[len(prefix):])[0])) 66 replica = len(ls) 67 68 if 'r_start' in kwargs: 69 r_start = kwargs.get('r_start') 70 if len(r_start) != replica: 71 raise Exception('r_start does not match number of replicas') 72 r_start = [o if o else None for o in r_start] 73 else: 74 r_start = [None] * replica 75 76 if 'r_stop' in kwargs: 77 r_stop = kwargs.get('r_stop') 78 if len(r_stop) != replica: 79 raise Exception('r_stop does not match number of replicas') 80 else: 81 r_stop = [None] * replica 82 83 if 'r_step' in kwargs: 84 r_step = kwargs.get('r_step') 85 else: 86 r_step = 1 87 88 print('Read reweighting factors from', prefix[:-1], ',', 89 replica, 'replica', end='') 90 91 if names is None: 92 rep_names = [] 93 for entry in ls: 94 truncated_entry = entry 95 suffixes = [".dat", ".rwms", ".ms1"] 96 for suffix in suffixes: 97 if truncated_entry.endswith(suffix): 98 truncated_entry = truncated_entry[0:-len(suffix)] 99 idx = truncated_entry.index('r') 100 rep_names.append(truncated_entry[:idx] + '|' + truncated_entry[idx:]) 101 else: 102 rep_names = names 103 104 print_err = 0 105 if 'print_err' in kwargs: 106 print_err = 1 107 print() 108 109 deltas = [] 110 111 configlist = [] 112 r_start_index = [] 113 r_stop_index = [] 114 115 for rep in range(replica): 116 tmp_array = [] 117 with open(path + '/' + ls[rep], 'rb') as fp: 118 119 t = fp.read(4) # number of reweighting factors 120 if rep == 0: 121 nrw = struct.unpack('i', t)[0] 122 if version == '2.0': 123 nrw = int(nrw / 2) 124 for k in range(nrw): 125 deltas.append([]) 126 else: 127 if ((nrw != struct.unpack('i', t)[0] and (not version == '2.0')) or (nrw != struct.unpack('i', t)[0] / 2 and version == '2.0')): 128 raise Exception('Error: different number of reweighting factors for replicum', rep) 129 130 for k in range(nrw): 131 tmp_array.append([]) 132 133 # This block is necessary for openQCD1.6 and openQCD2.0 ms1 files 134 nfct = [] 135 if version in ['1.6', '2.0']: 136 for i in range(nrw): 137 t = fp.read(4) 138 nfct.append(struct.unpack('i', t)[0]) 139 else: 140 for i in range(nrw): 141 nfct.append(1) 142 143 nsrc = [] 144 for i in range(nrw): 145 t = fp.read(4) 146 nsrc.append(struct.unpack('i', t)[0]) 147 if version == '2.0': 148 if not struct.unpack('i', fp.read(4))[0] == 0: 149 print('something is wrong!') 150 151 configlist.append([]) 152 while 0 < 1: 153 t = fp.read(4) 154 if len(t) < 4: 155 break 156 config_no = struct.unpack('i', t)[0] 157 configlist[-1].append(config_no) 158 for i in range(nrw): 159 if(version == '2.0'): 160 tmpd = _read_array_openQCD2(fp) 161 tmpd = _read_array_openQCD2(fp) 162 tmp_rw = tmpd['arr'] 163 tmp_nfct = 1.0 164 for j in range(tmpd['n'][0]): 165 tmp_nfct *= np.mean(np.exp(-np.asarray(tmp_rw[j]))) 166 if print_err: 167 print(config_no, i, j, 168 np.mean(np.exp(-np.asarray(tmp_rw[j]))), 169 np.std(np.exp(-np.asarray(tmp_rw[j])))) 170 print('Sources:', 171 np.exp(-np.asarray(tmp_rw[j]))) 172 print('Partial factor:', tmp_nfct) 173 elif version == '1.6' or version == '1.4': 174 tmp_nfct = 1.0 175 for j in range(nfct[i]): 176 t = fp.read(8 * nsrc[i]) 177 t = fp.read(8 * nsrc[i]) 178 tmp_rw = struct.unpack('d' * nsrc[i], t) 179 tmp_nfct *= np.mean(np.exp(-np.asarray(tmp_rw))) 180 if print_err: 181 print(config_no, i, j, 182 np.mean(np.exp(-np.asarray(tmp_rw))), 183 np.std(np.exp(-np.asarray(tmp_rw)))) 184 print('Sources:', np.exp(-np.asarray(tmp_rw))) 185 print('Partial factor:', tmp_nfct) 186 tmp_array[i].append(tmp_nfct) 187 188 diffmeas = configlist[-1][-1] - configlist[-1][-2] 189 configlist[-1] = [item // diffmeas for item in configlist[-1]] 190 if configlist[-1][0] > 1 and diffmeas > 1: 191 warnings.warn('Assume thermalization and that the first measurement belongs to the first config.') 192 offset = configlist[-1][0] - 1 193 configlist[-1] = [item - offset for item in configlist[-1]] 194 195 if r_start[rep] is None: 196 r_start_index.append(0) 197 else: 198 try: 199 r_start_index.append(configlist[-1].index(r_start[rep])) 200 except ValueError: 201 raise Exception('Config %d not in file with range [%d, %d]' % ( 202 r_start[rep], configlist[-1][0], configlist[-1][-1])) from None 203 204 if r_stop[rep] is None: 205 r_stop_index.append(len(configlist[-1]) - 1) 206 else: 207 try: 208 r_stop_index.append(configlist[-1].index(r_stop[rep])) 209 except ValueError: 210 raise Exception('Config %d not in file with range [%d, %d]' % ( 211 r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None 212 213 for k in range(nrw): 214 deltas[k].append(tmp_array[k][r_start_index[rep]:r_stop_index[rep] + 1][::r_step]) 215 216 if np.any([len(np.unique(np.diff(cl))) != 1 for cl in configlist]): 217 raise Exception('Irregular spaced data in input file!', [len(np.unique(np.diff(cl))) for cl in configlist]) 218 stepsizes = [list(np.unique(np.diff(cl)))[0] for cl in configlist] 219 if np.any([step != 1 for step in stepsizes]): 220 warnings.warn('Stepsize between configurations is greater than one!' + str(stepsizes), RuntimeWarning) 221 222 print(',', nrw, 'reweighting factors with', nsrc, 'sources') 223 result = [] 224 idl = [range(configlist[rep][r_start_index[rep]], configlist[rep][r_stop_index[rep]] + 1, r_step) for rep in range(replica)] 225 226 for t in range(nrw): 227 result.append(Obs(deltas[t], rep_names, idl=idl)) 228 return result 229 230 231def extract_t0(path, prefix, dtr_read, xmin, spatial_extent, fit_range=5, **kwargs): 232 """Extract t0 from given .ms.dat files. Returns t0 as Obs. 233 234 It is assumed that all boundary effects have 235 sufficiently decayed at x0=xmin. 236 The data around the zero crossing of t^2<E> - 0.3 237 is fitted with a linear function 238 from which the exact root is extracted. 239 240 It is assumed that one measurement is performed for each config. 241 If this is not the case, the resulting idl, as well as the handling 242 of r_start, r_stop and r_step is wrong and the user has to correct 243 this in the resulting observable. 244 245 Parameters 246 ---------- 247 path : str 248 Path to .ms.dat files 249 prefix : str 250 Ensemble prefix 251 dtr_read : int 252 Determines how many trajectories should be skipped 253 when reading the ms.dat files. 254 Corresponds to dtr_cnfg / dtr_ms in the openQCD input file. 255 xmin : int 256 First timeslice where the boundary 257 effects have sufficiently decayed. 258 spatial_extent : int 259 spatial extent of the lattice, required for normalization. 260 fit_range : int 261 Number of data points left and right of the zero 262 crossing to be included in the linear fit. (Default: 5) 263 r_start : list 264 list which contains the first config to be read for each replicum. 265 r_stop : list 266 list which contains the last config to be read for each replicum. 267 r_step : int 268 integer that defines a fixed step size between two measurements (in units of configs) 269 If not given, r_step=1 is assumed. 270 plaquette : bool 271 If true extract the plaquette estimate of t0 instead. 272 names : list 273 list of names that is assigned to the data according according 274 to the order in the file list. Use careful, if you do not provide file names! 275 files : list 276 list which contains the filenames to be read. No automatic detection of 277 files performed if given. 278 plot_fit : bool 279 If true, the fit for the extraction of t0 is shown together with the data. 280 assume_thermalization : bool 281 If True: If the first record divided by the distance between two measurements is larger than 282 1, it is assumed that this is due to thermalization and the first measurement belongs 283 to the first config (default). 284 If False: The config numbers are assumed to be traj_number // difference 285 """ 286 287 ls = [] 288 for (dirpath, dirnames, filenames) in os.walk(path): 289 ls.extend(filenames) 290 break 291 292 if not ls: 293 raise Exception('Error, directory not found') 294 295 if 'files' in kwargs: 296 ls = kwargs.get('files') 297 else: 298 for exc in ls: 299 if not fnmatch.fnmatch(exc, prefix + '*.ms.dat'): 300 ls = list(set(ls) - set([exc])) 301 if len(ls) > 1: 302 ls.sort(key=lambda x: int(re.findall(r'\d+', x[len(prefix):])[0])) 303 replica = len(ls) 304 305 if 'r_start' in kwargs: 306 r_start = kwargs.get('r_start') 307 if len(r_start) != replica: 308 raise Exception('r_start does not match number of replicas') 309 r_start = [o if o else None for o in r_start] 310 else: 311 r_start = [None] * replica 312 313 if 'r_stop' in kwargs: 314 r_stop = kwargs.get('r_stop') 315 if len(r_stop) != replica: 316 raise Exception('r_stop does not match number of replicas') 317 else: 318 r_stop = [None] * replica 319 320 if 'r_step' in kwargs: 321 r_step = kwargs.get('r_step') 322 else: 323 r_step = 1 324 325 print('Extract t0 from', prefix, ',', replica, 'replica') 326 327 if 'names' in kwargs: 328 rep_names = kwargs.get('names') 329 else: 330 rep_names = [] 331 for entry in ls: 332 truncated_entry = entry.split('.')[0] 333 idx = truncated_entry.index('r') 334 rep_names.append(truncated_entry[:idx] + '|' + truncated_entry[idx:]) 335 336 Ysum = [] 337 338 configlist = [] 339 r_start_index = [] 340 r_stop_index = [] 341 342 for rep in range(replica): 343 344 with open(path + '/' + ls[rep], 'rb') as fp: 345 t = fp.read(12) 346 header = struct.unpack('iii', t) 347 if rep == 0: 348 dn = header[0] 349 nn = header[1] 350 tmax = header[2] 351 elif dn != header[0] or nn != header[1] or tmax != header[2]: 352 raise Exception('Replica parameters do not match.') 353 354 t = fp.read(8) 355 if rep == 0: 356 eps = struct.unpack('d', t)[0] 357 print('Step size:', eps, ', Maximal t value:', dn * (nn) * eps) 358 elif eps != struct.unpack('d', t)[0]: 359 raise Exception('Values for eps do not match among replica.') 360 361 Ysl = [] 362 363 configlist.append([]) 364 while 0 < 1: 365 t = fp.read(4) 366 if(len(t) < 4): 367 break 368 nc = struct.unpack('i', t)[0] 369 configlist[-1].append(nc) 370 371 t = fp.read(8 * tmax * (nn + 1)) 372 if kwargs.get('plaquette'): 373 if nc % dtr_read == 0: 374 Ysl.append(struct.unpack('d' * tmax * (nn + 1), t)) 375 t = fp.read(8 * tmax * (nn + 1)) 376 if not kwargs.get('plaquette'): 377 if nc % dtr_read == 0: 378 Ysl.append(struct.unpack('d' * tmax * (nn + 1), t)) 379 t = fp.read(8 * tmax * (nn + 1)) 380 381 Ysum.append([]) 382 for i, item in enumerate(Ysl): 383 Ysum[-1].append([np.mean(item[current + xmin: 384 current + tmax - xmin]) 385 for current in range(0, len(item), tmax)]) 386 387 diffmeas = configlist[-1][-1] - configlist[-1][-2] 388 configlist[-1] = [item // diffmeas for item in configlist[-1]] 389 if kwargs.get('assume_thermalization', True) and configlist[-1][0] > 1: 390 warnings.warn('Assume thermalization and that the first measurement belongs to the first config.') 391 offset = configlist[-1][0] - 1 392 configlist[-1] = [item - offset for item in configlist[-1]] 393 394 if r_start[rep] is None: 395 r_start_index.append(0) 396 else: 397 try: 398 r_start_index.append(configlist[-1].index(r_start[rep])) 399 except ValueError: 400 raise Exception('Config %d not in file with range [%d, %d]' % ( 401 r_start[rep], configlist[-1][0], configlist[-1][-1])) from None 402 403 if r_stop[rep] is None: 404 r_stop_index.append(len(configlist[-1]) - 1) 405 else: 406 try: 407 r_stop_index.append(configlist[-1].index(r_stop[rep])) 408 except ValueError: 409 raise Exception('Config %d not in file with range [%d, %d]' % ( 410 r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None 411 412 if np.any([len(np.unique(np.diff(cl))) != 1 for cl in configlist]): 413 raise Exception('Irregular spaced data in input file!', [len(np.unique(np.diff(cl))) for cl in configlist]) 414 stepsizes = [list(np.unique(np.diff(cl)))[0] for cl in configlist] 415 if np.any([step != 1 for step in stepsizes]): 416 warnings.warn('Stepsize between configurations is greater than one!' + str(stepsizes), RuntimeWarning) 417 418 idl = [range(configlist[rep][r_start_index[rep]], configlist[rep][r_stop_index[rep]] + 1, r_step) for rep in range(replica)] 419 t2E_dict = {} 420 for n in range(nn + 1): 421 samples = [] 422 for nrep, rep in enumerate(Ysum): 423 samples.append([]) 424 for cnfg in rep: 425 samples[-1].append(cnfg[n]) 426 samples[-1] = samples[-1][r_start_index[nrep]:r_stop_index[nrep] + 1][::r_step] 427 new_obs = Obs(samples, rep_names, idl=idl) 428 t2E_dict[n * dn * eps] = (n * dn * eps) ** 2 * new_obs / (spatial_extent ** 3) - 0.3 429 430 zero_crossing = np.argmax(np.array( 431 [o.value for o in t2E_dict.values()]) > 0.0) 432 433 x = list(t2E_dict.keys())[zero_crossing - fit_range: 434 zero_crossing + fit_range] 435 y = list(t2E_dict.values())[zero_crossing - fit_range: 436 zero_crossing + fit_range] 437 [o.gamma_method() for o in y] 438 439 fit_result = fit_lin(x, y) 440 441 if kwargs.get('plot_fit'): 442 plt.figure() 443 gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1], wspace=0.0, hspace=0.0) 444 ax0 = plt.subplot(gs[0]) 445 xmore = list(t2E_dict.keys())[zero_crossing - fit_range - 2: zero_crossing + fit_range + 2] 446 ymore = list(t2E_dict.values())[zero_crossing - fit_range - 2: zero_crossing + fit_range + 2] 447 [o.gamma_method() for o in ymore] 448 ax0.errorbar(xmore, [yi.value for yi in ymore], yerr=[yi.dvalue for yi in ymore], fmt='x') 449 xplot = np.linspace(np.min(x), np.max(x)) 450 yplot = [fit_result[0] + fit_result[1] * xi for xi in xplot] 451 [yi.gamma_method() for yi in yplot] 452 ax0.fill_between(xplot, y1=[yi.value - yi.dvalue for yi in yplot], y2=[yi.value + yi.dvalue for yi in yplot]) 453 retval = (-fit_result[0] / fit_result[1]) 454 retval.gamma_method() 455 ylim = ax0.get_ylim() 456 ax0.fill_betweenx(ylim, x1=retval.value - retval.dvalue, x2=retval.value + retval.dvalue, color='gray', alpha=0.4) 457 ax0.set_ylim(ylim) 458 ax0.set_ylabel(r'$t^2 \langle E(t) \rangle - 0.3 $') 459 xlim = ax0.get_xlim() 460 461 fit_res = [fit_result[0] + fit_result[1] * xi for xi in x] 462 residuals = (np.asarray([o.value for o in y]) - [o.value for o in fit_res]) / np.asarray([o.dvalue for o in y]) 463 ax1 = plt.subplot(gs[1]) 464 ax1.plot(x, residuals, 'ko', ls='none', markersize=5) 465 ax1.tick_params(direction='out') 466 ax1.tick_params(axis="x", bottom=True, top=True, labelbottom=True) 467 ax1.axhline(y=0.0, ls='--', color='k') 468 ax1.fill_between(xlim, -1.0, 1.0, alpha=0.1, facecolor='k') 469 ax1.set_xlim(xlim) 470 ax1.set_ylabel('Residuals') 471 ax1.set_xlabel(r'$t/a^2$') 472 473 plt.draw() 474 return -fit_result[0] / fit_result[1] 475 476 477def _parse_array_openQCD2(d, n, size, wa, quadrupel=False): 478 arr = [] 479 if d == 2: 480 for i in range(n[0]): 481 tmp = wa[i * n[1]:(i + 1) * n[1]] 482 if quadrupel: 483 tmp2 = [] 484 for j in range(0, len(tmp), 2): 485 tmp2.append(tmp[j]) 486 arr.append(tmp2) 487 else: 488 arr.append(np.asarray(tmp)) 489 490 else: 491 raise Exception('Only two-dimensional arrays supported!') 492 493 return arr 494 495 496def _read_array_openQCD2(fp): 497 t = fp.read(4) 498 d = struct.unpack('i', t)[0] 499 t = fp.read(4 * d) 500 n = struct.unpack('%di' % (d), t) 501 t = fp.read(4) 502 size = struct.unpack('i', t)[0] 503 if size == 4: 504 types = 'i' 505 elif size == 8: 506 types = 'd' 507 elif size == 16: 508 types = 'dd' 509 else: 510 raise Exception("Type for size '" + str(size) + "' not known.") 511 m = n[0] 512 for i in range(1, d): 513 m *= n[i] 514 515 t = fp.read(m * size) 516 tmp = struct.unpack('%d%s' % (m, types), t) 517 518 arr = _parse_array_openQCD2(d, n, size, tmp, quadrupel=True) 519 return {'d': d, 'n': n, 'size': size, 'arr': arr} 520 521 522def read_qtop(path, prefix, c, dtr_cnfg=1, version="openQCD", **kwargs): 523 """Read the topologial charge based on openQCD gradient flow measurements. 524 525 Parameters 526 ---------- 527 path : str 528 path of the measurement files 529 prefix : str 530 prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat. 531 Ignored if file names are passed explicitly via keyword files. 532 c : double 533 Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L. 534 dtr_cnfg : int 535 (optional) parameter that specifies the number of measurements 536 between two configs. 537 If it is not set, the distance between two measurements 538 in the file is assumed to be the distance between two configurations. 539 steps : int 540 (optional) Distance between two configurations in units of trajectories / 541 cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given 542 version : str 543 Either openQCD or sfqcd, depending on the data. 544 L : int 545 spatial length of the lattice in L/a. 546 HAS to be set if version != sfqcd, since openQCD does not provide 547 this in the header 548 r_start : list 549 list which contains the first config to be read for each replicum. 550 r_stop : list 551 list which contains the last config to be read for each replicum. 552 files : list 553 specify the exact files that need to be read 554 from path, practical if e.g. only one replicum is needed 555 names : list 556 Alternative labeling for replicas/ensembles. 557 Has to have the appropriate length. 558 Zeuthen_flow : bool 559 (optional) If True, the Zeuthen flow is used for Qtop. Only possible 560 for version=='sfqcd' If False, the Wilson flow is used. 561 integer_charge : bool 562 If True, the charge is rounded towards the nearest integer on each config. 563 """ 564 known_versions = ["openQCD", "sfqcd"] 565 566 if version not in known_versions: 567 raise Exception("Unknown openQCD version.") 568 if "steps" in kwargs: 569 steps = kwargs.get("steps") 570 if version == "sfqcd": 571 if "L" in kwargs: 572 supposed_L = kwargs.get("L") 573 else: 574 supposed_L = None 575 postfix = ".gfms.dat" 576 else: 577 if "L" not in kwargs: 578 raise Exception("This version of openQCD needs you to provide the spatial length of the lattice as parameter 'L'.") 579 else: 580 L = kwargs.get("L") 581 postfix = ".ms.dat" 582 583 if "files" in kwargs: 584 files = kwargs.get("files") 585 postfix = '' 586 else: 587 found = [] 588 files = [] 589 for (dirpath, dirnames, filenames) in os.walk(path + "/"): 590 found.extend(filenames) 591 break 592 for f in found: 593 if fnmatch.fnmatch(f, prefix + "*" + postfix): 594 files.append(f) 595 596 if 'r_start' in kwargs: 597 r_start = kwargs.get('r_start') 598 if len(r_start) != len(files): 599 raise Exception('r_start does not match number of replicas') 600 r_start = [o if o else None for o in r_start] 601 else: 602 r_start = [None] * len(files) 603 604 if 'r_stop' in kwargs: 605 r_stop = kwargs.get('r_stop') 606 if len(r_stop) != len(files): 607 raise Exception('r_stop does not match number of replicas') 608 else: 609 r_stop = [None] * len(files) 610 rep_names = [] 611 612 zeuthen = kwargs.get('Zeuthen_flow', False) 613 if zeuthen and version not in ['sfqcd']: 614 raise Exception('Zeuthen flow can only be used for version==sfqcd') 615 616 r_start_index = [] 617 r_stop_index = [] 618 deltas = [] 619 configlist = [] 620 for rep, file in enumerate(files): 621 with open(path + "/" + file, "rb") as fp: 622 623 Q = [] 624 traj_list = [] 625 if version in ['sfqcd']: 626 if zeuthen: 627 obspos = 0 628 else: 629 obspos = 8 630 t = fp.read(12) 631 header = struct.unpack('<iii', t) 632 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) 633 ncs = header[1] # number of different values for c in t_flow=1/8 c² L² -> measurements done for ncs c's 634 tmax = header[2] # lattice T/a 635 636 t = fp.read(12) 637 Ls = struct.unpack('<iii', t) 638 if(Ls[0] == Ls[1] and Ls[1] == Ls[2]): 639 L = Ls[0] 640 if not (supposed_L == L) and supposed_L: 641 raise Exception("It seems the length given in the header and by you contradict each other") 642 else: 643 raise Exception("Found more than one spatial length in header!") 644 645 t = fp.read(16) 646 header2 = struct.unpack('<dd', t) 647 tol = header2[0] 648 cmax = header2[1] # highest value of c used 649 650 if c > cmax: 651 raise Exception('Flow has been determined between c=0 and c=%lf with tolerance %lf' % (cmax, tol)) 652 653 if(zthfl == 2): 654 nfl = 2 # number of flows 655 else: 656 nfl = 1 657 iobs = 8 * nfl # number of flow observables calculated 658 659 while 0 < 1: 660 t = fp.read(4) 661 if(len(t) < 4): 662 break 663 traj_list.append(struct.unpack('i', t)[0]) # trajectory number when measurement was done 664 665 for j in range(ncs + 1): 666 for i in range(iobs): 667 t = fp.read(8 * tmax) 668 if (i == obspos): # determines the flow observable -> i=0 <-> Zeuthen flow 669 Q.append(struct.unpack('d' * tmax, t)) 670 671 else: 672 t = fp.read(12) 673 header = struct.unpack('<iii', t) 674 # step size in integration steps "dnms" 675 dn = header[0] 676 # number of measurements, so "ntot"/dn 677 nn = header[1] 678 # lattice T/a 679 tmax = header[2] 680 681 t = fp.read(8) 682 eps = struct.unpack('d', t)[0] 683 684 while 0 < 1: 685 t = fp.read(4) 686 if(len(t) < 4): 687 break 688 traj_list.append(struct.unpack('i', t)[0]) 689 # Wsl 690 t = fp.read(8 * tmax * (nn + 1)) 691 # Ysl 692 t = fp.read(8 * tmax * (nn + 1)) 693 # Qsl, which is asked for in this method 694 t = fp.read(8 * tmax * (nn + 1)) 695 # unpack the array of Qtops, 696 # on each timeslice t=0,...,tmax-1 and the 697 # measurement number in = 0...nn (see README.qcd1) 698 tmpd = struct.unpack('d' * tmax * (nn + 1), t) 699 Q.append(tmpd) 700 701 if len(np.unique(np.diff(traj_list))) != 1: 702 raise Exception("Irregularities in stepsize found") 703 else: 704 if 'steps' in kwargs: 705 if steps != traj_list[1] - traj_list[0]: 706 raise Exception("steps and the found stepsize are not the same") 707 else: 708 steps = traj_list[1] - traj_list[0] 709 710 configlist.append([tr // steps // dtr_cnfg for tr in traj_list]) 711 if configlist[-1][0] > 1: 712 offset = configlist[-1][0] - 1 713 warnings.warn('Assume thermalization and that the first measurement belongs to the first config. Offset = %d configs (%d trajectories / cycles)' % ( 714 offset, offset * steps)) 715 configlist[-1] = [item - offset for item in configlist[-1]] 716 717 if r_start[rep] is None: 718 r_start_index.append(0) 719 else: 720 try: 721 r_start_index.append(configlist[-1].index(r_start[rep])) 722 except ValueError: 723 raise Exception('Config %d not in file with range [%d, %d]' % ( 724 r_start[rep], configlist[-1][0], configlist[-1][-1])) from None 725 726 if r_stop[rep] is None: 727 r_stop_index.append(len(configlist[-1]) - 1) 728 else: 729 try: 730 r_stop_index.append(configlist[-1].index(r_stop[rep])) 731 except ValueError: 732 raise Exception('Config %d not in file with range [%d, %d]' % ( 733 r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None 734 735 if version in ['sfqcd']: 736 cstepsize = cmax / ncs 737 index_aim = round(c / cstepsize) 738 else: 739 t_aim = (c * L) ** 2 / 8 740 index_aim = round(t_aim / eps / dn) 741 742 Q_sum = [] 743 for i, item in enumerate(Q): 744 Q_sum.append([sum(item[current:current + tmax]) 745 for current in range(0, len(item), tmax)]) 746 Q_top = [] 747 if version in ['sfqcd']: 748 for i in range(len(Q_sum) // (ncs + 1)): 749 Q_top.append(Q_sum[i * (ncs + 1) + index_aim][0]) 750 else: 751 for i in range(len(Q) // dtr_cnfg): 752 Q_top.append(Q_sum[dtr_cnfg * i][index_aim]) 753 if len(Q_top) != len(traj_list) // dtr_cnfg: 754 raise Exception("qtops and traj_list dont have the same length") 755 756 if kwargs.get('integer_charge', False): 757 Q_top = [round(q) for q in Q_top] 758 759 truncated_file = file[:-len(postfix)] 760 761 if "names" not in kwargs: 762 try: 763 idx = truncated_file.index('r') 764 except Exception: 765 if "names" not in kwargs: 766 raise Exception("Automatic recognition of replicum failed, please enter the key word 'names'.") 767 ens_name = truncated_file[:idx] 768 rep_names.append(ens_name + '|' + truncated_file[idx:]) 769 else: 770 names = kwargs.get("names") 771 rep_names = names 772 deltas.append(Q_top) 773 774 idl = [range(int(configlist[rep][r_start_index[rep]]), int(configlist[rep][r_stop_index[rep]]) + 1, 1) for rep in range(len(deltas))] 775 deltas = [deltas[nrep][r_start_index[nrep]:r_stop_index[nrep] + 1] for nrep in range(len(deltas))] 776 result = Obs(deltas, rep_names, idl=idl) 777 return result 778 779 780def qtop_projection(qtop, target=0): 781 """Returns the projection to the topological charge sector defined by target. 782 783 Parameters 784 ---------- 785 path : Obs 786 Topological charge. 787 target : int 788 Specifies the topological sector to be reweighted to (default 0) 789 """ 790 if qtop.reweighted: 791 raise Exception('You can not use a reweighted observable for reweighting!') 792 793 proj_qtop = [] 794 for n in qtop.deltas: 795 proj_qtop.append(np.array([1 if round(qtop.value + q) == target else 0 for q in qtop.deltas[n]])) 796 797 reto = Obs(proj_qtop, qtop.names, idl=[qtop.idl[name] for name in qtop.names]) 798 reto.is_merged = qtop.is_merged 799 return reto 800 801 802def read_qtop_sector(path, prefix, c, target=0, **kwargs): 803 """Constructs reweighting factors to a specified topological sector. 804 805 Parameters 806 ---------- 807 path : str 808 path of the measurement files 809 prefix : str 810 prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat 811 c : double 812 Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L 813 target : int 814 Specifies the topological sector to be reweighted to (default 0) 815 dtr_cnfg : int 816 (optional) parameter that specifies the number of trajectories 817 between two configs. 818 if it is not set, the distance between two measurements 819 in the file is assumed to be the distance between two configurations. 820 steps : int 821 (optional) Distance between two configurations in units of trajectories / 822 cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given 823 version : str 824 version string of the openQCD (sfqcd) version used to create 825 the ensemble. Default is 2.0. May also be set to sfqcd. 826 L : int 827 spatial length of the lattice in L/a. 828 HAS to be set if version != sfqcd, since openQCD does not provide 829 this in the header 830 r_start : list 831 offset of the first ensemble, making it easier to match 832 later on with other Obs 833 r_stop : list 834 last configurations that need to be read (per replicum) 835 files : list 836 specify the exact files that need to be read 837 from path, practical if e.g. only one replicum is needed 838 names : list 839 Alternative labeling for replicas/ensembles. 840 Has to have the appropriate length 841 Zeuthen_flow : bool 842 (optional) If True, the Zeuthen flow is used for Qtop. Only possible 843 for version=='sfqcd' If False, the Wilson flow is used. 844 """ 845 846 if not isinstance(target, int): 847 raise Exception("'target' has to be an integer.") 848 849 kwargs['integer_charge'] = True 850 qtop = read_qtop(path, prefix, c, **kwargs) 851 852 return qtop_projection(qtop, target=target)
View Source
13def read_rwms(path, prefix, version='2.0', names=None, **kwargs): 14 """Read rwms format from given folder structure. Returns a list of length nrw 15 16 Parameters 17 ---------- 18 path : str 19 path that contains the data files 20 prefix : str 21 all files in path that start with prefix are considered as input files. 22 May be used together postfix to consider only special file endings. 23 Prefix is ignored, if the keyword 'files' is used. 24 version : str 25 version of openQCD, default 2.0 26 names : list 27 list of names that is assigned to the data according according 28 to the order in the file list. Use careful, if you do not provide file names! 29 r_start : list 30 list which contains the first config to be read for each replicum 31 r_stop : list 32 list which contains the last config to be read for each replicum 33 r_step : int 34 integer that defines a fixed step size between two measurements (in units of configs) 35 If not given, r_step=1 is assumed. 36 postfix : str 37 postfix of the file to read, e.g. '.ms1' for openQCD-files 38 files : list 39 list which contains the filenames to be read. No automatic detection of 40 files performed if given. 41 print_err : bool 42 Print additional information that is useful for debugging. 43 """ 44 known_oqcd_versions = ['1.4', '1.6', '2.0'] 45 if not (version in known_oqcd_versions): 46 raise Exception('Unknown openQCD version defined!') 47 print("Working with openQCD version " + version) 48 if 'postfix' in kwargs: 49 postfix = kwargs.get('postfix') 50 else: 51 postfix = '' 52 ls = [] 53 for (dirpath, dirnames, filenames) in os.walk(path): 54 ls.extend(filenames) 55 break 56 57 if not ls: 58 raise Exception('Error, directory not found') 59 if 'files' in kwargs: 60 ls = kwargs.get('files') 61 else: 62 for exc in ls: 63 if not fnmatch.fnmatch(exc, prefix + '*' + postfix + '.dat'): 64 ls = list(set(ls) - set([exc])) 65 if len(ls) > 1: 66 ls.sort(key=lambda x: int(re.findall(r'\d+', x[len(prefix):])[0])) 67 replica = len(ls) 68 69 if 'r_start' in kwargs: 70 r_start = kwargs.get('r_start') 71 if len(r_start) != replica: 72 raise Exception('r_start does not match number of replicas') 73 r_start = [o if o else None for o in r_start] 74 else: 75 r_start = [None] * replica 76 77 if 'r_stop' in kwargs: 78 r_stop = kwargs.get('r_stop') 79 if len(r_stop) != replica: 80 raise Exception('r_stop does not match number of replicas') 81 else: 82 r_stop = [None] * replica 83 84 if 'r_step' in kwargs: 85 r_step = kwargs.get('r_step') 86 else: 87 r_step = 1 88 89 print('Read reweighting factors from', prefix[:-1], ',', 90 replica, 'replica', end='') 91 92 if names is None: 93 rep_names = [] 94 for entry in ls: 95 truncated_entry = entry 96 suffixes = [".dat", ".rwms", ".ms1"] 97 for suffix in suffixes: 98 if truncated_entry.endswith(suffix): 99 truncated_entry = truncated_entry[0:-len(suffix)] 100 idx = truncated_entry.index('r') 101 rep_names.append(truncated_entry[:idx] + '|' + truncated_entry[idx:]) 102 else: 103 rep_names = names 104 105 print_err = 0 106 if 'print_err' in kwargs: 107 print_err = 1 108 print() 109 110 deltas = [] 111 112 configlist = [] 113 r_start_index = [] 114 r_stop_index = [] 115 116 for rep in range(replica): 117 tmp_array = [] 118 with open(path + '/' + ls[rep], 'rb') as fp: 119 120 t = fp.read(4) # number of reweighting factors 121 if rep == 0: 122 nrw = struct.unpack('i', t)[0] 123 if version == '2.0': 124 nrw = int(nrw / 2) 125 for k in range(nrw): 126 deltas.append([]) 127 else: 128 if ((nrw != struct.unpack('i', t)[0] and (not version == '2.0')) or (nrw != struct.unpack('i', t)[0] / 2 and version == '2.0')): 129 raise Exception('Error: different number of reweighting factors for replicum', rep) 130 131 for k in range(nrw): 132 tmp_array.append([]) 133 134 # This block is necessary for openQCD1.6 and openQCD2.0 ms1 files 135 nfct = [] 136 if version in ['1.6', '2.0']: 137 for i in range(nrw): 138 t = fp.read(4) 139 nfct.append(struct.unpack('i', t)[0]) 140 else: 141 for i in range(nrw): 142 nfct.append(1) 143 144 nsrc = [] 145 for i in range(nrw): 146 t = fp.read(4) 147 nsrc.append(struct.unpack('i', t)[0]) 148 if version == '2.0': 149 if not struct.unpack('i', fp.read(4))[0] == 0: 150 print('something is wrong!') 151 152 configlist.append([]) 153 while 0 < 1: 154 t = fp.read(4) 155 if len(t) < 4: 156 break 157 config_no = struct.unpack('i', t)[0] 158 configlist[-1].append(config_no) 159 for i in range(nrw): 160 if(version == '2.0'): 161 tmpd = _read_array_openQCD2(fp) 162 tmpd = _read_array_openQCD2(fp) 163 tmp_rw = tmpd['arr'] 164 tmp_nfct = 1.0 165 for j in range(tmpd['n'][0]): 166 tmp_nfct *= np.mean(np.exp(-np.asarray(tmp_rw[j]))) 167 if print_err: 168 print(config_no, i, j, 169 np.mean(np.exp(-np.asarray(tmp_rw[j]))), 170 np.std(np.exp(-np.asarray(tmp_rw[j])))) 171 print('Sources:', 172 np.exp(-np.asarray(tmp_rw[j]))) 173 print('Partial factor:', tmp_nfct) 174 elif version == '1.6' or version == '1.4': 175 tmp_nfct = 1.0 176 for j in range(nfct[i]): 177 t = fp.read(8 * nsrc[i]) 178 t = fp.read(8 * nsrc[i]) 179 tmp_rw = struct.unpack('d' * nsrc[i], t) 180 tmp_nfct *= np.mean(np.exp(-np.asarray(tmp_rw))) 181 if print_err: 182 print(config_no, i, j, 183 np.mean(np.exp(-np.asarray(tmp_rw))), 184 np.std(np.exp(-np.asarray(tmp_rw)))) 185 print('Sources:', np.exp(-np.asarray(tmp_rw))) 186 print('Partial factor:', tmp_nfct) 187 tmp_array[i].append(tmp_nfct) 188 189 diffmeas = configlist[-1][-1] - configlist[-1][-2] 190 configlist[-1] = [item // diffmeas for item in configlist[-1]] 191 if configlist[-1][0] > 1 and diffmeas > 1: 192 warnings.warn('Assume thermalization and that the first measurement belongs to the first config.') 193 offset = configlist[-1][0] - 1 194 configlist[-1] = [item - offset for item in configlist[-1]] 195 196 if r_start[rep] is None: 197 r_start_index.append(0) 198 else: 199 try: 200 r_start_index.append(configlist[-1].index(r_start[rep])) 201 except ValueError: 202 raise Exception('Config %d not in file with range [%d, %d]' % ( 203 r_start[rep], configlist[-1][0], configlist[-1][-1])) from None 204 205 if r_stop[rep] is None: 206 r_stop_index.append(len(configlist[-1]) - 1) 207 else: 208 try: 209 r_stop_index.append(configlist[-1].index(r_stop[rep])) 210 except ValueError: 211 raise Exception('Config %d not in file with range [%d, %d]' % ( 212 r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None 213 214 for k in range(nrw): 215 deltas[k].append(tmp_array[k][r_start_index[rep]:r_stop_index[rep] + 1][::r_step]) 216 217 if np.any([len(np.unique(np.diff(cl))) != 1 for cl in configlist]): 218 raise Exception('Irregular spaced data in input file!', [len(np.unique(np.diff(cl))) for cl in configlist]) 219 stepsizes = [list(np.unique(np.diff(cl)))[0] for cl in configlist] 220 if np.any([step != 1 for step in stepsizes]): 221 warnings.warn('Stepsize between configurations is greater than one!' + str(stepsizes), RuntimeWarning) 222 223 print(',', nrw, 'reweighting factors with', nsrc, 'sources') 224 result = [] 225 idl = [range(configlist[rep][r_start_index[rep]], configlist[rep][r_stop_index[rep]] + 1, r_step) for rep in range(replica)] 226 227 for t in range(nrw): 228 result.append(Obs(deltas[t], rep_names, idl=idl)) 229 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.
View Source
232def extract_t0(path, prefix, dtr_read, xmin, spatial_extent, fit_range=5, **kwargs): 233 """Extract t0 from given .ms.dat files. Returns t0 as Obs. 234 235 It is assumed that all boundary effects have 236 sufficiently decayed at x0=xmin. 237 The data around the zero crossing of t^2<E> - 0.3 238 is fitted with a linear function 239 from which the exact root is extracted. 240 241 It is assumed that one measurement is performed for each config. 242 If this is not the case, the resulting idl, as well as the handling 243 of r_start, r_stop and r_step is wrong and the user has to correct 244 this in the resulting observable. 245 246 Parameters 247 ---------- 248 path : str 249 Path to .ms.dat files 250 prefix : str 251 Ensemble prefix 252 dtr_read : int 253 Determines how many trajectories should be skipped 254 when reading the ms.dat files. 255 Corresponds to dtr_cnfg / dtr_ms in the openQCD input file. 256 xmin : int 257 First timeslice where the boundary 258 effects have sufficiently decayed. 259 spatial_extent : int 260 spatial extent of the lattice, required for normalization. 261 fit_range : int 262 Number of data points left and right of the zero 263 crossing to be included in the linear fit. (Default: 5) 264 r_start : list 265 list which contains the first config to be read for each replicum. 266 r_stop : list 267 list which contains the last config to be read for each replicum. 268 r_step : int 269 integer that defines a fixed step size between two measurements (in units of configs) 270 If not given, r_step=1 is assumed. 271 plaquette : bool 272 If true extract the plaquette estimate of t0 instead. 273 names : list 274 list of names that is assigned to the data according according 275 to the order in the file list. Use careful, if you do not provide file names! 276 files : list 277 list which contains the filenames to be read. No automatic detection of 278 files performed if given. 279 plot_fit : bool 280 If true, the fit for the extraction of t0 is shown together with the data. 281 assume_thermalization : bool 282 If True: If the first record divided by the distance between two measurements is larger than 283 1, it is assumed that this is due to thermalization and the first measurement belongs 284 to the first config (default). 285 If False: The config numbers are assumed to be traj_number // difference 286 """ 287 288 ls = [] 289 for (dirpath, dirnames, filenames) in os.walk(path): 290 ls.extend(filenames) 291 break 292 293 if not ls: 294 raise Exception('Error, directory not found') 295 296 if 'files' in kwargs: 297 ls = kwargs.get('files') 298 else: 299 for exc in ls: 300 if not fnmatch.fnmatch(exc, prefix + '*.ms.dat'): 301 ls = list(set(ls) - set([exc])) 302 if len(ls) > 1: 303 ls.sort(key=lambda x: int(re.findall(r'\d+', x[len(prefix):])[0])) 304 replica = len(ls) 305 306 if 'r_start' in kwargs: 307 r_start = kwargs.get('r_start') 308 if len(r_start) != replica: 309 raise Exception('r_start does not match number of replicas') 310 r_start = [o if o else None for o in r_start] 311 else: 312 r_start = [None] * replica 313 314 if 'r_stop' in kwargs: 315 r_stop = kwargs.get('r_stop') 316 if len(r_stop) != replica: 317 raise Exception('r_stop does not match number of replicas') 318 else: 319 r_stop = [None] * replica 320 321 if 'r_step' in kwargs: 322 r_step = kwargs.get('r_step') 323 else: 324 r_step = 1 325 326 print('Extract t0 from', prefix, ',', replica, 'replica') 327 328 if 'names' in kwargs: 329 rep_names = kwargs.get('names') 330 else: 331 rep_names = [] 332 for entry in ls: 333 truncated_entry = entry.split('.')[0] 334 idx = truncated_entry.index('r') 335 rep_names.append(truncated_entry[:idx] + '|' + truncated_entry[idx:]) 336 337 Ysum = [] 338 339 configlist = [] 340 r_start_index = [] 341 r_stop_index = [] 342 343 for rep in range(replica): 344 345 with open(path + '/' + ls[rep], 'rb') as fp: 346 t = fp.read(12) 347 header = struct.unpack('iii', t) 348 if rep == 0: 349 dn = header[0] 350 nn = header[1] 351 tmax = header[2] 352 elif dn != header[0] or nn != header[1] or tmax != header[2]: 353 raise Exception('Replica parameters do not match.') 354 355 t = fp.read(8) 356 if rep == 0: 357 eps = struct.unpack('d', t)[0] 358 print('Step size:', eps, ', Maximal t value:', dn * (nn) * eps) 359 elif eps != struct.unpack('d', t)[0]: 360 raise Exception('Values for eps do not match among replica.') 361 362 Ysl = [] 363 364 configlist.append([]) 365 while 0 < 1: 366 t = fp.read(4) 367 if(len(t) < 4): 368 break 369 nc = struct.unpack('i', t)[0] 370 configlist[-1].append(nc) 371 372 t = fp.read(8 * tmax * (nn + 1)) 373 if kwargs.get('plaquette'): 374 if nc % dtr_read == 0: 375 Ysl.append(struct.unpack('d' * tmax * (nn + 1), t)) 376 t = fp.read(8 * tmax * (nn + 1)) 377 if not kwargs.get('plaquette'): 378 if nc % dtr_read == 0: 379 Ysl.append(struct.unpack('d' * tmax * (nn + 1), t)) 380 t = fp.read(8 * tmax * (nn + 1)) 381 382 Ysum.append([]) 383 for i, item in enumerate(Ysl): 384 Ysum[-1].append([np.mean(item[current + xmin: 385 current + tmax - xmin]) 386 for current in range(0, len(item), tmax)]) 387 388 diffmeas = configlist[-1][-1] - configlist[-1][-2] 389 configlist[-1] = [item // diffmeas for item in configlist[-1]] 390 if kwargs.get('assume_thermalization', True) and configlist[-1][0] > 1: 391 warnings.warn('Assume thermalization and that the first measurement belongs to the first config.') 392 offset = configlist[-1][0] - 1 393 configlist[-1] = [item - offset for item in configlist[-1]] 394 395 if r_start[rep] is None: 396 r_start_index.append(0) 397 else: 398 try: 399 r_start_index.append(configlist[-1].index(r_start[rep])) 400 except ValueError: 401 raise Exception('Config %d not in file with range [%d, %d]' % ( 402 r_start[rep], configlist[-1][0], configlist[-1][-1])) from None 403 404 if r_stop[rep] is None: 405 r_stop_index.append(len(configlist[-1]) - 1) 406 else: 407 try: 408 r_stop_index.append(configlist[-1].index(r_stop[rep])) 409 except ValueError: 410 raise Exception('Config %d not in file with range [%d, %d]' % ( 411 r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None 412 413 if np.any([len(np.unique(np.diff(cl))) != 1 for cl in configlist]): 414 raise Exception('Irregular spaced data in input file!', [len(np.unique(np.diff(cl))) for cl in configlist]) 415 stepsizes = [list(np.unique(np.diff(cl)))[0] for cl in configlist] 416 if np.any([step != 1 for step in stepsizes]): 417 warnings.warn('Stepsize between configurations is greater than one!' + str(stepsizes), RuntimeWarning) 418 419 idl = [range(configlist[rep][r_start_index[rep]], configlist[rep][r_stop_index[rep]] + 1, r_step) for rep in range(replica)] 420 t2E_dict = {} 421 for n in range(nn + 1): 422 samples = [] 423 for nrep, rep in enumerate(Ysum): 424 samples.append([]) 425 for cnfg in rep: 426 samples[-1].append(cnfg[n]) 427 samples[-1] = samples[-1][r_start_index[nrep]:r_stop_index[nrep] + 1][::r_step] 428 new_obs = Obs(samples, rep_names, idl=idl) 429 t2E_dict[n * dn * eps] = (n * dn * eps) ** 2 * new_obs / (spatial_extent ** 3) - 0.3 430 431 zero_crossing = np.argmax(np.array( 432 [o.value for o in t2E_dict.values()]) > 0.0) 433 434 x = list(t2E_dict.keys())[zero_crossing - fit_range: 435 zero_crossing + fit_range] 436 y = list(t2E_dict.values())[zero_crossing - fit_range: 437 zero_crossing + fit_range] 438 [o.gamma_method() for o in y] 439 440 fit_result = fit_lin(x, y) 441 442 if kwargs.get('plot_fit'): 443 plt.figure() 444 gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1], wspace=0.0, hspace=0.0) 445 ax0 = plt.subplot(gs[0]) 446 xmore = list(t2E_dict.keys())[zero_crossing - fit_range - 2: zero_crossing + fit_range + 2] 447 ymore = list(t2E_dict.values())[zero_crossing - fit_range - 2: zero_crossing + fit_range + 2] 448 [o.gamma_method() for o in ymore] 449 ax0.errorbar(xmore, [yi.value for yi in ymore], yerr=[yi.dvalue for yi in ymore], fmt='x') 450 xplot = np.linspace(np.min(x), np.max(x)) 451 yplot = [fit_result[0] + fit_result[1] * xi for xi in xplot] 452 [yi.gamma_method() for yi in yplot] 453 ax0.fill_between(xplot, y1=[yi.value - yi.dvalue for yi in yplot], y2=[yi.value + yi.dvalue for yi in yplot]) 454 retval = (-fit_result[0] / fit_result[1]) 455 retval.gamma_method() 456 ylim = ax0.get_ylim() 457 ax0.fill_betweenx(ylim, x1=retval.value - retval.dvalue, x2=retval.value + retval.dvalue, color='gray', alpha=0.4) 458 ax0.set_ylim(ylim) 459 ax0.set_ylabel(r'$t^2 \langle E(t) \rangle - 0.3 $') 460 xlim = ax0.get_xlim() 461 462 fit_res = [fit_result[0] + fit_result[1] * xi for xi in x] 463 residuals = (np.asarray([o.value for o in y]) - [o.value for o in fit_res]) / np.asarray([o.dvalue for o in y]) 464 ax1 = plt.subplot(gs[1]) 465 ax1.plot(x, residuals, 'ko', ls='none', markersize=5) 466 ax1.tick_params(direction='out') 467 ax1.tick_params(axis="x", bottom=True, top=True, labelbottom=True) 468 ax1.axhline(y=0.0, ls='--', color='k') 469 ax1.fill_between(xlim, -1.0, 1.0, alpha=0.1, facecolor='k') 470 ax1.set_xlim(xlim) 471 ax1.set_ylabel('Residuals') 472 ax1.set_xlabel(r'$t/a^2$') 473 474 plt.draw() 475 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
View Source
523def read_qtop(path, prefix, c, dtr_cnfg=1, version="openQCD", **kwargs): 524 """Read the topologial charge based on openQCD gradient flow measurements. 525 526 Parameters 527 ---------- 528 path : str 529 path of the measurement files 530 prefix : str 531 prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat. 532 Ignored if file names are passed explicitly via keyword files. 533 c : double 534 Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L. 535 dtr_cnfg : int 536 (optional) parameter that specifies the number of measurements 537 between two configs. 538 If it is not set, the distance between two measurements 539 in the file is assumed to be the distance between two configurations. 540 steps : int 541 (optional) Distance between two configurations in units of trajectories / 542 cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given 543 version : str 544 Either openQCD or sfqcd, depending on the data. 545 L : int 546 spatial length of the lattice in L/a. 547 HAS to be set if version != sfqcd, since openQCD does not provide 548 this in the header 549 r_start : list 550 list which contains the first config to be read for each replicum. 551 r_stop : list 552 list which contains the last config to be read for each replicum. 553 files : list 554 specify the exact files that need to be read 555 from path, practical if e.g. only one replicum is needed 556 names : list 557 Alternative labeling for replicas/ensembles. 558 Has to have the appropriate length. 559 Zeuthen_flow : bool 560 (optional) If True, the Zeuthen flow is used for Qtop. Only possible 561 for version=='sfqcd' If False, the Wilson flow is used. 562 integer_charge : bool 563 If True, the charge is rounded towards the nearest integer on each config. 564 """ 565 known_versions = ["openQCD", "sfqcd"] 566 567 if version not in known_versions: 568 raise Exception("Unknown openQCD version.") 569 if "steps" in kwargs: 570 steps = kwargs.get("steps") 571 if version == "sfqcd": 572 if "L" in kwargs: 573 supposed_L = kwargs.get("L") 574 else: 575 supposed_L = None 576 postfix = ".gfms.dat" 577 else: 578 if "L" not in kwargs: 579 raise Exception("This version of openQCD needs you to provide the spatial length of the lattice as parameter 'L'.") 580 else: 581 L = kwargs.get("L") 582 postfix = ".ms.dat" 583 584 if "files" in kwargs: 585 files = kwargs.get("files") 586 postfix = '' 587 else: 588 found = [] 589 files = [] 590 for (dirpath, dirnames, filenames) in os.walk(path + "/"): 591 found.extend(filenames) 592 break 593 for f in found: 594 if fnmatch.fnmatch(f, prefix + "*" + postfix): 595 files.append(f) 596 597 if 'r_start' in kwargs: 598 r_start = kwargs.get('r_start') 599 if len(r_start) != len(files): 600 raise Exception('r_start does not match number of replicas') 601 r_start = [o if o else None for o in r_start] 602 else: 603 r_start = [None] * len(files) 604 605 if 'r_stop' in kwargs: 606 r_stop = kwargs.get('r_stop') 607 if len(r_stop) != len(files): 608 raise Exception('r_stop does not match number of replicas') 609 else: 610 r_stop = [None] * len(files) 611 rep_names = [] 612 613 zeuthen = kwargs.get('Zeuthen_flow', False) 614 if zeuthen and version not in ['sfqcd']: 615 raise Exception('Zeuthen flow can only be used for version==sfqcd') 616 617 r_start_index = [] 618 r_stop_index = [] 619 deltas = [] 620 configlist = [] 621 for rep, file in enumerate(files): 622 with open(path + "/" + file, "rb") as fp: 623 624 Q = [] 625 traj_list = [] 626 if version in ['sfqcd']: 627 if zeuthen: 628 obspos = 0 629 else: 630 obspos = 8 631 t = fp.read(12) 632 header = struct.unpack('<iii', t) 633 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) 634 ncs = header[1] # number of different values for c in t_flow=1/8 c² L² -> measurements done for ncs c's 635 tmax = header[2] # lattice T/a 636 637 t = fp.read(12) 638 Ls = struct.unpack('<iii', t) 639 if(Ls[0] == Ls[1] and Ls[1] == Ls[2]): 640 L = Ls[0] 641 if not (supposed_L == L) and supposed_L: 642 raise Exception("It seems the length given in the header and by you contradict each other") 643 else: 644 raise Exception("Found more than one spatial length in header!") 645 646 t = fp.read(16) 647 header2 = struct.unpack('<dd', t) 648 tol = header2[0] 649 cmax = header2[1] # highest value of c used 650 651 if c > cmax: 652 raise Exception('Flow has been determined between c=0 and c=%lf with tolerance %lf' % (cmax, tol)) 653 654 if(zthfl == 2): 655 nfl = 2 # number of flows 656 else: 657 nfl = 1 658 iobs = 8 * nfl # number of flow observables calculated 659 660 while 0 < 1: 661 t = fp.read(4) 662 if(len(t) < 4): 663 break 664 traj_list.append(struct.unpack('i', t)[0]) # trajectory number when measurement was done 665 666 for j in range(ncs + 1): 667 for i in range(iobs): 668 t = fp.read(8 * tmax) 669 if (i == obspos): # determines the flow observable -> i=0 <-> Zeuthen flow 670 Q.append(struct.unpack('d' * tmax, t)) 671 672 else: 673 t = fp.read(12) 674 header = struct.unpack('<iii', t) 675 # step size in integration steps "dnms" 676 dn = header[0] 677 # number of measurements, so "ntot"/dn 678 nn = header[1] 679 # lattice T/a 680 tmax = header[2] 681 682 t = fp.read(8) 683 eps = struct.unpack('d', t)[0] 684 685 while 0 < 1: 686 t = fp.read(4) 687 if(len(t) < 4): 688 break 689 traj_list.append(struct.unpack('i', t)[0]) 690 # Wsl 691 t = fp.read(8 * tmax * (nn + 1)) 692 # Ysl 693 t = fp.read(8 * tmax * (nn + 1)) 694 # Qsl, which is asked for in this method 695 t = fp.read(8 * tmax * (nn + 1)) 696 # unpack the array of Qtops, 697 # on each timeslice t=0,...,tmax-1 and the 698 # measurement number in = 0...nn (see README.qcd1) 699 tmpd = struct.unpack('d' * tmax * (nn + 1), t) 700 Q.append(tmpd) 701 702 if len(np.unique(np.diff(traj_list))) != 1: 703 raise Exception("Irregularities in stepsize found") 704 else: 705 if 'steps' in kwargs: 706 if steps != traj_list[1] - traj_list[0]: 707 raise Exception("steps and the found stepsize are not the same") 708 else: 709 steps = traj_list[1] - traj_list[0] 710 711 configlist.append([tr // steps // dtr_cnfg for tr in traj_list]) 712 if configlist[-1][0] > 1: 713 offset = configlist[-1][0] - 1 714 warnings.warn('Assume thermalization and that the first measurement belongs to the first config. Offset = %d configs (%d trajectories / cycles)' % ( 715 offset, offset * steps)) 716 configlist[-1] = [item - offset for item in configlist[-1]] 717 718 if r_start[rep] is None: 719 r_start_index.append(0) 720 else: 721 try: 722 r_start_index.append(configlist[-1].index(r_start[rep])) 723 except ValueError: 724 raise Exception('Config %d not in file with range [%d, %d]' % ( 725 r_start[rep], configlist[-1][0], configlist[-1][-1])) from None 726 727 if r_stop[rep] is None: 728 r_stop_index.append(len(configlist[-1]) - 1) 729 else: 730 try: 731 r_stop_index.append(configlist[-1].index(r_stop[rep])) 732 except ValueError: 733 raise Exception('Config %d not in file with range [%d, %d]' % ( 734 r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None 735 736 if version in ['sfqcd']: 737 cstepsize = cmax / ncs 738 index_aim = round(c / cstepsize) 739 else: 740 t_aim = (c * L) ** 2 / 8 741 index_aim = round(t_aim / eps / dn) 742 743 Q_sum = [] 744 for i, item in enumerate(Q): 745 Q_sum.append([sum(item[current:current + tmax]) 746 for current in range(0, len(item), tmax)]) 747 Q_top = [] 748 if version in ['sfqcd']: 749 for i in range(len(Q_sum) // (ncs + 1)): 750 Q_top.append(Q_sum[i * (ncs + 1) + index_aim][0]) 751 else: 752 for i in range(len(Q) // dtr_cnfg): 753 Q_top.append(Q_sum[dtr_cnfg * i][index_aim]) 754 if len(Q_top) != len(traj_list) // dtr_cnfg: 755 raise Exception("qtops and traj_list dont have the same length") 756 757 if kwargs.get('integer_charge', False): 758 Q_top = [round(q) for q in Q_top] 759 760 truncated_file = file[:-len(postfix)] 761 762 if "names" not in kwargs: 763 try: 764 idx = truncated_file.index('r') 765 except Exception: 766 if "names" not in kwargs: 767 raise Exception("Automatic recognition of replicum failed, please enter the key word 'names'.") 768 ens_name = truncated_file[:idx] 769 rep_names.append(ens_name + '|' + truncated_file[idx:]) 770 else: 771 names = kwargs.get("names") 772 rep_names = names 773 deltas.append(Q_top) 774 775 idl = [range(int(configlist[rep][r_start_index[rep]]), int(configlist[rep][r_stop_index[rep]]) + 1, 1) for rep in range(len(deltas))] 776 deltas = [deltas[nrep][r_start_index[nrep]:r_stop_index[nrep] + 1] for nrep in range(len(deltas))] 777 result = Obs(deltas, rep_names, idl=idl) 778 return result
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
- 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.
View Source
781def qtop_projection(qtop, target=0): 782 """Returns the projection to the topological charge sector defined by target. 783 784 Parameters 785 ---------- 786 path : Obs 787 Topological charge. 788 target : int 789 Specifies the topological sector to be reweighted to (default 0) 790 """ 791 if qtop.reweighted: 792 raise Exception('You can not use a reweighted observable for reweighting!') 793 794 proj_qtop = [] 795 for n in qtop.deltas: 796 proj_qtop.append(np.array([1 if round(qtop.value + q) == target else 0 for q in qtop.deltas[n]])) 797 798 reto = Obs(proj_qtop, qtop.names, idl=[qtop.idl[name] for name in qtop.names]) 799 reto.is_merged = qtop.is_merged 800 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)
View Source
803def read_qtop_sector(path, prefix, c, target=0, **kwargs): 804 """Constructs reweighting factors to a specified topological sector. 805 806 Parameters 807 ---------- 808 path : str 809 path of the measurement files 810 prefix : str 811 prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat 812 c : double 813 Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L 814 target : int 815 Specifies the topological sector to be reweighted to (default 0) 816 dtr_cnfg : int 817 (optional) parameter that specifies the number of trajectories 818 between two configs. 819 if it is not set, the distance between two measurements 820 in the file is assumed to be the distance between two configurations. 821 steps : int 822 (optional) Distance between two configurations in units of trajectories / 823 cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given 824 version : str 825 version string of the openQCD (sfqcd) version used to create 826 the ensemble. Default is 2.0. May also be set to sfqcd. 827 L : int 828 spatial length of the lattice in L/a. 829 HAS to be set if version != sfqcd, since openQCD does not provide 830 this in the header 831 r_start : list 832 offset of the first ensemble, making it easier to match 833 later on with other Obs 834 r_stop : list 835 last configurations that need to be read (per replicum) 836 files : list 837 specify the exact files that need to be read 838 from path, practical if e.g. only one replicum is needed 839 names : list 840 Alternative labeling for replicas/ensembles. 841 Has to have the appropriate length 842 Zeuthen_flow : bool 843 (optional) If True, the Zeuthen flow is used for Qtop. Only possible 844 for version=='sfqcd' If False, the Wilson flow is used. 845 """ 846 847 if not isinstance(target, int): 848 raise Exception("'target' has to be an integer.") 849 850 kwargs['integer_charge'] = True 851 qtop = read_qtop(path, prefix, c, **kwargs) 852 853 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.