From e2abce203c2d07e85f3bca4f86d52112fb90696a Mon Sep 17 00:00:00 2001 From: Justus Kuhlmann Date: Tue, 11 Aug 2026 15:10:36 +0200 Subject: [PATCH 1/6] add a first working version of a extra plot for the extraction of t0 --- pyerrors/input/misc.py | 99 +++++++++++++++++++++++++++++++++++++++ pyerrors/input/openQCD.py | 15 ++++-- 2 files changed, 109 insertions(+), 5 deletions(-) diff --git a/pyerrors/input/misc.py b/pyerrors/input/misc.py index c7c7f5b2..f868e785 100644 --- a/pyerrors/input/misc.py +++ b/pyerrors/input/misc.py @@ -8,10 +8,109 @@ import numpy as np # Thinly-wrapped numpy from matplotlib import gridspec +from ..correlators import Corr from ..fits import fit_lin from ..obs import Obs +def plot_Ysl(Ysl, nn, dn, eps, tmax, xmin, r_start, r_stop, r_step, names, spatial_extent): + """ + Plot the plateaus needed for the fit of t0. + """ + + Ysl_samples = [] + for rep in range(len(Ysl)): + Ysl_array = [] + for cnfg in range(len(Ysl[rep])): + yarr = [] + for x0 in range(tmax): + yarr.append([]) + for flow in range(nn + 1): + ind = flow*tmax+x0 + yarr[x0].append(Ysl[rep][cnfg][ind]) + Ysl_array.append(yarr) + Ysl_samples.append(Ysl_array) + t2E_arr = [] + + for t_flow in range(nn + 1): + corr_obs =[] + for x0 in range(tmax): + samples = [] + rep_idls = [] + for rep in range(len(Ysl)): + samples.append([]) + rep_idls.append([]) + for cnfg in range(len(Ysl_samples[rep])): + if cnfg>=r_start[rep] and cnfg<=r_stop[rep] and (cnfg-r_start[rep]) % r_step==0: + samples[rep].append(Ysl_samples[rep][cnfg][x0][t_flow]) + rep_idls[rep].append(cnfg+1) + o = Obs(samples, names, rep_idls) + corr_obs.append(o) + E = Corr(corr_obs)/(spatial_extent**3) + t = (dn*eps*t_flow) + t2E = t**2*E + t2E.gm() + t2E_arr.append(t2E) + + ts = [] + t2expE_arr = [] + for t_flow in range(nn + 1): + samples = [] + rep_idls = [] + for rep in range(len(Ysl)): + samples.append([]) + rep_idls.append([]) + for cnfg in range(len(Ysl_samples[rep])): + if cnfg>=r_start[rep] and cnfg<=r_stop[rep] and (cnfg-r_start[rep]) % r_step==0: + myls = [Ysl_samples[rep][cnfg][x0][t_flow] for x0 in range(xmin, tmax - xmin)] + samples[rep].append(np.mean(myls)) + rep_idls[rep].append(cnfg+1) + o = Obs(samples, names, rep_idls) + expE = o/(spatial_extent**3) + t = (dn*eps*t_flow) + ts.append(t) + t2expE = t**2*expE + t2expE.gm() + t2expE_arr.append(t2expE) + + exp_vals = [] + ts = [] + + def compat(val1, val2): + return bool(np.abs(val1.value - val2.value) < np.abs(val1.dvalue + val2.dvalue)) + + # find t2expE element closest to .3 + closest_item = 0 + closest_value = 10 + for item, value in enumerate(t2expE_arr): + if abs(value.value-0.3) < abs(closest_value-0.3): + closest_item = item + closest_value = value + + if closest_value < 0.3: + closest_item += 1 + + range_max = min(len(t2expE_arr)-1, closest_item+5) + range_min = max(0, closest_item-5) + prange = [xmin, tmax-xmin] + print(range_min, range_max) + for i in range(range_min, range_max): + t = dn*eps*i + t2expE_arr[i].gm() + exp_vals.append(t2expE_arr[i]) + ts.append(t) + x,y,yerr = t2E_arr[i].plottable() + plt.errorbar(x, y, yerr, alpha = 0.7, color = f"C{i:02d}", linestyle="None", marker=".") + plt.fill_between(prange, t2expE_arr[i].value - t2expE_arr[i].dvalue, t2expE_arr[i].value + t2expE_arr[i].dvalue, alpha = 0.3, color = f"C{i:02d}") + plt.hlines(t2expE_arr[i], prange[0], prange[1], linestyle = "dashed", colors=f"C{i:02d}", label = r"$t^2\langle E(t)\rangle$") + + plt.ylabel(r"$t^2E(t)$") + plt.xlabel("$x_{0}/a$") + plt.xticks([i * tmax/4 for i in range(5)]) + plt.xlim(0,tmax) + plt.show() + + def fit_t0(t2E_dict, fit_range, plot_fit=False, observable='t0'): """Compute the root of (flow-based) data based on a dictionary that contains the necessary information in key-value pairs a la (flow time: observable at flow time). diff --git a/pyerrors/input/openQCD.py b/pyerrors/input/openQCD.py index 6405a98f..44dd82d0 100644 --- a/pyerrors/input/openQCD.py +++ b/pyerrors/input/openQCD.py @@ -7,7 +7,7 @@ from ..correlators import Corr from ..obs import CObs, Obs -from .misc import fit_t0 +from .misc import fit_t0, plot_Ysl from .utils import sort_names @@ -329,6 +329,7 @@ def _extract_flowed_energy_density(path, prefix, dtr_read, xmin, spatial_extent, idx = truncated_entry.index('r') rep_names.append(truncated_entry[:idx] + '|' + truncated_entry[idx:]) + Ysl = [] Ysum = [] configlist = [] @@ -354,7 +355,7 @@ def _extract_flowed_energy_density(path, prefix, dtr_read, xmin, spatial_extent, elif eps != struct.unpack('d', t)[0]: raise Exception('Values for eps do not match among replica.') - Ysl = [] + Ysl_rep = [] configlist.append([]) while True: @@ -367,15 +368,16 @@ def _extract_flowed_energy_density(path, prefix, dtr_read, xmin, spatial_extent, t = fp.read(8 * tmax * (nn + 1)) if kwargs.get('plaquette'): if nc % dtr_read == 0: - Ysl.append(struct.unpack('d' * tmax * (nn + 1), t)) + Ysl_rep.append(struct.unpack('d' * tmax * (nn + 1), t)) t = fp.read(8 * tmax * (nn + 1)) if not kwargs.get('plaquette'): if nc % dtr_read == 0: - Ysl.append(struct.unpack('d' * tmax * (nn + 1), t)) + Ysl_rep.append(struct.unpack('d' * tmax * (nn + 1), t)) t = fp.read(8 * tmax * (nn + 1)) + Ysl.append(Ysl_rep) Ysum.append([]) - for _i, item in enumerate(Ysl): + for _i, item in enumerate(Ysl_rep): Ysum[-1].append([np.mean(item[current + xmin: current + tmax - xmin]) for current in range(0, len(item), tmax)]) @@ -409,6 +411,9 @@ def _extract_flowed_energy_density(path, prefix, dtr_read, xmin, spatial_extent, f'Config {r_stop[rep]} not in file with range [{configlist[-1][0]}, {configlist[-1][-1]}]' ) from None + if kwargs.get("plot_Ysl", False): + plot_Ysl(Ysl, nn, dn, eps, tmax, xmin, r_start_index, r_stop_index, r_step, rep_names, spatial_extent) + 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 = [next(iter(np.unique(np.diff(cl)))) for cl in configlist] From 9d1d02c70edb95a21895873b539cd2a31fdf6b39 Mon Sep 17 00:00:00 2001 From: Justus Kuhlmann Date: Tue, 11 Aug 2026 15:14:03 +0200 Subject: [PATCH 2/6] fix ylim --- pyerrors/input/misc.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pyerrors/input/misc.py b/pyerrors/input/misc.py index f868e785..5995c81d 100644 --- a/pyerrors/input/misc.py +++ b/pyerrors/input/misc.py @@ -93,7 +93,7 @@ def compat(val1, val2): range_max = min(len(t2expE_arr)-1, closest_item+5) range_min = max(0, closest_item-5) prange = [xmin, tmax-xmin] - print(range_min, range_max) + for i in range(range_min, range_max): t = dn*eps*i t2expE_arr[i].gm() @@ -104,6 +104,7 @@ def compat(val1, val2): plt.fill_between(prange, t2expE_arr[i].value - t2expE_arr[i].dvalue, t2expE_arr[i].value + t2expE_arr[i].dvalue, alpha = 0.3, color = f"C{i:02d}") plt.hlines(t2expE_arr[i], prange[0], prange[1], linestyle = "dashed", colors=f"C{i:02d}", label = r"$t^2\langle E(t)\rangle$") + plt.ylim([t2expE_arr[range_min].value-.03,t2expE_arr[range_max].value+.01]) plt.ylabel(r"$t^2E(t)$") plt.xlabel("$x_{0}/a$") plt.xticks([i * tmax/4 for i in range(5)]) From 1868b2abba91e5786a9ebf115c1890d4768e71d6 Mon Sep 17 00:00:00 2001 From: Justus Kuhlmann Date: Fri, 14 Aug 2026 15:06:39 +0200 Subject: [PATCH 3/6] reuse Ysum, only use Ysl when plotting, plot all plateaus at once --- pyerrors/input/misc.py | 70 ++++++++++++--------------------------- pyerrors/input/openQCD.py | 15 +++++---- 2 files changed, 29 insertions(+), 56 deletions(-) diff --git a/pyerrors/input/misc.py b/pyerrors/input/misc.py index 5995c81d..6fd9ab74 100644 --- a/pyerrors/input/misc.py +++ b/pyerrors/input/misc.py @@ -13,11 +13,7 @@ from ..obs import Obs -def plot_Ysl(Ysl, nn, dn, eps, tmax, xmin, r_start, r_stop, r_step, names, spatial_extent): - """ - Plot the plateaus needed for the fit of t0. - """ - +def _reorder_Ysl(Ysl, tmax, nn): Ysl_samples = [] for rep in range(len(Ysl)): Ysl_array = [] @@ -30,10 +26,23 @@ def plot_Ysl(Ysl, nn, dn, eps, tmax, xmin, r_start, r_stop, r_step, names, spati yarr[x0].append(Ysl[rep][cnfg][ind]) Ysl_array.append(yarr) Ysl_samples.append(Ysl_array) - t2E_arr = [] + return Ysl_samples + + +def plot_Ysl(Ysl, expE_dict, nn, dn, eps, tmax, xmin, r_start, r_stop, r_step, names, spatial_extent): + """ + Plot the plateaus needed for the fit of t0. + """ + + prange = [xmin, tmax-xmin] - for t_flow in range(nn + 1): - corr_obs =[] + Ysl_samples = _reorder_Ysl(Ysl, tmax, nn) + ts = list(expE_dict.keys()) + + t2E_arr = [] + t2expE_arr = [] + for n, t in zip(range(nn + 1), ts, strict=True): + corr_obs = [] for x0 in range(tmax): samples = [] rep_idls = [] @@ -42,59 +51,22 @@ def plot_Ysl(Ysl, nn, dn, eps, tmax, xmin, r_start, r_stop, r_step, names, spati rep_idls.append([]) for cnfg in range(len(Ysl_samples[rep])): if cnfg>=r_start[rep] and cnfg<=r_stop[rep] and (cnfg-r_start[rep]) % r_step==0: - samples[rep].append(Ysl_samples[rep][cnfg][x0][t_flow]) + samples[rep].append(Ysl_samples[rep][cnfg][x0][n]) rep_idls[rep].append(cnfg+1) o = Obs(samples, names, rep_idls) corr_obs.append(o) E = Corr(corr_obs)/(spatial_extent**3) - t = (dn*eps*t_flow) t2E = t**2*E t2E.gm() t2E_arr.append(t2E) - - ts = [] - t2expE_arr = [] - for t_flow in range(nn + 1): - samples = [] - rep_idls = [] - for rep in range(len(Ysl)): - samples.append([]) - rep_idls.append([]) - for cnfg in range(len(Ysl_samples[rep])): - if cnfg>=r_start[rep] and cnfg<=r_stop[rep] and (cnfg-r_start[rep]) % r_step==0: - myls = [Ysl_samples[rep][cnfg][x0][t_flow] for x0 in range(xmin, tmax - xmin)] - samples[rep].append(np.mean(myls)) - rep_idls[rep].append(cnfg+1) - o = Obs(samples, names, rep_idls) - expE = o/(spatial_extent**3) - t = (dn*eps*t_flow) - ts.append(t) - t2expE = t**2*expE + t2expE = t**2*expE_dict[t] t2expE.gm() t2expE_arr.append(t2expE) exp_vals = [] ts = [] - def compat(val1, val2): - return bool(np.abs(val1.value - val2.value) < np.abs(val1.dvalue + val2.dvalue)) - - # find t2expE element closest to .3 - closest_item = 0 - closest_value = 10 - for item, value in enumerate(t2expE_arr): - if abs(value.value-0.3) < abs(closest_value-0.3): - closest_item = item - closest_value = value - - if closest_value < 0.3: - closest_item += 1 - - range_max = min(len(t2expE_arr)-1, closest_item+5) - range_min = max(0, closest_item-5) - prange = [xmin, tmax-xmin] - - for i in range(range_min, range_max): + for i in range(len(t2expE_arr)): t = dn*eps*i t2expE_arr[i].gm() exp_vals.append(t2expE_arr[i]) @@ -104,7 +76,7 @@ def compat(val1, val2): plt.fill_between(prange, t2expE_arr[i].value - t2expE_arr[i].dvalue, t2expE_arr[i].value + t2expE_arr[i].dvalue, alpha = 0.3, color = f"C{i:02d}") plt.hlines(t2expE_arr[i], prange[0], prange[1], linestyle = "dashed", colors=f"C{i:02d}", label = r"$t^2\langle E(t)\rangle$") - plt.ylim([t2expE_arr[range_min].value-.03,t2expE_arr[range_max].value+.01]) + plt.ylim([t2expE_arr[0].value-.03,t2expE_arr[-1].value+.01]) plt.ylabel(r"$t^2E(t)$") plt.xlabel("$x_{0}/a$") plt.xticks([i * tmax/4 for i in range(5)]) diff --git a/pyerrors/input/openQCD.py b/pyerrors/input/openQCD.py index 44dd82d0..8a191340 100644 --- a/pyerrors/input/openQCD.py +++ b/pyerrors/input/openQCD.py @@ -374,7 +374,8 @@ def _extract_flowed_energy_density(path, prefix, dtr_read, xmin, spatial_extent, if nc % dtr_read == 0: Ysl_rep.append(struct.unpack('d' * tmax * (nn + 1), t)) t = fp.read(8 * tmax * (nn + 1)) - Ysl.append(Ysl_rep) + if kwargs.get("plot_Ysl", False): + Ysl.append(Ysl_rep) Ysum.append([]) for _i, item in enumerate(Ysl_rep): @@ -411,9 +412,6 @@ def _extract_flowed_energy_density(path, prefix, dtr_read, xmin, spatial_extent, f'Config {r_stop[rep]} not in file with range [{configlist[-1][0]}, {configlist[-1][-1]}]' ) from None - if kwargs.get("plot_Ysl", False): - plot_Ysl(Ysl, nn, dn, eps, tmax, xmin, r_start_index, r_stop_index, r_step, rep_names, spatial_extent) - 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 = [next(iter(np.unique(np.diff(cl)))) for cl in configlist] @@ -421,7 +419,7 @@ def _extract_flowed_energy_density(path, prefix, dtr_read, xmin, spatial_extent, warnings.warn('Stepsize between configurations is greater than one!' + str(stepsizes), RuntimeWarning, stacklevel=2) idl = [range(configlist[rep][r_start_index[rep]], configlist[rep][r_stop_index[rep]] + 1, r_step) for rep in range(replica)] - E_dict = {} + expE_dict = {} for n in range(nn + 1): samples = [] for nrep, rep in enumerate(Ysum): @@ -430,9 +428,12 @@ def _extract_flowed_energy_density(path, prefix, dtr_read, xmin, spatial_extent, samples[-1].append(cnfg[n]) samples[-1] = samples[-1][r_start_index[nrep]:r_stop_index[nrep] + 1][::r_step] new_obs = Obs(samples, rep_names, idl=idl) - E_dict[n * dn * eps] = new_obs / (spatial_extent ** 3) + expE_dict[n * dn * eps] = new_obs / (spatial_extent ** 3) + + if kwargs.get("plot_Ysl", False): + plot_Ysl(Ysl, expE_dict, nn, dn, eps, tmax, xmin, r_start_index, r_stop_index, r_step, rep_names, spatial_extent) - return E_dict + return expE_dict def extract_t0(path, prefix, dtr_read, xmin, spatial_extent, fit_range=5, postfix='ms', c=0.3, **kwargs): From 3ecd33e8d9db5abebe98b1a2599d2af12cf393b8 Mon Sep 17 00:00:00 2001 From: Justus Kuhlmann Date: Fri, 14 Aug 2026 15:40:02 +0200 Subject: [PATCH 4/6] change internal function signature, evaluate zero_crossing for plot --- pyerrors/input/misc.py | 8 +++++--- pyerrors/input/openQCD.py | 8 ++++---- 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/pyerrors/input/misc.py b/pyerrors/input/misc.py index 6fd9ab74..08e497da 100644 --- a/pyerrors/input/misc.py +++ b/pyerrors/input/misc.py @@ -29,7 +29,7 @@ def _reorder_Ysl(Ysl, tmax, nn): return Ysl_samples -def plot_Ysl(Ysl, expE_dict, nn, dn, eps, tmax, xmin, r_start, r_stop, r_step, names, spatial_extent): +def plot_Ysl(Ysl, expE_dict, nn, dn, eps, tmax, xmin, r_start, r_stop, r_step, names, spatial_extent, fit_range, c): """ Plot the plateaus needed for the fit of t0. """ @@ -66,7 +66,9 @@ def plot_Ysl(Ysl, expE_dict, nn, dn, eps, tmax, xmin, r_start, r_stop, r_step, n exp_vals = [] ts = [] - for i in range(len(t2expE_arr)): + zero_crossing = np.argmax(np.array([o.value - c for o in t2expE_arr]) > 0.0) + + for i in range(zero_crossing-fit_range, zero_crossing+fit_range): t = dn*eps*i t2expE_arr[i].gm() exp_vals.append(t2expE_arr[i]) @@ -76,7 +78,7 @@ def plot_Ysl(Ysl, expE_dict, nn, dn, eps, tmax, xmin, r_start, r_stop, r_step, n plt.fill_between(prange, t2expE_arr[i].value - t2expE_arr[i].dvalue, t2expE_arr[i].value + t2expE_arr[i].dvalue, alpha = 0.3, color = f"C{i:02d}") plt.hlines(t2expE_arr[i], prange[0], prange[1], linestyle = "dashed", colors=f"C{i:02d}", label = r"$t^2\langle E(t)\rangle$") - plt.ylim([t2expE_arr[0].value-.03,t2expE_arr[-1].value+.01]) + plt.ylim([t2expE_arr[zero_crossing-fit_range].value-.03,t2expE_arr[zero_crossing+fit_range].value+.01]) plt.ylabel(r"$t^2E(t)$") plt.xlabel("$x_{0}/a$") plt.xticks([i * tmax/4 for i in range(5)]) diff --git a/pyerrors/input/openQCD.py b/pyerrors/input/openQCD.py index 8a191340..a5504094 100644 --- a/pyerrors/input/openQCD.py +++ b/pyerrors/input/openQCD.py @@ -232,7 +232,7 @@ def read_rwms(path, prefix, version='2.0', names=None, **kwargs): return result -def _extract_flowed_energy_density(path, prefix, dtr_read, xmin, spatial_extent, postfix='ms', **kwargs): +def _extract_flowed_energy_density(path, prefix, dtr_read, xmin, spatial_extent, fit_range, c, postfix='ms', **kwargs): """Extract a dictionary with the flowed Yang-Mills action density from given .ms.dat files. Returns a dictionary with Obs as values and flow times as keys. @@ -431,7 +431,7 @@ def _extract_flowed_energy_density(path, prefix, dtr_read, xmin, spatial_extent, expE_dict[n * dn * eps] = new_obs / (spatial_extent ** 3) if kwargs.get("plot_Ysl", False): - plot_Ysl(Ysl, expE_dict, nn, dn, eps, tmax, xmin, r_start_index, r_stop_index, r_step, rep_names, spatial_extent) + plot_Ysl(Ysl, expE_dict, nn, dn, eps, tmax, xmin, r_start_index, r_stop_index, r_step, rep_names, spatial_extent, fit_range, c) return expE_dict @@ -502,7 +502,7 @@ def extract_t0(path, prefix, dtr_read, xmin, spatial_extent, fit_range=5, postfi Extracted t0 """ - E_dict = _extract_flowed_energy_density(path, prefix, dtr_read, xmin, spatial_extent, postfix, **kwargs) + E_dict = _extract_flowed_energy_density(path, prefix, dtr_read, xmin, spatial_extent, fit_range, c, postfix, **kwargs) t2E_dict = {} for t in sorted(E_dict.keys()): t2E_dict[t] = t ** 2 * E_dict[t] - c @@ -575,7 +575,7 @@ def extract_w0(path, prefix, dtr_read, xmin, spatial_extent, fit_range=5, postfi Extracted w0 """ - E_dict = _extract_flowed_energy_density(path, prefix, dtr_read, xmin, spatial_extent, postfix, **kwargs) + E_dict = _extract_flowed_energy_density(path, prefix, dtr_read, xmin, spatial_extent, fit_range, c, postfix, **kwargs) ftimes = sorted(E_dict.keys()) From 3b9753a138121665d6df83ae909c0b1f2a20ba94 Mon Sep 17 00:00:00 2001 From: Justus Kuhlmann Date: Thu, 20 Aug 2026 09:56:59 +0200 Subject: [PATCH 5/6] Revert "change internal function signature, evaluate zero_crossing for plot" This reverts commit 3ecd33e8d9db5abebe98b1a2599d2af12cf393b8. --- pyerrors/input/misc.py | 8 +++----- pyerrors/input/openQCD.py | 8 ++++---- 2 files changed, 7 insertions(+), 9 deletions(-) diff --git a/pyerrors/input/misc.py b/pyerrors/input/misc.py index 08e497da..6fd9ab74 100644 --- a/pyerrors/input/misc.py +++ b/pyerrors/input/misc.py @@ -29,7 +29,7 @@ def _reorder_Ysl(Ysl, tmax, nn): return Ysl_samples -def plot_Ysl(Ysl, expE_dict, nn, dn, eps, tmax, xmin, r_start, r_stop, r_step, names, spatial_extent, fit_range, c): +def plot_Ysl(Ysl, expE_dict, nn, dn, eps, tmax, xmin, r_start, r_stop, r_step, names, spatial_extent): """ Plot the plateaus needed for the fit of t0. """ @@ -66,9 +66,7 @@ def plot_Ysl(Ysl, expE_dict, nn, dn, eps, tmax, xmin, r_start, r_stop, r_step, n exp_vals = [] ts = [] - zero_crossing = np.argmax(np.array([o.value - c for o in t2expE_arr]) > 0.0) - - for i in range(zero_crossing-fit_range, zero_crossing+fit_range): + for i in range(len(t2expE_arr)): t = dn*eps*i t2expE_arr[i].gm() exp_vals.append(t2expE_arr[i]) @@ -78,7 +76,7 @@ def plot_Ysl(Ysl, expE_dict, nn, dn, eps, tmax, xmin, r_start, r_stop, r_step, n plt.fill_between(prange, t2expE_arr[i].value - t2expE_arr[i].dvalue, t2expE_arr[i].value + t2expE_arr[i].dvalue, alpha = 0.3, color = f"C{i:02d}") plt.hlines(t2expE_arr[i], prange[0], prange[1], linestyle = "dashed", colors=f"C{i:02d}", label = r"$t^2\langle E(t)\rangle$") - plt.ylim([t2expE_arr[zero_crossing-fit_range].value-.03,t2expE_arr[zero_crossing+fit_range].value+.01]) + plt.ylim([t2expE_arr[0].value-.03,t2expE_arr[-1].value+.01]) plt.ylabel(r"$t^2E(t)$") plt.xlabel("$x_{0}/a$") plt.xticks([i * tmax/4 for i in range(5)]) diff --git a/pyerrors/input/openQCD.py b/pyerrors/input/openQCD.py index a5504094..8a191340 100644 --- a/pyerrors/input/openQCD.py +++ b/pyerrors/input/openQCD.py @@ -232,7 +232,7 @@ def read_rwms(path, prefix, version='2.0', names=None, **kwargs): return result -def _extract_flowed_energy_density(path, prefix, dtr_read, xmin, spatial_extent, fit_range, c, postfix='ms', **kwargs): +def _extract_flowed_energy_density(path, prefix, dtr_read, xmin, spatial_extent, postfix='ms', **kwargs): """Extract a dictionary with the flowed Yang-Mills action density from given .ms.dat files. Returns a dictionary with Obs as values and flow times as keys. @@ -431,7 +431,7 @@ def _extract_flowed_energy_density(path, prefix, dtr_read, xmin, spatial_extent, expE_dict[n * dn * eps] = new_obs / (spatial_extent ** 3) if kwargs.get("plot_Ysl", False): - plot_Ysl(Ysl, expE_dict, nn, dn, eps, tmax, xmin, r_start_index, r_stop_index, r_step, rep_names, spatial_extent, fit_range, c) + plot_Ysl(Ysl, expE_dict, nn, dn, eps, tmax, xmin, r_start_index, r_stop_index, r_step, rep_names, spatial_extent) return expE_dict @@ -502,7 +502,7 @@ def extract_t0(path, prefix, dtr_read, xmin, spatial_extent, fit_range=5, postfi Extracted t0 """ - E_dict = _extract_flowed_energy_density(path, prefix, dtr_read, xmin, spatial_extent, fit_range, c, postfix, **kwargs) + E_dict = _extract_flowed_energy_density(path, prefix, dtr_read, xmin, spatial_extent, postfix, **kwargs) t2E_dict = {} for t in sorted(E_dict.keys()): t2E_dict[t] = t ** 2 * E_dict[t] - c @@ -575,7 +575,7 @@ def extract_w0(path, prefix, dtr_read, xmin, spatial_extent, fit_range=5, postfi Extracted w0 """ - E_dict = _extract_flowed_energy_density(path, prefix, dtr_read, xmin, spatial_extent, fit_range, c, postfix, **kwargs) + E_dict = _extract_flowed_energy_density(path, prefix, dtr_read, xmin, spatial_extent, postfix, **kwargs) ftimes = sorted(E_dict.keys()) From 36d8b301709b81e83e1ff399e4fae35c38d08224 Mon Sep 17 00:00:00 2001 From: Justus Kuhlmann Date: Thu, 20 Aug 2026 11:28:49 +0200 Subject: [PATCH 6/6] Only show 7 equidistant flow times in plot, refactor with helper functions --- pyerrors/input/misc.py | 89 ++++++++++++++++++++++++++++++------------ 1 file changed, 63 insertions(+), 26 deletions(-) diff --git a/pyerrors/input/misc.py b/pyerrors/input/misc.py index 6fd9ab74..79e9f38f 100644 --- a/pyerrors/input/misc.py +++ b/pyerrors/input/misc.py @@ -13,40 +13,61 @@ from ..obs import Obs -def _reorder_Ysl(Ysl, tmax, nn): - Ysl_samples = [] - for rep in range(len(Ysl)): - Ysl_array = [] - for cnfg in range(len(Ysl[rep])): - yarr = [] - for x0 in range(tmax): - yarr.append([]) - for flow in range(nn + 1): - ind = flow*tmax+x0 - yarr[x0].append(Ysl[rep][cnfg][ind]) - Ysl_array.append(yarr) - Ysl_samples.append(Ysl_array) - return Ysl_samples - - def plot_Ysl(Ysl, expE_dict, nn, dn, eps, tmax, xmin, r_start, r_stop, r_step, names, spatial_extent): """ - Plot the plateaus needed for the fit of t0. - """ - - prange = [xmin, tmax-xmin] + Plot the plateaux of 7 equidistant times in the flow time. Can be used to evaluate the quality of the plateau in x0 that is used + in a subsequent fit of t^2 to find t0. - Ysl_samples = _reorder_Ysl(Ysl, tmax, nn) - ts = list(expE_dict.keys()) + Parameters + ---------- + Ysl: list[list[list[float]]] + Array of read values, which are ordered as Ysl[rep][cnfg][flow*tmax+x0]. + expE_dict: dict + Already transformed array of action densities. + nn: int + Number of measurements. + dn: int + Steps of the integrator in flow time that are taken between measurements. + eps: float + Integrator stepsize in flow time. + tmax: int + Time extent of the inspected lattice. + xmin: int + Start of the (symmetrically chosen) plateau in terms ofthe euclidean time. + r_start: list[int] + Minimum index of the configuration to consider per replica. + r_stop: list[int] + Maximum index of the configuration to consider per replica. + r_step: int + Stepsize of the configurations in each replica. + names: list[str] + Replica names for the construction of observables. + spatial_extent: int + Spatial extent L/a of the lattice in units of the lattice spacing. + """ - t2E_arr = [] - t2expE_arr = [] - for n, t in zip(range(nn + 1), ts, strict=True): + def _disentangle_Ysl_inds(Ysl, tmax, ind_list): + Ysl_samples = [] + for rep in range(len(Ysl)): + Ysl_array = [] + for cnfg in range(len(Ysl[rep])): + yarr = [] + for x0 in range(tmax): + yarr.append([]) + for flow in ind_list: + ind = flow*tmax+x0 + yarr[x0].append(Ysl[rep][cnfg][ind]) + Ysl_array.append(yarr) + Ysl_samples.append(Ysl_array) + return Ysl_samples + + + def _Ecorr(Ysl_samples, tmax, n, r_start, r_stop, r_step, names, spatial_extent): corr_obs = [] for x0 in range(tmax): samples = [] rep_idls = [] - for rep in range(len(Ysl)): + for rep in range(len(Ysl_samples)): samples.append([]) rep_idls.append([]) for cnfg in range(len(Ysl_samples[rep])): @@ -56,9 +77,25 @@ def plot_Ysl(Ysl, expE_dict, nn, dn, eps, tmax, xmin, r_start, r_stop, r_step, n o = Obs(samples, names, rep_idls) corr_obs.append(o) E = Corr(corr_obs)/(spatial_extent**3) + return E + + num_t_shown = 7 + delta_n = int(nn/num_t_shown) # index stepsize in flow time + ind_list = [(i+1)*delta_n for i in range(num_t_shown)] + ts = [list(expE_dict.keys())[n] for n in ind_list] + prange = [xmin, tmax-xmin] + + Ysl_samples = _disentangle_Ysl_inds(Ysl, tmax, ind_list) + + t2E_arr = [] + for n, t in zip(range(len(ind_list)), ts, strict=True): + E=_Ecorr(Ysl_samples, tmax, n, r_start, r_stop, r_step, names, spatial_extent) t2E = t**2*E t2E.gm() t2E_arr.append(t2E) + + t2expE_arr = [] + for t in ts: t2expE = t**2*expE_dict[t] t2expE.gm() t2expE_arr.append(t2expE)