Update plots

anvil/config.py

@@ -185,15 +185,36 @@ def parse_bootstrap_sample_size(self, n_boot: int):
         log.warning(f"Using user specified bootstrap sample size: {n_boot}")
         return n_boot
 
-    def produce_bootstrap_seed(
-        self, manual_bootstrap_seed: (int, type(None)) = None):
+    def produce_bootstrap_seed(self, manual_bootstrap_seed: (int, type(None)) = None):
         if manual_bootstrap_seed is None:
             return randint(0, maxsize)
         # numpy is actually this strict but let's keep it sensible.
-        if (manual_bootstrap_seed < 0) or (manual_bootstrap_seed > 2**32):
+        if (manual_bootstrap_seed < 0) or (manual_bootstrap_seed > 2 ** 32):
             raise ConfigError("Seed is outside of appropriate range: [0, 2 ** 32]")
         return manual_bootstrap_seed
 
+    def parse_cosh_fit_min_separation(self, n: int, training_geometry):
+        """The smallest lattice separation to include in when fitting a cosh function
+        to the correlator, so as to the extract the correlation length.
+
+        See also: ``produce_cosh_fit_window``.
+        """
+        if n > training_geometry.length // 2 - 2:
+            raise ConfigError("Not enough points to for a three-parameter fit.")
+        return n
+
+    def produce_cosh_fit_window(self, training_geometry, cosh_fit_min_separation=None):
+        """Window of values corresponding to lattice separations, within which to fit
+        a cosh function to the correlator, so as to the extract the correlation length.
+        """
+        if cosh_fit_min_separation is None:
+            return slice(1, None)  # include all but (0, 0) separation by default
+
+        return slice(
+            cosh_fit_min_separation,
+            training_geometry.length - cosh_fit_min_separation + 1,
+        )
+
     @element_of("windows")
     def parse_window(self, window: float):
         """A numerical factor featuring in the calculation of the optimal 'window'
@@ -215,4 +236,4 @@ def produce_use_multiprocessing(self):
         """Don't use mp on MacOS"""
         if platform.system() == "Darwin":
             return False
-        return True
+        return True

anvil/observables.py

@@ -14,38 +14,46 @@
 
 log = logging.getLogger(__name__)
 
+
 def cosh_shift(x, xi, A, c):
     return A * np.cosh(-x / xi) + c
 
 
-def fit_zero_momentum_correlator(zero_momentum_correlator, training_geometry):
-    # TODO should I bootstrap this whole process...?
-
-    T = training_geometry.length
-    # TODO: would be good to specify this in runcard
-    t0 = T // 4
-    window = slice(t0, T - t0 + 1)
-
-    t = np.arange(T)
-    y = zero_momentum_correlator.mean(axis=-1)
-    yerr = zero_momentum_correlator.std(axis=-1)
-
-    try:
-        popt, pcov = optim.curve_fit(
-            cosh_shift,
-            xdata=t[window] - T // 2,
-            ydata=y[window],
-            sigma=yerr[window],
+def fit_zero_momentum_correlator(
+    zero_momentum_correlator, training_geometry, cosh_fit_window=slice(1, None)
+):
+    t = np.arange(training_geometry.length) - training_geometry.length // 2
+
+    # fit for each correlation func in the bootstrap ensemble
+    optimised_parameters = []
+    for correlator in zero_momentum_correlator.transpose():
+        try:
+            popt, pcov = optim.curve_fit(
+                cosh_shift,
+                xdata=t[cosh_fit_window],
+                ydata=correlator[cosh_fit_window],
+            )
+            optimised_parameters.append(popt)
+        except RuntimeError:
+            pass
+
+    n_boot = zero_momentum_correlator.shape[-1]
+    failures = n_boot - len(optimised_parameters)
+    if failures > 0:
+        log.warning(
+            f"Failed to fit cosh to correlation function for {failures}/{n_boot} members of the bootstrap ensemble."
         )
-        return (popt, pcov, t0)
-    except RuntimeError:
-        log.warning("Failed to fit cosh to correlation function.")
-        return None
+        if failures >= n_boot - 1:
+            log.warning("Too many failures: no fit parameters will be returned.")
+            return None
+
+    xi, A, c = np.array(optimised_parameters).transpose()
+    return xi, A, c
 
 
 def correlation_length_from_fit(fit_zero_momentum_correlator):
-    popt, pcov, _ = fit_zero_momentum_correlator
-    return popt[0], np.sqrt(pcov[0, 0])
+    xi, _, _ = fit_zero_momentum_correlator
+    return xi
 
 
 def autocorrelation(chain):
@@ -111,7 +119,7 @@ def magnetic_susceptibility(magnetization, abs_magnetization_squared):
 
 
 def magnetization_series(configs):
-    return configs.sum(axis=1)
+    return configs.sum(axis=1).numpy()
 
 
 def magnetization_autocorr(magnetization_series):
@@ -144,7 +152,9 @@ def __two_point_correlator(
             axis=-1  # sample average
         )
 
-    return correlator.reshape((training_geometry.length, training_geometry.length, -1))
+    return correlator.reshape(
+        (training_geometry.length, training_geometry.length, -1)
+    ).numpy()
 
 
 def two_point_correlator(
@@ -210,16 +220,6 @@ def ising_energy(two_point_correlator):
     return (two_point_correlator[1, 0] + two_point_correlator[0, 1]) / 2
 
 
-def inverse_pole_mass(effective_pole_mass, training_geometry):
-    T = training_geometry.length
-    t0 = T // 4
-    window = slice(t0, T - t0 + 1)
-
-    xi = np.reciprocal(effective_pole_mass)[window]
-
-    return np.nanmean(xi, axis=0)  # average over "large" t points
-
-
 def second_moment_correlation_length(two_point_correlator, susceptibility):
     """Second moment correlation length, defined as the normalised second
     moment of the two point correlator."""