diff --git a/pyerrors/input/misc.py b/pyerrors/input/misc.py index c7c7f5b2..5995c81d 100644 --- a/pyerrors/input/misc.py +++ b/pyerrors/input/misc.py @@ -8,10 +8,110 @@ 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] + + 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.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)]) + 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]