From ff738de6af1e76d299de096b40aa9661c934ea4e Mon Sep 17 00:00:00 2001 From: Simon Kuberski Date: Fri, 4 Feb 2022 09:49:17 +0100 Subject: [PATCH 1/2] Ensure correctness of idl in RWFs and t0 --- pyerrors/input/openQCD.py | 90 ++++++++++++++++++++++++++++++--------- 1 file changed, 70 insertions(+), 20 deletions(-) diff --git a/pyerrors/input/openQCD.py b/pyerrors/input/openQCD.py index 4d2eb5fa..91311672 100644 --- a/pyerrors/input/openQCD.py +++ b/pyerrors/input/openQCD.py @@ -3,6 +3,7 @@ import fnmatch import re import struct import numpy as np # Thinly-wrapped numpy +import warnings import matplotlib.pyplot as plt from matplotlib import gridspec from ..obs import Obs @@ -69,8 +70,7 @@ def read_rwms(path, prefix, version='2.0', names=None, **kwargs): r_start = kwargs.get('r_start') if len(r_start) != replica: raise Exception('r_start does not match number of replicas') - # Adjust Configuration numbering to python index - r_start = [o - 1 if o else None for o in r_start] + r_start = [o if o else None for o in r_start] else: r_start = [None] * replica @@ -105,6 +105,10 @@ def read_rwms(path, prefix, version='2.0', names=None, **kwargs): deltas = [] + configlist = [] + r_start_index = [] + r_stop_index = [] + for rep in range(replica): tmp_array = [] with open(path + '/' + ls[rep], 'rb') as fp: @@ -141,15 +145,13 @@ def read_rwms(path, prefix, version='2.0', names=None, **kwargs): if not struct.unpack('i', fp.read(4))[0] == 0: print('something is wrong!') - if r_start[rep] is None: - r_start[rep] = 0 - + configlist.append([]) while 0 < 1: t = fp.read(4) if len(t) < 4: break - if print_err: - config_no = struct.unpack('i', t) + config_no = struct.unpack('i', t)[0] + configlist[-1].append(config_no) for i in range(nrw): if(version == '2.0'): tmpd = _read_array_openQCD2(fp) @@ -179,14 +181,37 @@ def read_rwms(path, prefix, version='2.0', names=None, **kwargs): print('Sources:', np.exp(-np.asarray(tmp_rw))) print('Partial factor:', tmp_nfct) tmp_array[i].append(tmp_nfct) + + if r_start[rep] is None: + r_start_index.append(0) + else: + try: + r_start_index.append(configlist[-1].index(r_start[rep])) + except ValueError: + raise Exception('Config %d not in file with range [%d, %d]' % ( + r_start[rep], configlist[-1][0], configlist[-1][-1])) from None + if r_stop[rep] is None: - r_stop[rep] = len(tmp_array[0]) + r_stop_index.append(len(configlist[-1]) - 1) + else: + try: + r_stop_index.append(configlist[-1].index(r_stop[rep])) + except ValueError: + raise Exception('Config %d not in file with range [%d, %d]' % ( + r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None + for k in range(nrw): - deltas[k].append(tmp_array[k][r_start[rep]:r_stop[rep]][::r_step]) + deltas[k].append(tmp_array[k][r_start_index[rep]:r_stop_index[rep]][::r_step]) + + if np.any([len(np.unique(np.diff(cl))) != 1 for cl in configlist]): + raise Exception('Irregular spaced data in input file!', [len(np.unique(np.diff(cl))) for cl in configlist]) + stepsizes = [list(np.unique(np.diff(cl)))[0] for cl in configlist] + if np.any([step != 1 for step in stepsizes]): + warnings.warn('Stepsize between configurations is greater than one!' + str(stepsizes), RuntimeWarning) print(',', nrw, 'reweighting factors with', nsrc, 'sources') result = [] - idl = [range(r_start[rep] + 1, r_stop[rep] + 1, r_step) for rep in range(replica)] + idl = [range(configlist[rep][r_start_index[rep]], configlist[rep][r_stop_index[rep]], r_step) for rep in range(replica)] for t in range(nrw): result.append(Obs(deltas[t], rep_names, idl=idl)) return result @@ -262,7 +287,7 @@ def extract_t0(path, prefix, dtr_read, xmin, r_start = kwargs.get('r_start') if len(r_start) != replica: raise Exception('r_start does not match number of replicas') - r_start = [o - 1 if o else None for o in r_start] + r_start = [o if o else None for o in r_start] else: r_start = [None] * replica @@ -291,6 +316,10 @@ def extract_t0(path, prefix, dtr_read, xmin, Ysum = [] + configlist = [] + r_start_index = [] + r_stop_index = [] + for rep in range(replica): with open(path + '/' + ls[rep], 'rb') as fp: @@ -312,15 +341,13 @@ def extract_t0(path, prefix, dtr_read, xmin, Ysl = [] - if r_start[rep] is None: - r_start[rep] = 0 - - cfgcount = -1 + configlist.append([]) while 0 < 1: t = fp.read(4) if(len(t) < 4): break nc = struct.unpack('i', t)[0] + configlist[-1].append(nc) t = fp.read(8 * tmax * (nn + 1)) if kwargs.get('plaquette'): @@ -332,16 +359,39 @@ def extract_t0(path, prefix, dtr_read, xmin, Ysl.append(struct.unpack('d' * tmax * (nn + 1), t)) t = fp.read(8 * tmax * (nn + 1)) - cfgcount += 1 - if r_stop[rep] is None: - r_stop[rep] = cfgcount 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)]) - idl = [range(r_start[rep] + 1, r_stop[rep] + 1, r_step) for rep in range(len(r_start))] + diffmeas = configlist[-1][-1] - configlist[-1][-2] + configlist[-1] = [item // diffmeas for item in configlist[-1]] + if r_start[rep] is None: + r_start_index.append(0) + else: + try: + r_start_index.append(configlist[-1].index(r_start[rep])) + except ValueError: + raise Exception('Config %d not in file with range [%d, %d]' % ( + r_start[rep], configlist[-1][0], configlist[-1][-1])) from None + + if r_stop[rep] is None: + r_stop_index.append(len(configlist[-1]) - 1) + else: + try: + r_stop_index.append(configlist[-1].index(r_stop[rep])) + except ValueError: + raise Exception('Config %d not in file with range [%d, %d]' % ( + r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None + + if np.any([len(np.unique(np.diff(cl))) != 1 for cl in configlist]): + raise Exception('Irregular spaced data in input file!', [len(np.unique(np.diff(cl))) for cl in configlist]) + stepsizes = [list(np.unique(np.diff(cl)))[0] for cl in configlist] + if np.any([step != 1 for step in stepsizes]): + warnings.warn('Stepsize between configurations is greater than one!' + str(stepsizes), RuntimeWarning) + + idl = [range(configlist[rep][r_start_index[rep]], configlist[rep][r_stop_index[rep]], r_step) for rep in range(replica)] t2E_dict = {} for n in range(nn + 1): samples = [] @@ -349,7 +399,7 @@ def extract_t0(path, prefix, dtr_read, xmin, samples.append([]) for cnfg in rep: samples[-1].append(cnfg[n]) - samples[-1] = samples[-1][r_start[nrep]:r_stop[nrep]][::r_step] + samples[-1] = samples[-1][r_start_index[nrep]:r_stop_index[nrep]][::r_step] new_obs = Obs(samples, rep_names, idl=idl) t2E_dict[n * dn * eps] = (n * dn * eps) ** 2 * new_obs / (spatial_extent ** 3) - 0.3 From 65568c84a41d8e25dcb30353ee72d213c7d06f60 Mon Sep 17 00:00:00 2001 From: Simon Kuberski Date: Fri, 4 Feb 2022 10:15:43 +0100 Subject: [PATCH 2/2] Handling thermalization correctly in t0 --- pyerrors/input/openQCD.py | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/pyerrors/input/openQCD.py b/pyerrors/input/openQCD.py index 91311672..3c50a3d2 100644 --- a/pyerrors/input/openQCD.py +++ b/pyerrors/input/openQCD.py @@ -226,7 +226,12 @@ def extract_t0(path, prefix, dtr_read, 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. + Only works with openQCD + + 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 ---------- @@ -263,6 +268,11 @@ def extract_t0(path, prefix, dtr_read, xmin, 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 """ ls = [] @@ -367,6 +377,11 @@ def extract_t0(path, prefix, dtr_read, xmin, diffmeas = configlist[-1][-1] - configlist[-1][-2] configlist[-1] = [item // diffmeas for item in configlist[-1]] + if kwargs.get('assume_thermalization', True) and configlist[-1][0] > 1: + warnings.warn('Assume thermalization and that the first measurement belongs to the first config.') + offset = configlist[-1][0] - 1 + configlist[-1] = [item - offset for item in configlist[-1]] + if r_start[rep] is None: r_start_index.append(0) else: