Merge pull request #76 from wilsonmr/plots-update

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

@@ -5,47 +5,97 @@
 """
 import numpy as np
 from scipy.signal import correlate
-from math import ceil, pi, sin
+import scipy.optimize as optim
 import logging
 
 from anvil.utils import bootstrap_sample, Multiprocessing
 
-import scipy.optimize as optim
 
 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)
+def fit_zero_momentum_correlator(
+    zero_momentum_correlator, training_geometry, cosh_fit_window
+):
+    r"""Uses scipy.optimize.curve_fit to fit a cosh function (i.e. exponential decay
+    with periodicity) to each correlator in the bootrap ensemble.
+
+    The correlator decays as a pure exponential in the limit of large separations,
+    and the characteristic scale of this decay is the correlation length, whose
+    reciprocal is a.k.a the (effective) pole mass.
+
+    Parameters
+    ----------
+    zero_momentum_correlator
+        The two point correlation function at zero spatial momentum, i.e. the
+        correlation between 1-d 'slices'.
+    training_geometry
+        The anvil.geometry object defining the lattice.
+    cosh_fit_window: slice object
+        A slice object which selects the points (i.e. separations) to include in the
+        fit. In general the signal at short separations will be contaminated by
+        shorter modes and should not be included in the fit.
+
+    Returns
+    -------
+    xi: list
+        List of optimal correlation lengths for each member of the bootstrap ensemble
+        for whom the fitting process converged successfully.
+    A: list
+        Same as above, but for the amplitude of the cosh function.
+    c: list
+        Same as above, but for the global shift in the fit (which should correspond
+        to the absolute value of the magnetization, squared.
+
+    See also
+    --------
+    :py:func:`anvil.observables.cosh_shift` : the function being fit to the data.
+    """
+    t = np.arange(training_geometry.length) - training_geometry.length // 2
+
+    # fit for each correlation func in the bootstrap ensemble
+    xi, A, c = [], [], []
+    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],
+            )
+            xi.append(popt[0])
+            A.append(popt[1])
+            c.append(popt[2])
+        except RuntimeError:
+            pass
+
+    n_boot = zero_momentum_correlator.shape[-1]
+    n_fits = len(xi)
+    log.info(
+        f"Cosh fit succeeded for {n_fits}/{n_boot} members of the bootstrap ensemble."
+    )
+    return xi, A, c
 
-    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],
-        )
-        return (popt, pcov, t0)
-    except RuntimeError:
-        log.warning("Failed to fit cosh to correlation function.")
-        return None
+def correlation_length_from_fit(fit_zero_momentum_correlator):
+    """Returns numpy array containing a value for the  correlation length for each member
+    of the bootstrap ensemble for whom :py:func:`fit_zero_momentum_correlator` successfully
+    converged.
+    """
+    xi, _, _ = fit_zero_momentum_correlator
+    return np.array(xi)
 
 
-def correlation_length_from_fit(fit_zero_momentum_correlator):
-    popt, pcov, _ = fit_zero_momentum_correlator
-    return popt[0], np.sqrt(pcov[0, 0])
+def abs_magnetization_sq_from_fit(fit_zero_momentum_correlator):
+    """Returns numpy array containing a value for the absolute magnetization squared
+    for each member of the bootstrap ensemble for whom :py:func:`fit_zero_momentum_correlator`
+    successfully converged.
+    """
+    _, _, c = fit_zero_momentum_correlator
+    return np.array(c)
 
 
 def autocorrelation(chain):
@@ -102,16 +152,16 @@ def magnetization(configs, bootstrap_sample_size, bootstrap_seed):
     )
 
 
-def abs_magnetization_squared(magnetization):
+def abs_magnetization_sq(magnetization):
     return np.abs(magnetization).mean(axis=-1) ** 2  # <|m|>^2
 
 
-def magnetic_susceptibility(magnetization, abs_magnetization_squared):
-    return (magnetization ** 2).mean(axis=-1) - abs_magnetization_squared
+def magnetic_susceptibility(magnetization, abs_magnetization_sq):
+    return (magnetization ** 2).mean(axis=-1) - abs_magnetization_sq
 
 
 def magnetization_series(configs):
-    return configs.sum(axis=1)
+    return configs.sum(axis=1).numpy()
 
 
 def magnetization_autocorr(magnetization_series):
@@ -144,7 +194,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(
@@ -173,8 +225,8 @@ def two_point_correlator(
     return correlator.reshape((training_geometry.length, training_geometry.length, -1))
 
 
-def two_point_connected_correlator(two_point_correlator, abs_magnetization_squared):
-    return two_point_correlator - abs_magnetization_squared.view(1, 1, -1)
+def two_point_connected_correlator(two_point_correlator, abs_magnetization_sq):
+    return two_point_correlator - abs_magnetization_sq.view(1, 1, -1)
 
 
 def zero_momentum_correlator(two_point_correlator):
@@ -210,16 +262,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."""
@@ -239,11 +281,11 @@ def second_moment_correlation_length(two_point_correlator, susceptibility):
 def low_momentum_correlation_length(two_point_correlator, susceptibility):
     """A low-momentum estimate for the correlation length."""
     L = two_point_correlator.shape[0]
-    kernel = np.cos(2 * pi / L * np.arange(L)).reshape(L, 1, 1)
+    kernel = np.cos(2 * np.pi / L * np.arange(L)).reshape(L, 1, 1)
 
     g_tilde_00 = susceptibility
     g_tilde_10 = (kernel * two_point_correlator).sum(axis=(0, 1))
 
-    xi_sq = (g_tilde_00 / g_tilde_10 - 1) / (4 * sin(pi / L) ** 2)
+    xi_sq = (g_tilde_00 / g_tilde_10 - 1) / (4 * np.sin(np.pi / L) ** 2)
 
     return np.sqrt(xi_sq)