From a7d1e46fecb42cbb514ba275a3bd707904ccb19e Mon Sep 17 00:00:00 2001
From: marshrossney <17361029+marshrossney@users.noreply.github.com>
Date: Mon, 17 May 2021 20:07:04 +0100
Subject: [PATCH 01/17] update magnetization series with rows

---
 anvil/observables.py |  7 +++++--
 anvil/plot.py        | 32 ++++++++++++++++++++------------
 2 files changed, 25 insertions(+), 14 deletions(-)

diff --git a/anvil/observables.py b/anvil/observables.py
index 02930ba..9be30b6 100644
--- a/anvil/observables.py
+++ b/anvil/observables.py
@@ -14,6 +14,7 @@
 
 log = logging.getLogger(__name__)
 
+
 def cosh_shift(x, xi, A, c):
     return A * np.cosh(-x / xi) + c
 
@@ -111,7 +112,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 +145,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(
diff --git a/anvil/plot.py b/anvil/plot.py
index bea99ab..4591a55 100644
--- a/anvil/plot.py
+++ b/anvil/plot.py
@@ -218,23 +218,31 @@ def plot_magnetization(magnetization_series):
     ax.set_title("Magnetization")
     ax.set_ylabel("Frequency")
     ax.set_xlabel("$M(t)$")
-    ax.hist(magnetization_series.numpy(), histtype="stepfilled", edgecolor="black")
+    ax.hist(magnetization_series, histtype="stepfilled", edgecolor="black")
     return fig
 
 
 @figure
 def plot_magnetization_series(magnetization_series, sample_interval):
-    chain_indices = np.arange(magnetization_series.shape[-1]) * sample_interval
-    fig, ax = plt.subplots()
-    ax.set_title("Magnetization")
-    ax.set_ylabel("$M(t)$")
-    ax.set_xlabel("$t$")
-    ax.plot(
-        chain_indices,
-        magnetization_series,
-        linestyle="-",
-        marker="",
-    )
+    n_rows = 5
+    if magnetization_series.size % n_rows != 0:
+        magnetization_series = np.pad(
+            magnetization_series,
+            (0, magnetization_series.size % n_rows),
+            "empty",
+        )
+    magnetization_series = magnetization_series.reshape(n_rows, -1)
+    t = (np.arange(magnetization_series.size) * sample_interval).reshape(n_rows, -1)
+
+    fig, axes = plt.subplots(n_rows, 1, sharey=True)
+    for ax, x, y in zip(axes, t, magnetization_series):
+        ax.plot(x, y, linestyle="-", marker="")
+        ax.margins(0, 0)
+
+    axes[0].set_title("Magnetization")
+    axes[n_rows // 2].set_ylabel("$M(t)$")
+    axes[-1].set_xlabel("$t$")
+    fig.tight_layout()
     return fig
 
 

From f1ef44708ccbc49e29874214b0944fd2bfd94f65 Mon Sep 17 00:00:00 2001
From: marshrossney <17361029+marshrossney@users.noreply.github.com>
Date: Mon, 17 May 2021 22:20:17 +0100
Subject: [PATCH 02/17] split correlator and error into two figures

---
 anvil/plot.py | 75 +++++++++++++++++++++++++++++++--------------------
 1 file changed, 46 insertions(+), 29 deletions(-)

diff --git a/anvil/plot.py b/anvil/plot.py
index 4591a55..e7225d3 100644
--- a/anvil/plot.py
+++ b/anvil/plot.py
@@ -9,7 +9,7 @@
 import torch
 import numpy as np
 import matplotlib.pyplot as plt
-from mpl_toolkits.axes_grid1 import make_axes_locatable
+import matplotlib as mpl
 from mpl_toolkits.axes_grid1.anchored_artists import AnchoredSizeBar
 from matplotlib.font_manager import FontProperties
 
@@ -167,47 +167,64 @@ def plot_effective_pole_mass(training_geometry, effective_pole_mass):
 
 @figure
 def plot_two_point_correlator(two_point_correlator):
-    """Represent the two point function and it's error in x and t as heatmaps
-    of the respective matrices. Returns a figure with two plots, the left plot
-    is the mean two point function across bootstrap and the right plot is the
-    standard deviation divide by the mean (fractional error)
-
+    """Represent the two point correlator as a heatmap. The data shown is the mean
+    of the bootstrap sample of correlation functions, and is normalised so that
+    G(0, 0) = 1. The colour axis is scaled using a symmetric log scale, with a linear
+    region spanning [-0.01, 0.01].
     """
     corr = two_point_correlator.mean(axis=-1)
-    error = two_point_correlator.std(axis=-1)
-    norm = corr[0, 0]
+    corr /= corr[0, 0]
 
     L = corr.shape[0]
     corr = np.roll(corr, (-L // 2 - 1, -L // 2 - 1), (0, 1))
-    error = np.roll(error, (-L // 2 - 1, -L // 2 - 1), (0, 1))
 
-    fig, (ax1, ax2) = plt.subplots(1, 2, sharex=True, sharey=True)
-    ax1.set_title("$G(x, t)$")
-    ax2.set_title(r"$| \sigma_G / G |$")
-    ax1.set_xlabel("$x_1$")
-    ax2.set_xlabel("$x_1$")
-    ax1.set_ylabel("$x_2$")
+    fig, ax = plt.subplots()
+    ax.set_title("$G(x)$")
+    ax.set_xlabel("$x_1$")
+    ax.set_ylabel("$x_2$")
 
     tick_positions = [0, L // 2 - 1, L - 1]
     tick_labels = [r"$-\frac{L + 1}{2}$", 0, r"$\frac{L}{2}$"]
-    ax1.set_xticks(tick_positions)
-    ax1.set_yticks(tick_positions)
-    ax1.set_xticklabels(tick_labels)
-    ax1.set_yticklabels(tick_labels)
-    ax2.tick_params(axis="y", width=0)
+    ax.set_xticks(tick_positions)
+    ax.set_yticks(tick_positions)
+    ax.set_xticklabels(tick_labels)
+    ax.set_yticklabels(tick_labels)
 
-    im1 = ax1.imshow(corr / norm)
-    im2 = ax2.imshow(np.abs(error / corr))
+    norm = mpl.colors.SymLogNorm(linthresh=0.01, base=10)
+    im = ax.imshow(corr, norm=norm)
+    fig.colorbar(im, ax=ax, pad=0.01)
 
-    div1 = make_axes_locatable(ax1)
-    cax1 = div1.append_axes("right", size="5%", pad=0.05)
-    fig.colorbar(im1, cax=cax1)
+    return fig
 
-    div2 = make_axes_locatable(ax2)
-    cax2 = div2.append_axes("right", size="7%", pad=0.05)
-    fig.colorbar(im2, cax=cax2)
 
-    fig.tight_layout()
+@figure
+def plot_two_point_correlator_error(two_point_correlator):
+    """Heatmap of the error in the two point correlator for each separation (x_1, x_2).
+    The error is computed as the standard deviation over the bootstrap sample. The data
+    shown is this error divided by the mean of the bootstrap sample, i.e. the
+    fractional error.
+    """
+    corr = two_point_correlator.mean(axis=-1)
+    error = two_point_correlator.std(axis=-1)
+
+    L = corr.shape[0]
+    corr = np.roll(corr, (-L // 2 - 1, -L // 2 - 1), (0, 1))
+    error = np.roll(error, (-L // 2 - 1, -L // 2 - 1), (0, 1))
+
+    fig, ax = plt.subplots()
+    ax.set_title(r"$| \sigma_G(x) / G(x) |$")
+    ax.set_xlabel("$x_1$")
+    ax.set_ylabel("$x_2$")
+
+    tick_positions = [0, L // 2 - 1, L - 1]
+    tick_labels = [r"$-\frac{L + 1}{2}$", 0, r"$\frac{L}{2}$"]
+    ax.set_xticks(tick_positions)
+    ax.set_yticks(tick_positions)
+    ax.set_xticklabels(tick_labels)
+    ax.set_yticklabels(tick_labels)
+
+    im = ax.imshow(np.abs(error / corr))
+    fig.colorbar(im, ax=ax, pad=0.01)
 
     return fig
 

From 738cad131f5ce144aed26eb25ee9a554f9db7bf0 Mon Sep 17 00:00:00 2001
From: marshrossney <17361029+marshrossney@users.noreply.github.com>
Date: Mon, 17 May 2021 23:23:47 +0100
Subject: [PATCH 03/17] plot with multiple estimates of correlation length

---
 anvil/observables.py | 12 +--------
 anvil/plot.py        | 59 ++++++++++++++++++++++++++++++++++++++++++++
 anvil/table.py       | 14 ++++++++---
 3 files changed, 70 insertions(+), 15 deletions(-)

diff --git a/anvil/observables.py b/anvil/observables.py
index 9be30b6..548f6e2 100644
--- a/anvil/observables.py
+++ b/anvil/observables.py
@@ -20,7 +20,7 @@ def cosh_shift(x, xi, A, c):
 
 
 def fit_zero_momentum_correlator(zero_momentum_correlator, training_geometry):
-    # TODO should I bootstrap this whole process...?
+    # Bootstrap this whole process
 
     T = training_geometry.length
     # TODO: would be good to specify this in runcard
@@ -213,16 +213,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."""
diff --git a/anvil/plot.py b/anvil/plot.py
index e7225d3..d8f1c0d 100644
--- a/anvil/plot.py
+++ b/anvil/plot.py
@@ -165,6 +165,65 @@ def plot_effective_pole_mass(training_geometry, effective_pole_mass):
     return fig
 
 
+@figure
+def plot_correlation_length(
+    effective_pole_mass,
+    low_momentum_correlation_length,
+    correlation_length_from_fit,
+):
+    """Plots three estimates of correlation length on the same figure.
+
+    These are:
+        1. Estimate from fitting a cosh function to the correlation between
+            1-dimensional slices, using py:func:`correlation_length_from_fit`
+        2. Reciprocal of the effective pole mass estimator, using
+            :py:func:`effective_pole_mass`
+        3. Low momentum estimate, using :py:func:`low_momentum_correlation_length`
+
+    2. is evaluated at a specific value of the separation, x_2.
+    """
+    xi_arcosh = np.reciprocal(effective_pole_mass)
+
+    fig, ax = plt.subplots()
+
+    arcosh_points = ax.errorbar(
+        x=range(1, xi_arcosh.shape[0] + 1),
+        y=xi_arcosh.mean(axis=-1),
+        yerr=xi_arcosh.std(axis=-1),
+        zorder=3,
+    )
+
+    xi_lm = low_momentum_correlation_length.mean()
+    e_lm = low_momentum_correlation_length.std()
+    lm_hline = ax.axhline(xi_lm, linestyle="-", marker="", color="grey", zorder=2)
+    lm_fill = ax.fill_between(
+        ax.get_xlim(),
+        xi_lm + e_lm,
+        xi_lm - e_lm,
+        color="grey",
+        alpha=0.3,
+        zorder=1,
+    )
+
+    xi_fit = correlation_length_from_fit[0]  # .mean(  # TODO: update when bootstrapped)
+    e_fit = correlation_length_from_fit[1]  # .std()
+    fit_hline = ax.axhline(xi_fit, linestyle="-", marker="", color="orange", zorder=2)
+    fit_fill = ax.fill_between(
+        ax.get_xlim(),
+        xi_fit + e_fit,
+        xi_fit - e_fit,
+        color="orange",
+        alpha=0.3,
+        zorder=1,
+    )
+
+    ax.legend(
+        handles=[arcosh_points, (lm_fill, lm_hline), (fit_fill, fit_hline)],
+        labels=["Estimate using arcosh", "Low momentum estimate", "Estimate from fit"],
+    )
+    return fig
+
+
 @figure
 def plot_two_point_correlator(two_point_correlator):
     """Represent the two point correlator as a heatmap. The data shown is the mean
diff --git a/anvil/table.py b/anvil/table.py
index 3346aa4..86f96b2 100644
--- a/anvil/table.py
+++ b/anvil/table.py
@@ -79,20 +79,26 @@ def table_magnetization(abs_magnetization_squared, magnetic_susceptibility):
 
 @table
 def table_correlation_length(
-    inverse_pole_mass,
+    effective_pole_mass,
     second_moment_correlation_length,
     low_momentum_correlation_length,
     correlation_length_from_fit,
     training_geometry,
 ):
-    """Tabulate four estimators of correlation length, with values and errors
+    """Tabulate four estimates of correlation length, with values and errors
     taken as the mean and standard deviation of the bootstrap sample.
     """
-    # TODO could add column with inverse (i.e. effective pole mass)
+    # Take the mean of the arcosh estimator over "large" separations
+    x0 = training_geometry.length // 4
+    window = slice(x0, training_geometry.length - x0 + 1)
+    xi_arcosh = np.nanmean(
+        np.reciprocal(effective_pole_mass)[window],
+        axis=0,
+    )
 
     res = [
         list(correlation_length_from_fit),
-        [inverse_pole_mass.mean(), inverse_pole_mass.std()],
+        [xi_arcosh.mean(), xi_arcosh.std()],
         [
             second_moment_correlation_length.mean(),
             second_moment_correlation_length.std(),

From dad6a8acdcc36fa5b9b57ac06e6224bb16accc4b Mon Sep 17 00:00:00 2001
From: marshrossney <17361029+marshrossney@users.noreply.github.com>
Date: Thu, 20 May 2021 14:13:45 +0100
Subject: [PATCH 04/17] bootstrap the cosh fit

---
 anvil/config.py      | 29 +++++++++++++++++++----
 anvil/observables.py | 55 +++++++++++++++++++++++++-------------------
 anvil/plot.py        | 37 ++++++++++++++++++-----------
 anvil/table.py       |  8 +++----
 4 files changed, 84 insertions(+), 45 deletions(-)

diff --git a/anvil/config.py b/anvil/config.py
index 7eeb68e..44e5d54 100644
--- a/anvil/config.py
+++ b/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
\ No newline at end of file
+        return True
diff --git a/anvil/observables.py b/anvil/observables.py
index 548f6e2..a7e580f 100644
--- a/anvil/observables.py
+++ b/anvil/observables.py
@@ -19,34 +19,41 @@ def cosh_shift(x, xi, A, c):
     return A * np.cosh(-x / xi) + c
 
 
-def fit_zero_momentum_correlator(zero_momentum_correlator, training_geometry):
-    # 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):
diff --git a/anvil/plot.py b/anvil/plot.py
index d8f1c0d..043f618 100644
--- a/anvil/plot.py
+++ b/anvil/plot.py
@@ -113,32 +113,43 @@ def plot_zero_momentum_correlator(
     zero_momentum_correlator,
     training_geometry,
     fit_zero_momentum_correlator,
+    cosh_fit_window,
 ):
     """Plot zero_momentum_2pf as a function of t. Points are means across bootstrap
     sample and errorbars are standard deviations across boostrap samples
     """
-    L = training_geometry.length
-
     fig, ax = plt.subplots()
 
     ax.errorbar(
-        x=np.arange(L),
+        x=np.arange(training_geometry.length),
         y=zero_momentum_correlator.mean(axis=-1),
         yerr=zero_momentum_correlator.std(axis=-1),
         linestyle="",
+        zorder=2,
         label="sample statistics",
     )
 
     if fit_zero_momentum_correlator is not None:
-        popt, pcov, x0 = fit_zero_momentum_correlator
-
-        x_2 = np.linspace(x0, L - x0, 100)
-        ax.plot(
-            x_2,
-            cosh_shift(x_2 - L // 2, *popt),
-            marker="",
-            label=r"fit: $A \cosh(-(x_2 - L/2) / \xi) + c$",
+        t = np.arange(training_geometry.length)[cosh_fit_window]
+        fit = []
+        for xi, A, c in zip(*fit_zero_momentum_correlator):
+            fit.append(
+                cosh_shift(t - training_geometry.length // 2, xi, A, c),
+            )
+        fit = np.array(fit).T  # (n_points, n_boot)
+
+        ax.fill_between(
+            t,
+            fit.mean(axis=1) - fit.std(axis=1),
+            fit.mean(axis=1) + fit.std(axis=1),
+            color="orange",
+            alpha=0.3,
+            zorder=1,
+            label=r"fit: $A \cosh(-(x_2 - L/2) / \xi) + c$"
+            + "\n"
+            + r"($1\sigma$ confidence)",
         )
+
     ax.set_yscale("log")
     ax.set_ylabel(r"$\sum_{x_1=0}^{L-1} G(x_1, x_2)$")
     ax.set_xlabel("$x_2$")
@@ -205,8 +216,8 @@ def plot_correlation_length(
         zorder=1,
     )
 
-    xi_fit = correlation_length_from_fit[0]  # .mean(  # TODO: update when bootstrapped)
-    e_fit = correlation_length_from_fit[1]  # .std()
+    xi_fit = correlation_length_from_fit.mean()
+    e_fit = correlation_length_from_fit.std()
     fit_hline = ax.axhline(xi_fit, linestyle="-", marker="", color="orange", zorder=2)
     fit_fill = ax.fill_between(
         ax.get_xlim(),
diff --git a/anvil/table.py b/anvil/table.py
index 86f96b2..683dec2 100644
--- a/anvil/table.py
+++ b/anvil/table.py
@@ -32,11 +32,11 @@ def table_autocorrelation(
 @table
 def table_fit(fit_zero_momentum_correlator, training_geometry):
     if fit_zero_momentum_correlator is not None:
-        popt, pcov, t0 = fit_zero_momentum_correlator
+        xi, A, c = fit_zero_momentum_correlator
 
         res = [
-            [popt[0], np.sqrt(pcov[0, 0])],
-            [popt[2], np.sqrt(pcov[2, 2])],
+            [xi.mean(), xi.std()],
+            [c.mean(), c.std()],
         ]
         df = pd.DataFrame(
             res,
@@ -97,7 +97,7 @@ def table_correlation_length(
     )
 
     res = [
-        list(correlation_length_from_fit),
+        [correlation_length_from_fit.mean(), correlation_length_from_fit.std()],
         [xi_arcosh.mean(), xi_arcosh.std()],
         [
             second_moment_correlation_length.mean(),

From 1eaac3e4c826fc14e39b605a50828bcb4e13fa9e Mon Sep 17 00:00:00 2001
From: marshrossney <17361029+marshrossney@users.noreply.github.com>
Date: Thu, 20 May 2021 14:21:36 +0100
Subject: [PATCH 05/17] standardise table row/col names

---
 anvil/table.py | 32 ++++++++++++++++----------------
 1 file changed, 16 insertions(+), 16 deletions(-)

diff --git a/anvil/table.py b/anvil/table.py
index 683dec2..7a9eda4 100644
--- a/anvil/table.py
+++ b/anvil/table.py
@@ -23,7 +23,7 @@ def table_autocorrelation(
 
     df = pd.DataFrame(
         [acceptance, tau_mag, tau_chain],
-        index=["acceptance", "tau_mag", "tau_chain"],
+        index=["acceptance", "tau_from_magnetization", "tau_from_chain"],
         columns=["value"],
     )
     return df
@@ -40,8 +40,8 @@ def table_fit(fit_zero_momentum_correlator, training_geometry):
         ]
         df = pd.DataFrame(
             res,
-            columns=["Mean", "Standard deviation"],
-            index=["xi_fit", "m_fit"],
+            columns=["mean", "error"],
+            index=["xi_from_fit", "magnetization_from_fit"],
         )
         return df
 
@@ -57,8 +57,8 @@ def table_two_point_scalars(ising_energy, susceptibility):
     ]
     df = pd.DataFrame(
         res,
-        columns=["Mean", "Standard deviation"],
-        index=["Ising energy", "susceptibility"],
+        columns=["mean", "error"],
+        index=["ising_energy", "susceptibility"],
     )
     return df
 
@@ -71,8 +71,8 @@ def table_magnetization(abs_magnetization_squared, magnetic_susceptibility):
     ]
     df = pd.DataFrame(
         res,
-        columns=["Mean", "Standard deviation"],
-        index=["<|m|>^2", "<m^2> - <|m|>^2"],
+        columns=["mean", "error"],
+        index=["abs_magnetization_squared", "magnetic_susceptibility"],
     )
     return df
 
@@ -108,15 +108,15 @@ def table_correlation_length(
 
     df = pd.DataFrame(
         res,
-        columns=["Mean", "Standard deviation"],
+        columns=["mean", "error"],
         index=[
-            "Estimate from fit",
-            "Estimate using arcosh",
-            "Second moment estimate",
-            "Low momentum estimate",
+            "xi_from_fit",
+            "xi_from_arcosh",
+            "xi_from_second_moment",
+            "xi_from_low_momentum",
         ],
     )
-    df["No. correlation lengths"] = training_geometry.length / df["Mean"]
+    df["n_correlation_lengths"] = training_geometry.length / df["mean"]
     return df
 
 
@@ -132,7 +132,7 @@ def table_zero_momentum_correlator(zero_momentum_correlator, training_geometry):
 
     df = pd.DataFrame(
         data,
-        columns=["Mean", "Standard deviation"],
+        columns=["mean", "error"],
         index=range(training_geometry.length),
     )
     return df
@@ -149,7 +149,7 @@ def table_effective_pole_mass(effective_pole_mass, training_geometry):
     data = np.concatenate((means, stds), axis=1)
     df = pd.DataFrame(
         data,
-        columns=["Mean", "Standard deviation"],
+        columns=["mean", "error"],
         index=range(1, training_geometry.length - 1),
     )
     return df
@@ -169,5 +169,5 @@ def table_two_point_correlator(training_geometry, two_point_correlator):
         for j in range(training_geometry.length):
             corr.append([float(means[i, j]), float(stds[i, j])])
             index.append((i, j))
-    df = pd.DataFrame(corr, columns=["Mean", "Standard deviation"], index=index)
+    df = pd.DataFrame(corr, columns=["mean", "error"], index=index)
     return df

From 8f582c6076d090085762c68303be77d8418d5d95 Mon Sep 17 00:00:00 2001
From: marshrossney <17361029+marshrossney@users.noreply.github.com>
Date: Thu, 20 May 2021 14:27:26 +0100
Subject: [PATCH 06/17] update report.md to include new plots

---
 examples/runcards/report.md | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/examples/runcards/report.md b/examples/runcards/report.md
index 9d71141..06f2d26 100644
--- a/examples/runcards/report.md
+++ b/examples/runcards/report.md
@@ -5,15 +5,18 @@
 {@table_two_point_scalars@}
 ## Two Point Correlator
 {@plot_two_point_correlator@}
+{@plot_two_point_correlator_error@}
 ## Zero momentum two point correlator
 {@table_fit@}
 {@plot_zero_momentum_correlator@}
 ## Correlation length
 {@plot_effective_pole_mass@}
 {@table_correlation_length@}
+{@plot_correlation_length@}
 ## Magnetisation
 {@table_magnetization@}
 {@plot_magnetization@}
 {@plot_magnetization_series@}
+## Autocorrelation
 {@plot_magnetization_autocorr@}
 {@plot_magnetization_integrated_autocorr@}

From 9c4ce28724f729bc897faeaa1e379bd4b479351a Mon Sep 17 00:00:00 2001
From: marshrossney <17361029+marshrossney@users.noreply.github.com>
Date: Thu, 20 May 2021 14:32:52 +0100
Subject: [PATCH 07/17] minor adjustment to autocorrelation plot

---
 anvil/plot.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/anvil/plot.py b/anvil/plot.py
index 043f618..36714de 100644
--- a/anvil/plot.py
+++ b/anvil/plot.py
@@ -395,11 +395,11 @@ def plot_magnetization_integrated_autocorr(
         ax.transData,
         sample_interval,
         f"sample interval: {sample_interval}",
-        "upper left",
+        "center right",
         pad=0.6,
         frameon=False,
         sep=4,
-        label_top=True,
+        label_top=False,
         fontproperties=FontProperties(size="x-large"),
     )
     ax.add_artist(scalebar)

From 30ba442570e43b23cdb90e109b184ef6d2eaf794 Mon Sep 17 00:00:00 2001
From: marshrossney <17361029+marshrossney@users.noreply.github.com>
Date: Thu, 20 May 2021 15:20:03 +0100
Subject: [PATCH 08/17] add docstrings to plot actions, remove plots due to be
 superseded

---
 anvil/plot.py | 138 +++++++++++++++-----------------------------------
 1 file changed, 41 insertions(+), 97 deletions(-)

diff --git a/anvil/plot.py b/anvil/plot.py
index 36714de..1ca4ca5 100644
--- a/anvil/plot.py
+++ b/anvil/plot.py
@@ -19,95 +19,6 @@
 from anvil.observables import cosh_shift
 
 
-def field_component(i, x_base, phi_model, base_neg, model_neg):
-    fig, ax = plt.subplots()
-
-    ax.hist(x_base, bins=50, density=True, histtype="step", label="base")
-    ax.hist(phi_model, bins=50, density=True, histtype="step", label="model, full")
-    ax.hist(
-        base_neg, bins=50, density=True, histtype="step", label="model, $M_{base} < 0$"
-    )
-    ax.hist(
-        model_neg, bins=50, density=True, histtype="step", label="model, $M_{mod} < 0$"
-    )
-    ax.set_title(f"Coordinate {i}")
-    ax.legend()
-    fig.tight_layout()
-    return fig
-
-
-def field_components(loaded_model, base_dist, lattice_size):
-    """Plot the distributions of base coordinates 'x' and output coordinates 'phi' and,
-    if known, plot the pdf of the target distribution."""
-    sample_size = 10000
-
-    # Generate a large sample from the base distribution and pass it through the trained model
-    with torch.no_grad():
-        x_base, _ = base_dist(sample_size)
-        sign = x_base.sum(dim=1).sign()
-        neg = (sign < 0).nonzero().squeeze()
-        phi_model, model_log_density = loaded_model(x_base, 0, neg)
-
-    base_neg = phi_model[neg]
-
-    sign = phi_model.sum(dim=1).sign()
-    neg = (sign < 0).nonzero().squeeze()
-    model_neg = phi_model[neg]
-
-    # Convert to shape (n_coords, sample_size * lattice_size)
-    # NOTE: this is all pointless for the 1-component scalar
-    x_base = x_base.reshape(sample_size * lattice_size, -1).transpose(0, 1)
-    phi_model = phi_model.reshape(sample_size * lattice_size, -1).transpose(0, 1)
-
-    base_neg = base_neg.reshape(1, -1)
-    model_neg = model_neg.reshape(1, -1)
-
-    for i in range(x_base.shape[0]):
-        yield field_component(i, x_base[i], phi_model[i], base_neg[i], model_neg[i])
-
-
-_plot_field_components = collect("field_components", ("training_context",))
-
-
-def example_configs(loaded_model, base_dist, training_geometry):
-    sample_size = 10
-
-    # Generate a large sample from the base distribution and pass it through the trained model
-    with torch.no_grad():
-        x_base, _ = base_dist(sample_size)
-        sign = x_base.sum(dim=1).sign()
-        neg = (sign < 0).nonzero().squeeze()
-        phi_model, model_log_density = loaded_model(x_base, 0, neg)
-
-    L = int(np.sqrt(phi_model.shape[1]))
-
-    phi_true = np.zeros((4, L, L))
-    phi_true[:, training_geometry.checkerboard] = phi_model[:4, : L ** 2 // 2]
-    phi_true[:, ~training_geometry.checkerboard] = phi_model[:4, L ** 2 // 2 :]
-
-    fig, axes = plt.subplots(2, 2, sharex=True, sharey=True)
-    for i, ax in enumerate(axes.flatten()):
-        conf = ax.imshow(phi_true[i])
-        fig.colorbar(conf, ax=ax)
-
-    fig.suptitle("Example configurations")
-
-    return fig
-
-
-_plot_example_configs = collect("example_configs", ("training_context",))
-
-
-@figure
-def plot_example_configs(_plot_example_configs):
-    return _plot_example_configs[0]
-
-
-@figuregen
-def plot_field_components(_plot_field_components):
-    yield from _plot_field_components[0]
-
-
 @figure
 def plot_zero_momentum_correlator(
     zero_momentum_correlator,
@@ -115,8 +26,14 @@ def plot_zero_momentum_correlator(
     fit_zero_momentum_correlator,
     cosh_fit_window,
 ):
-    """Plot zero_momentum_2pf as a function of t. Points are means across bootstrap
-    sample and errorbars are standard deviations across boostrap samples
+    """Plots the correlation function for pairs of one-dimensional 'slices', otherwise
+    referred to as the two point correlator at zero spatial momentum, as a function of
+    time. Points and errorbars are means and standard deviations across a boostrap
+    sample.
+
+    Also plots the :math:`1\sigma` confidence interval for a pure-exponential (cosh)
+    fit performed for each member of the bootstrap sample in
+    :py:func:`fit_zero_momentum_correlator`.
     """
     fig, ax = plt.subplots()
 
@@ -160,9 +77,8 @@ def plot_zero_momentum_correlator(
 
 @figure
 def plot_effective_pole_mass(training_geometry, effective_pole_mass):
-    """Plot effective pole mass as a function of x_2. The points are means
-    across bootstrap samples and the errorbars are standard deviation across
-    bootstrap.
+    """Plots the (effective) pole mass as a function of 'time' separation. Points and
+    error bars are means and standard deviations across a bootstrap sample.
     """
     fig, ax = plt.subplots()
     ax.errorbar(
@@ -182,7 +98,7 @@ def plot_correlation_length(
     low_momentum_correlation_length,
     correlation_length_from_fit,
 ):
-    """Plots three estimates of correlation length on the same figure.
+    """Plots three estimates of correlation length.
 
     These are:
         1. Estimate from fitting a cosh function to the correlation between
@@ -192,6 +108,9 @@ def plot_correlation_length(
         3. Low momentum estimate, using :py:func:`low_momentum_correlation_length`
 
     2. is evaluated at a specific value of the separation, x_2.
+
+    Points and error bars are means and standard deviations taken across a bootstrap
+    sample.
     """
     xi_arcosh = np.reciprocal(effective_pole_mass)
 
@@ -237,8 +156,8 @@ def plot_correlation_length(
 
 @figure
 def plot_two_point_correlator(two_point_correlator):
-    """Represent the two point correlator as a heatmap. The data shown is the mean
-    of the bootstrap sample of correlation functions, and is normalised so that
+    """Represents the two point correlator as a heatmap. The data shown is the mean
+    of a bootstrap sample of correlation functions, and is normalised so that
     G(0, 0) = 1. The colour axis is scaled using a symmetric log scale, with a linear
     region spanning [-0.01, 0.01].
     """
@@ -301,6 +220,14 @@ def plot_two_point_correlator_error(two_point_correlator):
 
 @figure
 def plot_magnetization(magnetization_series):
+    """Plots a histogram of the magnetization of each configuration in the Markov
+    chain resulting from the Metropolis-Hastings sampling phase.
+
+    See also
+    --------
+    :py:func:`plot_magnetization_series` displays the same data as a sequence,
+    i.e. in order of their appearance in the Markov chain.
+    """
     fig, ax = plt.subplots()
     ax.set_title("Magnetization")
     ax.set_ylabel("Frequency")
@@ -311,6 +238,13 @@ def plot_magnetization(magnetization_series):
 
 @figure
 def plot_magnetization_series(magnetization_series, sample_interval):
+    """Plots the magnetization of each configuration in the Markov chain over the
+    course of the Metropolis-Hastings sampling phase.
+
+    See also
+    --------
+    :py:func:`plot_magnetization` displays the same data as a histogram.
+    """
     n_rows = 5
     if magnetization_series.size % n_rows != 0:
         magnetization_series = np.pad(
@@ -337,6 +271,11 @@ def plot_magnetization_series(magnetization_series, sample_interval):
 def plot_magnetization_autocorr(
     magnetization_autocorr, magnetization_optimal_window, sample_interval
 ):
+    """Plots the autocorrelation function for the magnetization of the sequence of
+    configurations generated in the Metropolis-Hastings sampling phase. The x-axis
+    corresponds to a number of steps separating pairs of configurations in the
+    sequence.
+    """
     cut = max(10, 2 * magnetization_optimal_window)
     chain_indices = np.arange(cut) * sample_interval
 
@@ -377,6 +316,11 @@ def plot_magnetization_integrated_autocorr(
     magnetization_optimal_window,
     sample_interval,
 ):
+    """Plots the integrated autocorrelation function for the magnetization of the
+    sequence of configurations generated in the Metropolis-Hastings sampling phase.
+    The x axis represents the size of the 'window' in which the summation is performed,
+    i.e. the point at which the autocorrelation function is truncated.
+    """
     cut = max(10, 2 * np.max(magnetization_optimal_window))
     chain_indices = np.arange(cut) * sample_interval
     tau = magnetization_integrated_autocorr[magnetization_optimal_window]

From c530e0b832c7eab78602122e291f27d9d34c3f9d Mon Sep 17 00:00:00 2001
From: marshrossney <17361029+marshrossney@users.noreply.github.com>
Date: Thu, 20 May 2021 15:27:20 +0100
Subject: [PATCH 09/17] raw strings to get rid of pytest warnings

---
 anvil/plot.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/anvil/plot.py b/anvil/plot.py
index 1ca4ca5..fa7f1c5 100644
--- a/anvil/plot.py
+++ b/anvil/plot.py
@@ -26,7 +26,7 @@ def plot_zero_momentum_correlator(
     fit_zero_momentum_correlator,
     cosh_fit_window,
 ):
-    """Plots the correlation function for pairs of one-dimensional 'slices', otherwise
+    r"""Plots the correlation function for pairs of one-dimensional 'slices', otherwise
     referred to as the two point correlator at zero spatial momentum, as a function of
     time. Points and errorbars are means and standard deviations across a boostrap
     sample.
@@ -86,7 +86,7 @@ def plot_effective_pole_mass(training_geometry, effective_pole_mass):
         y=effective_pole_mass.mean(axis=-1),
         yerr=effective_pole_mass.std(axis=-1),
     )
-    ax.set_ylabel("$m_p^\mathrm{eff}$")
+    ax.set_ylabel(r"$m_p^\mathrm{eff}$")
     ax.set_xlabel("$x_2$")
     ax.set_title("Effective pole mass")
     return fig
@@ -282,7 +282,7 @@ def plot_magnetization_autocorr(
     fig, ax = plt.subplots()
     ax.set_title("Autocorrelation of magnetization")
     ax.set_ylabel(r"$\Gamma_M(\delta t)$")
-    ax.set_xlabel("$\delta t$")
+    ax.set_xlabel(r"$\delta t$")
 
     ax.plot(chain_indices, magnetization_autocorr[:cut])
 

From 7504922831fb09268493f6893feb8a465a73b157 Mon Sep 17 00:00:00 2001
From: marshrossney <17361029+marshrossney@users.noreply.github.com>
Date: Thu, 20 May 2021 16:29:06 +0100
Subject: [PATCH 10/17] minor change to counting number of successful cosh fits

---
 anvil/observables.py | 15 +++++++--------
 1 file changed, 7 insertions(+), 8 deletions(-)

diff --git a/anvil/observables.py b/anvil/observables.py
index a7e580f..25d8068 100644
--- a/anvil/observables.py
+++ b/anvil/observables.py
@@ -38,14 +38,13 @@ def fit_zero_momentum_correlator(
             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."
-        )
-        if failures >= n_boot - 1:
-            log.warning("Too many failures: no fit parameters will be returned.")
-            return None
+    n_fits = len(optimised_parameters)
+    log.info(
+        f"Cosh fit succeeded for {n_fits}/{n_boot} members of the bootstrap ensemble."
+    )
+    if n_fits < 2:
+        log.warning("Too few successful fits: no fit parameters will be returned.")
+        return None
 
     xi, A, c = np.array(optimised_parameters).transpose()
     return xi, A, c

From 7bb45a8db47d3fafb116df764c21a236ae881210 Mon Sep 17 00:00:00 2001
From: marshrossney <17361029+marshrossney@users.noreply.github.com>
Date: Thu, 20 May 2021 16:43:30 +0100
Subject: [PATCH 11/17] remove default cosh_fit_window in plot.py

---
 anvil/observables.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/anvil/observables.py b/anvil/observables.py
index 25d8068..b1259a4 100644
--- a/anvil/observables.py
+++ b/anvil/observables.py
@@ -20,7 +20,7 @@ def cosh_shift(x, xi, A, c):
 
 
 def fit_zero_momentum_correlator(
-    zero_momentum_correlator, training_geometry, cosh_fit_window=slice(1, None)
+    zero_momentum_correlator, training_geometry, cosh_fit_window
 ):
     t = np.arange(training_geometry.length) - training_geometry.length // 2
 

From c0a1df8c68c567fae2ebf5e2cef9fc81f6df5455 Mon Sep 17 00:00:00 2001
From: marshrossney <17361029+marshrossney@users.noreply.github.com>
Date: Thu, 20 May 2021 17:40:58 +0100
Subject: [PATCH 12/17] option to plot cosh fit overlay, add docstrings

---
 anvil/observables.py | 79 ++++++++++++++++++++++++++++++++++----------
 anvil/plot.py        |  5 +--
 anvil/table.py       | 31 ++++++++---------
 3 files changed, 78 insertions(+), 37 deletions(-)

diff --git a/anvil/observables.py b/anvil/observables.py
index b1259a4..6c7f7fb 100644
--- a/anvil/observables.py
+++ b/anvil/observables.py
@@ -5,12 +5,11 @@
 """
 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__)
 
@@ -22,10 +21,44 @@ def cosh_shift(x, xi, A, c):
 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: torch.Tensor
+        The two point correlation function at zero spatial momentum, i.e. the
+        correlation between 1-d 'slices'.
+    training_geometry: anvil.Geometry.
+        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
-    optimised_parameters = []
+    xi, A, c = [], [], []
     for correlator in zero_momentum_correlator.transpose():
         try:
             popt, pcov = optim.curve_fit(
@@ -33,26 +66,36 @@ def fit_zero_momentum_correlator(
                 xdata=t[cosh_fit_window],
                 ydata=correlator[cosh_fit_window],
             )
-            optimised_parameters.append(popt)
+            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(optimised_parameters)
+    n_fits = len(xi)
     log.info(
         f"Cosh fit succeeded for {n_fits}/{n_boot} members of the bootstrap ensemble."
     )
-    if n_fits < 2:
-        log.warning("Too few successful fits: 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):
+    """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 xi
+    return np.array(xi)
+
+
+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):
@@ -109,12 +152,12 @@ 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):
@@ -182,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):
@@ -238,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)
diff --git a/anvil/plot.py b/anvil/plot.py
index fa7f1c5..3b84444 100644
--- a/anvil/plot.py
+++ b/anvil/plot.py
@@ -25,13 +25,14 @@ def plot_zero_momentum_correlator(
     training_geometry,
     fit_zero_momentum_correlator,
     cosh_fit_window,
+    plot_cosh_fit=True,
 ):
     r"""Plots the correlation function for pairs of one-dimensional 'slices', otherwise
     referred to as the two point correlator at zero spatial momentum, as a function of
     time. Points and errorbars are means and standard deviations across a boostrap
     sample.
 
-    Also plots the :math:`1\sigma` confidence interval for a pure-exponential (cosh)
+    Optionally plots a :math:`1\sigma` confidence interval for a pure-exponential (cosh)
     fit performed for each member of the bootstrap sample in
     :py:func:`fit_zero_momentum_correlator`.
     """
@@ -46,7 +47,7 @@ def plot_zero_momentum_correlator(
         label="sample statistics",
     )
 
-    if fit_zero_momentum_correlator is not None:
+    if plot_cosh_fit:
         t = np.arange(training_geometry.length)[cosh_fit_window]
         fit = []
         for xi, A, c in zip(*fit_zero_momentum_correlator):
diff --git a/anvil/table.py b/anvil/table.py
index 7a9eda4..e2daab0 100644
--- a/anvil/table.py
+++ b/anvil/table.py
@@ -30,20 +30,17 @@ def table_autocorrelation(
 
 
 @table
-def table_fit(fit_zero_momentum_correlator, training_geometry):
-    if fit_zero_momentum_correlator is not None:
-        xi, A, c = fit_zero_momentum_correlator
-
-        res = [
-            [xi.mean(), xi.std()],
-            [c.mean(), c.std()],
-        ]
-        df = pd.DataFrame(
-            res,
-            columns=["mean", "error"],
-            index=["xi_from_fit", "magnetization_from_fit"],
-        )
-        return df
+def table_fit(correlation_length_from_fit, abs_magnetization_sq_from_fit):
+    res = [
+        [correlation_length_from_fit.mean(), correlation_length_from_fit.std()],
+        [abs_magnetization_sq_from_fit.mean(), abs_magnetization_sq_from_fit.std()],
+    ]
+    df = pd.DataFrame(
+        res,
+        columns=["mean", "error"],
+        index=["xi_from_fit", "abs_magnetization_sq_from_fit"],
+    )
+    return df
 
 
 @table
@@ -64,15 +61,15 @@ def table_two_point_scalars(ising_energy, susceptibility):
 
 
 @table
-def table_magnetization(abs_magnetization_squared, magnetic_susceptibility):
+def table_magnetization(abs_magnetization_sq, magnetic_susceptibility):
     res = [
-        [abs_magnetization_squared.mean(), abs_magnetization_squared.std()],
+        [abs_magnetization_sq.mean(), abs_magnetization_sq.std()],
         [magnetic_susceptibility.mean(), magnetic_susceptibility.std()],
     ]
     df = pd.DataFrame(
         res,
         columns=["mean", "error"],
-        index=["abs_magnetization_squared", "magnetic_susceptibility"],
+        index=["abs_magnetization_sq", "magnetic_susceptibility"],
     )
     return df
 

From 691a2829d8c524238921fad4ffcebc0fc851c276 Mon Sep 17 00:00:00 2001
From: marshrossney <17361029+marshrossney@users.noreply.github.com>
Date: Thu, 20 May 2021 17:46:44 +0100
Subject: [PATCH 13/17] minor typos in docstring

---
 anvil/observables.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/anvil/observables.py b/anvil/observables.py
index 6c7f7fb..c9657cf 100644
--- a/anvil/observables.py
+++ b/anvil/observables.py
@@ -30,11 +30,11 @@ def fit_zero_momentum_correlator(
 
     Parameters
     ----------
-    zero_momentum_correlator: torch.Tensor
+    zero_momentum_correlator
         The two point correlation function at zero spatial momentum, i.e. the
         correlation between 1-d 'slices'.
-    training_geometry: anvil.Geometry.
-        The lattice.
+    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

From a1850781aa622aea3b53026ec134ec2dd4beff3b Mon Sep 17 00:00:00 2001
From: marshrossney <17361029+marshrossney@users.noreply.github.com>
Date: Mon, 24 May 2021 11:16:06 +0100
Subject: [PATCH 14/17] update plot docstrings

---
 anvil/plot.py | 199 +++++++++++++++++++++++++++++++++++++++++++++-----
 1 file changed, 180 insertions(+), 19 deletions(-)

diff --git a/anvil/plot.py b/anvil/plot.py
index 3b84444..bc1aefa 100644
--- a/anvil/plot.py
+++ b/anvil/plot.py
@@ -35,6 +35,30 @@ def plot_zero_momentum_correlator(
     Optionally plots a :math:`1\sigma` confidence interval for a pure-exponential (cosh)
     fit performed for each member of the bootstrap sample in
     :py:func:`fit_zero_momentum_correlator`.
+
+    Parameters
+    ---------
+    zero_momentum_correlator
+        Array containing bootstrapped correlation function
+    training_geometry
+        Geometry object defining the lattice
+    fit_zero_momentum_correlator
+        The parameters resulting from a least-squares fit of a cosh function to the
+        correlator
+    cosh_fit_window
+        Slice object which indexes the lattice separations that were used to perform
+        the cosh fit to the correlation function
+    plot_cosh_fit
+        If False, only plot the correlation function, and not the result of the fit
+
+    Returns
+    -------
+    matplotlib.figure.Figure
+
+    See Also
+    --------
+    :py:func:`anvil.observables.zero_momentum_correlator`
+    :py:func:`anvil.observables.fit_zero_momentum_correlator`
     """
     fig, ax = plt.subplots()
 
@@ -80,6 +104,22 @@ def plot_zero_momentum_correlator(
 def plot_effective_pole_mass(training_geometry, effective_pole_mass):
     """Plots the (effective) pole mass as a function of 'time' separation. Points and
     error bars are means and standard deviations across a bootstrap sample.
+
+    Parameters
+    ----------
+    training_geometry
+        Geometry object defining the lattice.
+    effective_pole_mass
+        Array containing bootstrap ensemble of effective pole mass, for each
+        separation :math:`t = 1, \ldots, T - 1`
+
+    Returns
+    -------
+    matplotlib.figure.Figure
+
+    See Also
+    --------
+    :py:func:`anvil.observables.effective_pole_mass`
     """
     fig, ax = plt.subplots()
     ax.errorbar(
@@ -99,19 +139,39 @@ def plot_correlation_length(
     low_momentum_correlation_length,
     correlation_length_from_fit,
 ):
-    """Plots three estimates of correlation length.
+    """Plots three estimates of correlation length. Points and error bars are means
+    and standard deviations taken across a bootstrap ensemble.
 
     These are:
         1. Estimate from fitting a cosh function to the correlation between
             1-dimensional slices, using py:func:`correlation_length_from_fit`
         2. Reciprocal of the effective pole mass estimator, using
-            :py:func:`effective_pole_mass`
+            :py:func:`effective_pole_mass` (evaluated at each separation, :math:`x_2`. 
         3. Low momentum estimate, using :py:func:`low_momentum_correlation_length`
 
-    2. is evaluated at a specific value of the separation, x_2.
-
-    Points and error bars are means and standard deviations taken across a bootstrap
-    sample.
+    
+    Parameters
+    ----------
+    effective_pole_mass
+        Array containing estimate of the effective pole mass, for each separation
+        and each member of the bootstrap ensemble
+    low_momentum_correlation_length
+        Array containing a low-momentum estimate of the correlation length for
+        each member of the bootstrap ensemble.
+    correlation_length_from_fit
+        Array containing an estimate of the correlation length from a cosh fit
+        to the correlation function, for each member of the bootstrap
+        ensemble.
+
+    Returns
+    -------
+    matplotlib.figure.Figure
+
+    See Also
+    --------
+    :py:func:`anvil.observables.effective_pole_mass`
+    :py:func:`anvil.observables.low_momentum_correlation_length`
+    :py:func:`anvil.observables.fit_zero_momentum_correlator`
     """
     xi_arcosh = np.reciprocal(effective_pole_mass)
 
@@ -157,10 +217,26 @@ def plot_correlation_length(
 
 @figure
 def plot_two_point_correlator(two_point_correlator):
-    """Represents the two point correlator as a heatmap. The data shown is the mean
-    of a bootstrap sample of correlation functions, and is normalised so that
-    G(0, 0) = 1. The colour axis is scaled using a symmetric log scale, with a linear
-    region spanning [-0.01, 0.01].
+    """Represents the two point correlator as a heatmap.
+
+    The data shown is the mean of a bootstrap sample of correlation functions, and is
+    normalised so that :math:`G(0, 0) = 1`. The colour axis is scaled using a symmetric
+    log scale, with a linear region spanning :math:`[-0.01, 0.01]`.
+    
+    Parameters
+    ----------
+    two_point_correlator
+        Array containing two point correlation function for each two-dimensional
+        separation :math:`(x_1, x_2)`, for each member of a bootstrap ensemble.
+
+    Returns
+    -------
+    matplotlib.figure.Figure
+
+    See Also
+    --------
+    :py:func:`anvil.observables.two_point_correlator`
+    :py:func:`anvil.plot.plot_two_point_correlator_error`
     """
     corr = two_point_correlator.mean(axis=-1)
     corr /= corr[0, 0]
@@ -189,10 +265,26 @@ def plot_two_point_correlator(two_point_correlator):
 
 @figure
 def plot_two_point_correlator_error(two_point_correlator):
-    """Heatmap of the error in the two point correlator for each separation (x_1, x_2).
-    The error is computed as the standard deviation over the bootstrap sample. The data
-    shown is this error divided by the mean of the bootstrap sample, i.e. the
+    """Heatmap of the error in the two point correlator for each separation.
+
+    The error is computed as the standard deviation over the bootstrap sample. The
+    data shown is this error divided by the mean of the bootstrap sample, i.e. the
     fractional error.
+    
+    Parameters
+    ----------
+    two_point_correlator
+        Array containing two point correlation function for each two-dimensional
+        separation :math:`(x_1, x_2)`, for each member of a bootstrap ensemble.
+
+    Returns
+    -------
+    matplotlib.figure.Figure
+
+    See Also
+    --------
+    :py:func:`anvil.observables.two_point_correlator`
+    :py:func:`anvil.plot.plot_two_point_correlator`
     """
     corr = two_point_correlator.mean(axis=-1)
     error = two_point_correlator.std(axis=-1)
@@ -224,10 +316,19 @@ def plot_magnetization(magnetization_series):
     """Plots a histogram of the magnetization of each configuration in the Markov
     chain resulting from the Metropolis-Hastings sampling phase.
 
+    Parameters
+    ----------
+    magnetization_series
+        Array containing the magnetization for each configuration in the output
+        sample from the Metropolis-Hastings sampling phase.
+
+    Returns
+    -------
+    matplotlib.figure.Figure
+
     See also
     --------
-    :py:func:`plot_magnetization_series` displays the same data as a sequence,
-    i.e. in order of their appearance in the Markov chain.
+    :py:func:`plot_magnetization_series`
     """
     fig, ax = plt.subplots()
     ax.set_title("Magnetization")
@@ -241,10 +342,23 @@ def plot_magnetization(magnetization_series):
 def plot_magnetization_series(magnetization_series, sample_interval):
     """Plots the magnetization of each configuration in the Markov chain over the
     course of the Metropolis-Hastings sampling phase.
+    
+    Parameters
+    ----------
+    magnetization_series
+        Array containing the magnetization for each configuration in the output
+        sample from the Metropolis-Hastings sampling phase.
+    sample_interval
+        The number of Metropolis updates which were discarded between each
+        configuration appearing in the input series.
+
+    Returns
+    -------
+    matplotlib.figure.Figure
 
     See also
     --------
-    :py:func:`plot_magnetization` displays the same data as a histogram.
+    :py:func:`plot_magnetization`.
     """
     n_rows = 5
     if magnetization_series.size % n_rows != 0:
@@ -273,9 +387,32 @@ def plot_magnetization_autocorr(
     magnetization_autocorr, magnetization_optimal_window, sample_interval
 ):
     """Plots the autocorrelation function for the magnetization of the sequence of
-    configurations generated in the Metropolis-Hastings sampling phase. The x-axis
-    corresponds to a number of steps separating pairs of configurations in the
-    sequence.
+    configurations generated in the Metropolis-Hastings sampling phase.
+
+    The x-axis corresponds to a number of steps separating pairs of configurations
+    in the sequence.
+    
+    Parameters
+    ----------
+    magnetization_autocorr
+        Array containing the autocorrelation function of the magnetization for each
+        configuration in the output sample from the Metropolis-Hastings sampling phase.
+    magnetization_optimal_window
+        The size of the window in which the integrated autocorrelation time is to be
+        computed such that the total error is minimized.
+    sample_interval
+        The number of Metropolis updates which were discarded between each
+        configuration appearing in the input series.
+
+    Returns
+    -------
+    matplotlib.figure.Figure
+
+    See also
+    --------
+    :py:func:`anvil.observables.magnetization_autocorr`
+    :py:func:`anvil.observables.magnetization_optimal_window`
+    :py:func:`anvil.plot.plot_magnetization_integrated_autocorr`.
     """
     cut = max(10, 2 * magnetization_optimal_window)
     chain_indices = np.arange(cut) * sample_interval
@@ -319,8 +456,32 @@ def plot_magnetization_integrated_autocorr(
 ):
     """Plots the integrated autocorrelation function for the magnetization of the
     sequence of configurations generated in the Metropolis-Hastings sampling phase.
+    
     The x axis represents the size of the 'window' in which the summation is performed,
     i.e. the point at which the autocorrelation function is truncated.
+    
+    Parameters
+    ----------
+    magnetization_integrated_autocorr
+        Array containing the cumulative sum of the autocorrelation function of the
+        magnetization for each configuration in the output sample from the Metropolis-
+        Hastings sampling phase.
+    magnetization_optimal_window
+        The size of the window in which the integrated autocorrelation time is to be
+        computed such that the total error is minimized.
+    sample_interval
+        The number of Metropolis updates which were discarded between each
+        configuration appearing in the input series.
+
+    Returns
+    -------
+    matplotlib.figure.Figure
+
+    See also
+    --------
+    :py:func:`anvil.observables.magnetization_integrated_autocorr`
+    :py:func:`anvil.observables.magnetization_optimal_window`
+    :py:func:`anvil.plot.plot_magnetization_autocorr`.
     """
     cut = max(10, 2 * np.max(magnetization_optimal_window))
     chain_indices = np.arange(cut) * sample_interval

From 64ea16e1ab78aa566f907912ec9127975776f18f Mon Sep 17 00:00:00 2001
From: marshrossney <17361029+marshrossney@users.noreply.github.com>
Date: Mon, 24 May 2021 12:08:37 +0100
Subject: [PATCH 15/17] update table docstrings, minor alterations to plot
 docstrings

---
 anvil/plot.py  |  38 +++++-----
 anvil/table.py | 192 ++++++++++++++++++++++++++++++++++++++++++++++---
 2 files changed, 199 insertions(+), 31 deletions(-)

diff --git a/anvil/plot.py b/anvil/plot.py
index bc1aefa..cd0a6af 100644
--- a/anvil/plot.py
+++ b/anvil/plot.py
@@ -29,8 +29,9 @@ def plot_zero_momentum_correlator(
 ):
     r"""Plots the correlation function for pairs of one-dimensional 'slices', otherwise
     referred to as the two point correlator at zero spatial momentum, as a function of
-    time. Points and errorbars are means and standard deviations across a boostrap
-    sample.
+    time.
+
+    Points and errorbars are means and standard deviations across a boostrap ensemble.
 
     Optionally plots a :math:`1\sigma` confidence interval for a pure-exponential (cosh)
     fit performed for each member of the bootstrap sample in
@@ -57,8 +58,7 @@ def plot_zero_momentum_correlator(
 
     See Also
     --------
-    :py:func:`anvil.observables.zero_momentum_correlator`
-    :py:func:`anvil.observables.fit_zero_momentum_correlator`
+    :py:func:`anvil.table.table_zero_momentum_correlator`
     """
     fig, ax = plt.subplots()
 
@@ -102,8 +102,9 @@ def plot_zero_momentum_correlator(
 
 @figure
 def plot_effective_pole_mass(training_geometry, effective_pole_mass):
-    """Plots the (effective) pole mass as a function of 'time' separation. Points and
-    error bars are means and standard deviations across a bootstrap sample.
+    """Plots the (effective) pole mass as a function of 'time' separation.
+    
+    Points and errorbars are means and standard deviations across a boostrap ensemble.
 
     Parameters
     ----------
@@ -119,7 +120,7 @@ def plot_effective_pole_mass(training_geometry, effective_pole_mass):
 
     See Also
     --------
-    :py:func:`anvil.observables.effective_pole_mass`
+    :py:func:`anvil.table.table_effective_pole_mass`
     """
     fig, ax = plt.subplots()
     ax.errorbar(
@@ -139,9 +140,8 @@ def plot_correlation_length(
     low_momentum_correlation_length,
     correlation_length_from_fit,
 ):
-    """Plots three estimates of correlation length. Points and error bars are means
-    and standard deviations taken across a bootstrap ensemble.
-
+    """Plots three estimates of correlation length.
+    
     These are:
         1. Estimate from fitting a cosh function to the correlation between
             1-dimensional slices, using py:func:`correlation_length_from_fit`
@@ -149,6 +149,7 @@ def plot_correlation_length(
             :py:func:`effective_pole_mass` (evaluated at each separation, :math:`x_2`. 
         3. Low momentum estimate, using :py:func:`low_momentum_correlation_length`
 
+    Points and errorbars are means and standard deviations across a boostrap ensemble.
     
     Parameters
     ----------
@@ -169,9 +170,8 @@ def plot_correlation_length(
 
     See Also
     --------
-    :py:func:`anvil.observables.effective_pole_mass`
-    :py:func:`anvil.observables.low_momentum_correlation_length`
     :py:func:`anvil.observables.fit_zero_momentum_correlator`
+    :py:func:`anvil.table.table_correlation_length`
     """
     xi_arcosh = np.reciprocal(effective_pole_mass)
 
@@ -235,8 +235,8 @@ def plot_two_point_correlator(two_point_correlator):
 
     See Also
     --------
-    :py:func:`anvil.observables.two_point_correlator`
     :py:func:`anvil.plot.plot_two_point_correlator_error`
+    :py:func:`anvil.table.table_two_point_correlator`
     """
     corr = two_point_correlator.mean(axis=-1)
     corr /= corr[0, 0]
@@ -283,8 +283,8 @@ def plot_two_point_correlator_error(two_point_correlator):
 
     See Also
     --------
-    :py:func:`anvil.observables.two_point_correlator`
     :py:func:`anvil.plot.plot_two_point_correlator`
+    :py:func:`anvil.table.table_two_point_correlator`
     """
     corr = two_point_correlator.mean(axis=-1)
     error = two_point_correlator.std(axis=-1)
@@ -328,7 +328,8 @@ def plot_magnetization(magnetization_series):
 
     See also
     --------
-    :py:func:`plot_magnetization_series`
+    :py:func:`anvil.plot.plot_magnetization_series`
+    :py:func:`anvil.table.table_magnetization`
     """
     fig, ax = plt.subplots()
     ax.set_title("Magnetization")
@@ -358,7 +359,8 @@ def plot_magnetization_series(magnetization_series, sample_interval):
 
     See also
     --------
-    :py:func:`plot_magnetization`.
+    :py:func:`anvil.plot.plot_magnetization`.
+    :py:func:`anvil.table.table_magnetization`
     """
     n_rows = 5
     if magnetization_series.size % n_rows != 0:
@@ -410,8 +412,6 @@ def plot_magnetization_autocorr(
 
     See also
     --------
-    :py:func:`anvil.observables.magnetization_autocorr`
-    :py:func:`anvil.observables.magnetization_optimal_window`
     :py:func:`anvil.plot.plot_magnetization_integrated_autocorr`.
     """
     cut = max(10, 2 * magnetization_optimal_window)
@@ -479,8 +479,6 @@ def plot_magnetization_integrated_autocorr(
 
     See also
     --------
-    :py:func:`anvil.observables.magnetization_integrated_autocorr`
-    :py:func:`anvil.observables.magnetization_optimal_window`
     :py:func:`anvil.plot.plot_magnetization_autocorr`.
     """
     cut = max(10, 2 * np.max(magnetization_optimal_window))
diff --git a/anvil/table.py b/anvil/table.py
index e2daab0..fab78a4 100644
--- a/anvil/table.py
+++ b/anvil/table.py
@@ -19,6 +19,34 @@ def table_autocorrelation(
     tau_chain,
     acceptance,
 ):
+    """
+    Tabulate some information related to the statistical efficiency of the Metropolis-
+    Hastings sampling phase.
+
+    Parameters
+    ----------
+    magnetization_integrated_autocorr
+        Array containing the cumulative sum of the autocorrelation function of the
+        magnetization for each configuration in the sample output by the Metropolis-
+        Hastings sampling phase.
+    magnetization_optimal_window
+        Integer corresponding to a window size in which the autocorrelation function
+        should be summed, such that the resulting estimate of the integrated
+        autocorrelation has the smallest possible total error.
+    tau_chain
+        Estimate of the integrated autocorrelation using the accept-reject statistics
+        of the sampling phase.
+    acceptance
+        Fraction of proposals which were accepted in the sampling phase.
+
+    Returns
+    -------
+    pandas.core.frame.DataFrame
+
+    See Also
+    --------
+    :py:func:`anvil.sample.calc_tau_chain`
+    """
     tau_mag = magnetization_integrated_autocorr[magnetization_optimal_window]
 
     df = pd.DataFrame(
@@ -31,6 +59,24 @@ def table_autocorrelation(
 
 @table
 def table_fit(correlation_length_from_fit, abs_magnetization_sq_from_fit):
+    """Tabulate the correlation length and magnetization estimates resulting from the
+    fitting of a cosh to the correlation function.
+    
+    Values and errors are means and standard deviations over a bootstrap ensemble.
+
+    Parameters
+    ----------
+    correlation_length_from_fit
+    abs_magnetization_sq_from_fit
+
+    Returns
+    -------
+    pandas.core.frame.DataFrame
+
+    See Also
+    --------
+    :py:func:`anvil.observables.fit_zero_momentum_correlator`
+    """
     res = [
         [correlation_length_from_fit.mean(), correlation_length_from_fit.std()],
         [abs_magnetization_sq_from_fit.mean(), abs_magnetization_sq_from_fit.std()],
@@ -45,8 +91,25 @@ def table_fit(correlation_length_from_fit, abs_magnetization_sq_from_fit):
 
 @table
 def table_two_point_scalars(ising_energy, susceptibility):
-    """Table of the ising observables, with mean and standard deviation taken
-    across boostrap samples
+    """Table of scalar observables derived from the two point correlation function.
+    
+    Values and errors are means and standard deviations over a bootstrap ensemble.
+
+    Parameters
+    ----------
+    ising_energy
+        Nearest-neighbour iteraction energy.
+    susceptibility
+        Magnetic susceptibility defined by the sum of the correlation function.
+
+    Returns
+    -------
+    pandas.core.frame.DataFrame
+
+    See Also
+    --------
+    :py:func:`anvil.observables.ising_energy`
+    :py:func:`anvil.observables.susceptibility`
     """
     res = [
         [ising_energy.mean(), ising_energy.std()],
@@ -62,6 +125,28 @@ def table_two_point_scalars(ising_energy, susceptibility):
 
 @table
 def table_magnetization(abs_magnetization_sq, magnetic_susceptibility):
+    """Table containing quantities derived from the sample-averaged magnetization.
+
+    Values and errors are means and standard deviations over a bootstrap ensemble.
+
+    Parameters
+    ----------
+    abs_magnetization_sq
+        Array containing the sample mean of the absolute magnetization, squared, for each
+        member of the bootstrap ensemble.
+    magnetic_susceptibility
+        Array containing the susceptibility for each member of the bootstrap ensemble.
+
+    Returns
+    -------
+    pandas.core.frame.DataFrame
+
+    See Also
+    --------
+    :py:func:`anvil.tables.table_two_point_scalars`
+    :py:func:`anvil.tables.table_fit`
+    """
+
     res = [
         [abs_magnetization_sq.mean(), abs_magnetization_sq.std()],
         [magnetic_susceptibility.mean(), magnetic_susceptibility.std()],
@@ -82,8 +167,39 @@ def table_correlation_length(
     correlation_length_from_fit,
     training_geometry,
 ):
-    """Tabulate four estimates of correlation length, with values and errors
-    taken as the mean and standard deviation of the bootstrap sample.
+    """Table containing four estimates of correlation length.
+
+    Values and errors are means and standard deviations over a bootstrap ensemble.
+
+    Also displays the number of correlation lengths that can fit on the lattice, i.e.
+    :math:`\xi / L` where :math:`\xi` is the correlation length and :math:`L` is the
+    linear extent of the lattice.
+
+    Parameters
+    ----------
+    effective_pole_mass
+        Array containing estimate of the effective pole mass, for each separation
+        and each member of the bootstrap ensemble
+    second_moment_correlation_length
+        Estimate of the correlation length based on the second moment of the
+        two point correlation function, for each member of the bootstrap ensemble.
+    low_momentum_correlation_length
+        Array containing a low-momentum estimate of the correlation length for each
+        member of the bootstrap ensemble.
+    correlation_length_from_fit
+        Array containing an estimate of the correlation length from a cosh fit to
+        the correlation function, for each member of the bootstrap ensemble.
+    training_geometry
+        Geometry object defining the lattice.
+
+    Returns
+    -------
+    pandas.core.frame.DataFrame
+
+    See Also
+    --------
+    :py:func:`anvil.observables.fit_zero_momentum_correlator`
+    :py:func:`anvil.plot.plot_correlation_length`
     """
     # Take the mean of the arcosh estimator over "large" separations
     x0 = training_geometry.length // 4
@@ -119,8 +235,26 @@ def table_correlation_length(
 
 @table
 def table_zero_momentum_correlator(zero_momentum_correlator, training_geometry):
-    """Table of zero_momentum_correlator, with mean and standard deviation
-    from bootstrap
+    """Table containing values of the two point correlation function in time-momentum
+    representation at zero momentum, for each separation.
+
+    Values and errors are means and standard deviations over a bootstrap ensemble.
+
+    Parameters
+    ----------
+    zero_momentum_correlator
+        Array containing the correlation function for each 1-d separation, for each
+        member of the bootstrap ensemble.
+    training_geometry
+        Geometry object defining the lattice.
+
+    Returns
+    -------
+    pandas.core.frame.DataFrame
+
+    See Also
+    --------
+    :py:func:`anvil.plot.plot_zero_momentum_correlator`
     """
     means = zero_momentum_correlator.mean(axis=-1)[:, np.newaxis]
     stds = zero_momentum_correlator.std(axis=-1)[:, np.newaxis]
@@ -137,8 +271,25 @@ def table_zero_momentum_correlator(zero_momentum_correlator, training_geometry):
 
 @table
 def table_effective_pole_mass(effective_pole_mass, training_geometry):
-    """Table of effective_pole_mass, with mean and standard deviation
-    from bootstrap
+    """Table containing values of the effective pole mass for each separation.
+
+    Values and errors are means and standard deviations over a bootstrap ensemble.
+
+    Parameters
+    ----------
+    effective_pole_mass
+        Array containing the effective pole mass for each separation, for each
+        member of the bootstrap ensemble.
+    training_geometry
+        Geometry object defining the lattice.
+
+    Returns
+    -------
+    pandas.core.frame.DataFrame
+
+    See Also
+    --------
+    :py:func:`anvil.plot.plot_effective_pole_mass`
     """
     means = effective_pole_mass.mean(axis=-1)[:, np.newaxis]
     stds = effective_pole_mass.std(axis=-1)[:, np.newaxis]
@@ -153,9 +304,28 @@ def table_effective_pole_mass(effective_pole_mass, training_geometry):
 
 
 @table
-def table_two_point_correlator(training_geometry, two_point_correlator):
-    """For each x and t, tabulate the mean and standard deviation of the two
-    point function, estimated from bootstrap sample
+def table_two_point_correlator(two_point_correlator, training_geometry):
+    """Table containing values of the two point correlation function for each
+    two-dimensional separation.
+
+    Values and errors are means and standard deviations over a bootstrap ensemble.
+
+    Parameters
+    ----------
+    two_point_correlator
+        Array containing the correlation function for each 2-d separation, for each
+        member of the bootstrap ensemble.
+    training_geometry
+        Geometry object defining the lattice.
+
+    Returns
+    -------
+    pandas.core.frame.DataFrame
+
+    See Also
+    --------
+    :py:func:`anvil.plot.plot_two_point_correlator`
+    :py:func:`anvil.plot.plot_two_point_correlator_error`
     """
     corr = []
     index = []

From 35536e4c206c8c3c40fd0175afcc6c8fb12a6d60 Mon Sep 17 00:00:00 2001
From: marshrossney <17361029+marshrossney@users.noreply.github.com>
Date: Mon, 24 May 2021 12:16:24 +0100
Subject: [PATCH 16/17] indicate that bootstrap dim is assumed to be last

---
 anvil/plot.py  | 13 ++++++++++---
 anvil/table.py | 21 ++++++++++++++-------
 2 files changed, 24 insertions(+), 10 deletions(-)

diff --git a/anvil/plot.py b/anvil/plot.py
index cd0a6af..75e27a8 100644
--- a/anvil/plot.py
+++ b/anvil/plot.py
@@ -31,7 +31,8 @@ def plot_zero_momentum_correlator(
     referred to as the two point correlator at zero spatial momentum, as a function of
     time.
 
-    Points and errorbars are means and standard deviations across a boostrap ensemble.
+    Points and errorbars are means and standard deviations across a boostrap ensemble,
+    which is assumed to be the last (``-1``) dimension of input arrays.
 
     Optionally plots a :math:`1\sigma` confidence interval for a pure-exponential (cosh)
     fit performed for each member of the bootstrap sample in
@@ -104,7 +105,8 @@ def plot_zero_momentum_correlator(
 def plot_effective_pole_mass(training_geometry, effective_pole_mass):
     """Plots the (effective) pole mass as a function of 'time' separation.
     
-    Points and errorbars are means and standard deviations across a boostrap ensemble.
+    Points and errorbars are means and standard deviations across a boostrap ensemble,
+    which is assumed to be the last (``-1``) dimension of input arrays.
 
     Parameters
     ----------
@@ -149,7 +151,8 @@ def plot_correlation_length(
             :py:func:`effective_pole_mass` (evaluated at each separation, :math:`x_2`. 
         3. Low momentum estimate, using :py:func:`low_momentum_correlation_length`
 
-    Points and errorbars are means and standard deviations across a boostrap ensemble.
+    Points and errorbars are means and standard deviations across a boostrap ensemble,
+    which is assumed to be the last (``-1``) dimension of input arrays.
     
     Parameters
     ----------
@@ -223,6 +226,8 @@ def plot_two_point_correlator(two_point_correlator):
     normalised so that :math:`G(0, 0) = 1`. The colour axis is scaled using a symmetric
     log scale, with a linear region spanning :math:`[-0.01, 0.01]`.
     
+    The bootstrap dimension is assumed to be the last (``-1``) dimension of input arrays.
+    
     Parameters
     ----------
     two_point_correlator
@@ -271,6 +276,8 @@ def plot_two_point_correlator_error(two_point_correlator):
     data shown is this error divided by the mean of the bootstrap sample, i.e. the
     fractional error.
     
+    The bootstrap dimension is assumed to be the last (``-1``) dimension of input arrays.
+    
     Parameters
     ----------
     two_point_correlator
diff --git a/anvil/table.py b/anvil/table.py
index fab78a4..a92fbc0 100644
--- a/anvil/table.py
+++ b/anvil/table.py
@@ -62,7 +62,8 @@ def table_fit(correlation_length_from_fit, abs_magnetization_sq_from_fit):
     """Tabulate the correlation length and magnetization estimates resulting from the
     fitting of a cosh to the correlation function.
     
-    Values and errors are means and standard deviations over a bootstrap ensemble.
+    Values and errors are means and standard deviations over a bootstrap ensemble,
+    which is assumed to be the last (``-1``) dimension of input arrays.
 
     Parameters
     ----------
@@ -93,7 +94,8 @@ def table_fit(correlation_length_from_fit, abs_magnetization_sq_from_fit):
 def table_two_point_scalars(ising_energy, susceptibility):
     """Table of scalar observables derived from the two point correlation function.
     
-    Values and errors are means and standard deviations over a bootstrap ensemble.
+    Values and errors are means and standard deviations over a bootstrap ensemble,
+    which is assumed to be the last (``-1``) dimension of input arrays.
 
     Parameters
     ----------
@@ -127,7 +129,8 @@ def table_two_point_scalars(ising_energy, susceptibility):
 def table_magnetization(abs_magnetization_sq, magnetic_susceptibility):
     """Table containing quantities derived from the sample-averaged magnetization.
 
-    Values and errors are means and standard deviations over a bootstrap ensemble.
+    Values and errors are means and standard deviations over a bootstrap ensemble,
+    which is assumed to be the last (``-1``) dimension of input arrays.
 
     Parameters
     ----------
@@ -169,7 +172,8 @@ def table_correlation_length(
 ):
     """Table containing four estimates of correlation length.
 
-    Values and errors are means and standard deviations over a bootstrap ensemble.
+    Values and errors are means and standard deviations over a bootstrap ensemble,
+    which is assumed to be the last (``-1``) dimension of input arrays.
 
     Also displays the number of correlation lengths that can fit on the lattice, i.e.
     :math:`\xi / L` where :math:`\xi` is the correlation length and :math:`L` is the
@@ -238,7 +242,8 @@ def table_zero_momentum_correlator(zero_momentum_correlator, training_geometry):
     """Table containing values of the two point correlation function in time-momentum
     representation at zero momentum, for each separation.
 
-    Values and errors are means and standard deviations over a bootstrap ensemble.
+    Values and errors are means and standard deviations over a bootstrap ensemble,
+    which is assumed to be the last (``-1``) dimension of input arrays.
 
     Parameters
     ----------
@@ -273,7 +278,8 @@ def table_zero_momentum_correlator(zero_momentum_correlator, training_geometry):
 def table_effective_pole_mass(effective_pole_mass, training_geometry):
     """Table containing values of the effective pole mass for each separation.
 
-    Values and errors are means and standard deviations over a bootstrap ensemble.
+    Values and errors are means and standard deviations over a bootstrap ensemble,
+    which is assumed to be the last (``-1``) dimension of input arrays.
 
     Parameters
     ----------
@@ -308,7 +314,8 @@ def table_two_point_correlator(two_point_correlator, training_geometry):
     """Table containing values of the two point correlation function for each
     two-dimensional separation.
 
-    Values and errors are means and standard deviations over a bootstrap ensemble.
+    Values and errors are means and standard deviations over a bootstrap ensemble
+    which is assumed to be the last (``-1``) dimension of input arrays.
 
     Parameters
     ----------

From 0980f87b07aa96a7736108c02627a2bed0a8e840 Mon Sep 17 00:00:00 2001
From: marshrossney <17361029+marshrossney@users.noreply.github.com>
Date: Mon, 24 May 2021 12:28:10 +0100
Subject: [PATCH 17/17] docstrings need to be raw strings

---
 anvil/plot.py  | 28 +++++++++++++++-------------
 anvil/table.py | 23 ++++++++++++-----------
 2 files changed, 27 insertions(+), 24 deletions(-)

diff --git a/anvil/plot.py b/anvil/plot.py
index 75e27a8..f6d1294 100644
--- a/anvil/plot.py
+++ b/anvil/plot.py
@@ -24,8 +24,8 @@ def plot_zero_momentum_correlator(
     zero_momentum_correlator,
     training_geometry,
     fit_zero_momentum_correlator,
-    cosh_fit_window,
-    plot_cosh_fit=True,
+    cosh_fit_window: type(slice),
+    plot_cosh_fit: bool = True,
 ):
     r"""Plots the correlation function for pairs of one-dimensional 'slices', otherwise
     referred to as the two point correlator at zero spatial momentum, as a function of
@@ -103,7 +103,7 @@ def plot_zero_momentum_correlator(
 
 @figure
 def plot_effective_pole_mass(training_geometry, effective_pole_mass):
-    """Plots the (effective) pole mass as a function of 'time' separation.
+    r"""Plots the (effective) pole mass as a function of 'time' separation.
     
     Points and errorbars are means and standard deviations across a boostrap ensemble,
     which is assumed to be the last (``-1``) dimension of input arrays.
@@ -142,7 +142,7 @@ def plot_correlation_length(
     low_momentum_correlation_length,
     correlation_length_from_fit,
 ):
-    """Plots three estimates of correlation length.
+    r"""Plots three estimates of correlation length.
     
     These are:
         1. Estimate from fitting a cosh function to the correlation between
@@ -220,7 +220,7 @@ def plot_correlation_length(
 
 @figure
 def plot_two_point_correlator(two_point_correlator):
-    """Represents the two point correlator as a heatmap.
+    r"""Represents the two point correlator as a heatmap.
 
     The data shown is the mean of a bootstrap sample of correlation functions, and is
     normalised so that :math:`G(0, 0) = 1`. The colour axis is scaled using a symmetric
@@ -270,7 +270,7 @@ def plot_two_point_correlator(two_point_correlator):
 
 @figure
 def plot_two_point_correlator_error(two_point_correlator):
-    """Heatmap of the error in the two point correlator for each separation.
+    r"""Heatmap of the error in the two point correlator for each separation.
 
     The error is computed as the standard deviation over the bootstrap sample. The
     data shown is this error divided by the mean of the bootstrap sample, i.e. the
@@ -320,7 +320,7 @@ def plot_two_point_correlator_error(two_point_correlator):
 
 @figure
 def plot_magnetization(magnetization_series):
-    """Plots a histogram of the magnetization of each configuration in the Markov
+    r"""Plots a histogram of the magnetization of each configuration in the Markov
     chain resulting from the Metropolis-Hastings sampling phase.
 
     Parameters
@@ -348,7 +348,7 @@ def plot_magnetization(magnetization_series):
 
 @figure
 def plot_magnetization_series(magnetization_series, sample_interval):
-    """Plots the magnetization of each configuration in the Markov chain over the
+    r"""Plots the magnetization of each configuration in the Markov chain over the
     course of the Metropolis-Hastings sampling phase.
     
     Parameters
@@ -393,9 +393,9 @@ def plot_magnetization_series(magnetization_series, sample_interval):
 
 @figure
 def plot_magnetization_autocorr(
-    magnetization_autocorr, magnetization_optimal_window, sample_interval
+        magnetization_autocorr, magnetization_optimal_window: int, sample_interval: int
 ):
-    """Plots the autocorrelation function for the magnetization of the sequence of
+    r"""Plots the autocorrelation function for the magnetization of the sequence of
     configurations generated in the Metropolis-Hastings sampling phase.
 
     The x-axis corresponds to a number of steps separating pairs of configurations
@@ -419,6 +419,7 @@ def plot_magnetization_autocorr(
 
     See also
     --------
+    :py:func:`anvil.observables.optimal_window`
     :py:func:`anvil.plot.plot_magnetization_integrated_autocorr`.
     """
     cut = max(10, 2 * magnetization_optimal_window)
@@ -458,10 +459,10 @@ def plot_magnetization_autocorr(
 @figure
 def plot_magnetization_integrated_autocorr(
     magnetization_integrated_autocorr,
-    magnetization_optimal_window,
-    sample_interval,
+    magnetization_optimal_window: int,
+    sample_interval: int,
 ):
-    """Plots the integrated autocorrelation function for the magnetization of the
+    r"""Plots the integrated autocorrelation function for the magnetization of the
     sequence of configurations generated in the Metropolis-Hastings sampling phase.
     
     The x axis represents the size of the 'window' in which the summation is performed,
@@ -486,6 +487,7 @@ def plot_magnetization_integrated_autocorr(
 
     See also
     --------
+    :py:func:`anvil.observables.optimal_window`
     :py:func:`anvil.plot.plot_magnetization_autocorr`.
     """
     cut = max(10, 2 * np.max(magnetization_optimal_window))
diff --git a/anvil/table.py b/anvil/table.py
index a92fbc0..2e873dc 100644
--- a/anvil/table.py
+++ b/anvil/table.py
@@ -15,11 +15,11 @@
 @table
 def table_autocorrelation(
     magnetization_integrated_autocorr,
-    magnetization_optimal_window,
-    tau_chain,
-    acceptance,
+    magnetization_optimal_window: int,
+    tau_chain: float,
+    acceptance: float,
 ):
-    """
+    r"""
     Tabulate some information related to the statistical efficiency of the Metropolis-
     Hastings sampling phase.
 
@@ -45,6 +45,7 @@ def table_autocorrelation(
 
     See Also
     --------
+    :py:func:`anvil.observables.optimal_window`
     :py:func:`anvil.sample.calc_tau_chain`
     """
     tau_mag = magnetization_integrated_autocorr[magnetization_optimal_window]
@@ -59,7 +60,7 @@ def table_autocorrelation(
 
 @table
 def table_fit(correlation_length_from_fit, abs_magnetization_sq_from_fit):
-    """Tabulate the correlation length and magnetization estimates resulting from the
+    r"""Tabulate the correlation length and magnetization estimates resulting from the
     fitting of a cosh to the correlation function.
     
     Values and errors are means and standard deviations over a bootstrap ensemble,
@@ -92,7 +93,7 @@ def table_fit(correlation_length_from_fit, abs_magnetization_sq_from_fit):
 
 @table
 def table_two_point_scalars(ising_energy, susceptibility):
-    """Table of scalar observables derived from the two point correlation function.
+    r"""Table of scalar observables derived from the two point correlation function.
     
     Values and errors are means and standard deviations over a bootstrap ensemble,
     which is assumed to be the last (``-1``) dimension of input arrays.
@@ -127,7 +128,7 @@ def table_two_point_scalars(ising_energy, susceptibility):
 
 @table
 def table_magnetization(abs_magnetization_sq, magnetic_susceptibility):
-    """Table containing quantities derived from the sample-averaged magnetization.
+    r"""Table containing quantities derived from the sample-averaged magnetization.
 
     Values and errors are means and standard deviations over a bootstrap ensemble,
     which is assumed to be the last (``-1``) dimension of input arrays.
@@ -170,7 +171,7 @@ def table_correlation_length(
     correlation_length_from_fit,
     training_geometry,
 ):
-    """Table containing four estimates of correlation length.
+    r"""Table containing four estimates of correlation length.
 
     Values and errors are means and standard deviations over a bootstrap ensemble,
     which is assumed to be the last (``-1``) dimension of input arrays.
@@ -239,7 +240,7 @@ def table_correlation_length(
 
 @table
 def table_zero_momentum_correlator(zero_momentum_correlator, training_geometry):
-    """Table containing values of the two point correlation function in time-momentum
+    r"""Table containing values of the two point correlation function in time-momentum
     representation at zero momentum, for each separation.
 
     Values and errors are means and standard deviations over a bootstrap ensemble,
@@ -276,7 +277,7 @@ def table_zero_momentum_correlator(zero_momentum_correlator, training_geometry):
 
 @table
 def table_effective_pole_mass(effective_pole_mass, training_geometry):
-    """Table containing values of the effective pole mass for each separation.
+    r"""Table containing values of the effective pole mass for each separation.
 
     Values and errors are means and standard deviations over a bootstrap ensemble,
     which is assumed to be the last (``-1``) dimension of input arrays.
@@ -311,7 +312,7 @@ def table_effective_pole_mass(effective_pole_mass, training_geometry):
 
 @table
 def table_two_point_correlator(two_point_correlator, training_geometry):
-    """Table containing values of the two point correlation function for each
+    r"""Table containing values of the two point correlation function for each
     two-dimensional separation.
 
     Values and errors are means and standard deviations over a bootstrap ensemble