diff --git a/pyerrors/input/openQCD.py b/pyerrors/input/openQCD.py index 5c44fd2f..60ed64c7 100644 --- a/pyerrors/input/openQCD.py +++ b/pyerrors/input/openQCD.py @@ -8,7 +8,6 @@ import struct import numpy as np # Thinly-wrapped numpy from ..obs import Obs from ..fits import fit_lin -from . import utils def read_rwms(path, prefix, version='2.0', names=None, **kwargs): @@ -67,7 +66,8 @@ def read_rwms(path, prefix, version='2.0', names=None, **kwargs): else: r_stop = [None] * replica - print('Read reweighting factors from', prefix[:-1], ',', replica, 'replica', end='') + print('Read reweighting factors from', prefix[:-1], ',', + replica, 'replica', end='') # Adjust replica names to new bookmarking system if names is None: @@ -75,7 +75,8 @@ def read_rwms(path, prefix, version='2.0', names=None, **kwargs): for entry in ls: truncated_entry = entry.split('.')[0] idx = truncated_entry.index('r') - rep_names.append(truncated_entry[:idx] + '|' + truncated_entry[idx:]) + rep_names.append(truncated_entry[:idx] + '|' + + truncated_entry[idx:]) print_err = 0 if 'print_err' in kwargs: @@ -97,8 +98,13 @@ def read_rwms(path, prefix, version='2.0', names=None, **kwargs): for k in range(nrw): deltas.append([]) else: - if ((nrw != struct.unpack('i', t)[0] and (not version == '2.0')) or (nrw != struct.unpack('i', t)[0] / 2 and version == '2.0')): # little weird if-clause due to the /2 operation needed. - raise Exception('Error: different number of reweighting factors for replicum', rep) + # little weird if-clause due to the /2 operation needed. + if ((nrw != struct.unpack('i', t)[0] and + (not version == '2.0')) or + (nrw != struct.unpack('i', t)[0] / 2 and + version == '2.0')): + raise Exception('Error: different number of reweighting\ + factors for replicum', rep) for k in range(nrw): tmp_array.append([]) @@ -109,7 +115,8 @@ def read_rwms(path, prefix, version='2.0', names=None, **kwargs): for i in range(nrw): t = fp.read(4) nfct.append(struct.unpack('i', t)[0]) - # print('nfct: ', nfct) # Hasenbusch factor, 1 for rat reweighting + # print('nfct: ', nfct) # Hasenbusch factor, + # 1 for rat reweighting else: for i in range(nrw): nfct.append(1) @@ -138,8 +145,11 @@ def read_rwms(path, prefix, version='2.0', names=None, **kwargs): for j in range(tmpd['n'][0]): tmp_nfct *= np.mean(np.exp(-np.asarray(tmp_rw[j]))) if print_err: - print(config_no, i, j, np.mean(np.exp(-np.asarray(tmp_rw[j]))), np.std(np.exp(-np.asarray(tmp_rw[j])))) - print('Sources:', np.exp(-np.asarray(tmp_rw[j]))) + print(config_no, i, j, + np.mean(np.exp(-np.asarray(tmp_rw[j]))), + np.std(np.exp(-np.asarray(tmp_rw[j])))) + print('Sources:', + np.exp(-np.asarray(tmp_rw[j]))) print('Partial factor:', tmp_nfct) elif version == '1.6' or version == '1.4': tmp_nfct = 1.0 @@ -149,7 +159,9 @@ def read_rwms(path, prefix, version='2.0', names=None, **kwargs): tmp_rw = struct.unpack('d' * nsrc[i], t) tmp_nfct *= np.mean(np.exp(-np.asarray(tmp_rw))) if print_err: - print(config_no, i, j, np.mean(np.exp(-np.asarray(tmp_rw))), np.std(np.exp(-np.asarray(tmp_rw)))) + print(config_no, i, j, + np.mean(np.exp(-np.asarray(tmp_rw))), + np.std(np.exp(-np.asarray(tmp_rw)))) print('Sources:', np.exp(-np.asarray(tmp_rw))) print('Partial factor:', tmp_nfct) tmp_array[i].append(tmp_nfct) @@ -168,11 +180,14 @@ def read_rwms(path, prefix, version='2.0', names=None, **kwargs): return result -def extract_t0(path, prefix, dtr_read, xmin, spatial_extent, fit_range=5, **kwargs): +def extract_t0(path, prefix, dtr_read, xmin, + spatial_extent, fit_range=5, **kwargs): """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 - 0.3 is fitted with a linear function + It is assumed that all boundary effects have + sufficiently decayed at x0=xmin. + The data around the zero crossing of t^2 - 0.3 + is fitted with a linear function from which the exact root is extracted. Only works with openQCD v 1.2. @@ -183,14 +198,17 @@ def extract_t0(path, prefix, dtr_read, xmin, spatial_extent, fit_range=5, **kwar prefix : str Ensemble prefix dtr_read : int - Determines how many trajectories should be skipped when reading the ms.dat files. + 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. + 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) + 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 @@ -276,7 +294,9 @@ def extract_t0(path, prefix, dtr_read, xmin, spatial_extent, fit_range=5, **kwar Ysum.append([]) for i, item in enumerate(Ysl): - Ysum[-1].append([np.mean(item[current + xmin:current + tmax - xmin]) for current in range(0, len(item), tmax)]) + Ysum[-1].append([np.mean(item[current + xmin: + current + tmax - xmin]) + for current in range(0, len(item), tmax)]) t2E_dict = {} for n in range(nn + 1): @@ -287,12 +307,16 @@ def extract_t0(path, prefix, dtr_read, xmin, spatial_extent, fit_range=5, **kwar samples[-1].append(cnfg[n]) samples[-1] = samples[-1][r_start[nrep]:r_stop[nrep]] new_obs = Obs(samples, [(w.split('.'))[0] for w in ls]) - t2E_dict[n * dn * eps] = (n * dn * eps) ** 2 * new_obs / (spatial_extent ** 3) - 0.3 + t2E_dict[n * dn * eps] = (n * dn * eps) ** 2 * new_obs \ + / (spatial_extent ** 3) - 0.3 - zero_crossing = np.argmax(np.array([o.value for o in t2E_dict.values()]) > 0.0) + zero_crossing = np.argmax(np.array( + [o.value for o in t2E_dict.values()]) > 0.0) - x = list(t2E_dict.keys())[zero_crossing - fit_range: zero_crossing + fit_range] - y = list(t2E_dict.values())[zero_crossing - fit_range: zero_crossing + fit_range] + x = list(t2E_dict.keys())[zero_crossing - fit_range: + zero_crossing + fit_range] + y = list(t2E_dict.values())[zero_crossing - fit_range: + zero_crossing + fit_range] [o.gamma_method() for o in y] fit_result = fit_lin(x, y) @@ -348,7 +372,7 @@ def _read_array_openQCD2(fp): return {'d': d, 'n': n, 'size': size, 'arr': arr} -def read_qtop(path, prefix,c, dtr_cnfg = 1, version = "1.2",**kwargs): +def read_qtop(path, prefix, c, dtr_cnfg=1, version="1.2", **kwargs): """Read qtop format from given folder structure. Parameters @@ -360,144 +384,150 @@ def read_qtop(path, prefix,c, dtr_cnfg = 1, version = "1.2",**kwargs): 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 trajectories between two configs. - if it is not set, the distance between two measurements in the file is assumed to be + (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) (maybe only necessary for openQCD2.0) nt step size, guessed if not given version: str - version string of the openQCD (sfqcd) version used to create the ensemble + version string of the openQCD (sfqcd) version used to create + the ensemble 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 + 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 + 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) - r_meas_start: list - offset of the first measured ensemble, if there is any files: list - specify the exact files that need to be read from path, pratical if e.g. only one replicum is needed + specify the exact files that need to be read + from path, pratical if e.g. only one replicum is needed names: list - Alternative labeling for replicas/ensembles. Has to have the appropriate length + Alternative labeling for replicas/ensembles. + Has to have the appropriate length """ - #one could read L from the header in case of sfQCD - #c = 0.35 - known_versions = ["1.0","1.2","1.4","1.6","2.0", "sfqcd"] - - if not version in known_versions: + # one could read L from the header in case of sfQCD + # c = 0.35 + known_versions = ["1.0", "1.2", "1.4", "1.6", "2.0", "sfqcd"] + + if version not in known_versions: raise Exception("Unknown openQCD version.") - target = 0 if "steps" in kwargs: steps = kwargs.get("steps") - - if 'target' in kwargs: - target = kwargs.get('target') if version == "sfqcd": if "L" in kwargs: - supposed_L = kwargs.get("L") + supposed_L = kwargs.get("L") else: - if not "L" in kwargs: - raise Exception("This version of openQCD needs you to provide the spatial length of the lattice as parameter 'L'.") + if "L" not in kwargs: + raise Exception("This version of openQCD needs you \ + to provide the spatial length of the \ + lattice as parameter 'L'.") else: L = kwargs.get("L") r_start = 1 - r_meas_start = 1 - if "r_meas_start" in kwargs: - r_meas_start = kwargs.get("r_meas_start") if "r_start" in kwargs: r_start = kwargs.get("r_start") if "r_stop" in kwargs: r_stop = kwargs.get("r_stop") - #if one wants to read specific files with this method... + # if one wants to read specific files with this method... if "files" in kwargs: files = kwargs.get("files") else: - #find files in path + # find files in path found = [] files = [] for (dirpath, dirnames, filenames) in os.walk(path+"/"): - #print(filenames) + # print(filenames) found.extend(filenames) break for f in found: if fnmatch.fnmatch(f, prefix+"*"+".ms.dat"): files.append(f) print(files) - #now that we found our files, we dechiffer them... + # now that we found our files, we dechiffer them... rep_names = [] - + deltas = [] idl = [] - for rep,file in enumerate(files): - + for rep, file in enumerate(files): with open(path+"/"+file, "rb") as fp: - #this, for now, is for version 1.2,1.4,1.6 and 2.0, but needs to be tested for the last 3, isncethe doc says its the same - #header + # header t = fp.read(12) header = struct.unpack(' 1: + range(1, len(data_starts))])) > 1: raise Exception("Irregularities in file structure found,\ - not all runs have the same output length") + not all runs have the same output length") # first chunk of data chunk = content[:data_starts[1]] for linenumber, line in enumerate(chunk): @@ -397,7 +398,7 @@ def read_sfcf(path, prefix, name, quarks='.*', noffset=0, wf=0, wf2=0, found_pat += li if re.search(pattern, found_pat): for t, line in \ - enumerate(chunk[start_read:start_read+T]): + enumerate(chunk[start_read:start_read+T]): floats = list(map(float, line.split())) deltas[t][rep][cnfg] = floats[-2:][im] idl.append(rep_idl)