Feat/plateaux for t0 extraction

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).

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]