From 03083b5b8f06792f1991c1428d5fdb9257914e13 Mon Sep 17 00:00:00 2001
From: Finlay Clark <finlay.clark@newcastle.ac.uk>
Date: Thu, 6 Aug 2026 11:58:59 +0100
Subject: [PATCH 1/3] docs: fix pep257 issues

---
 devtools/scripts/create_conda_env.py |  2 +-
 pyproject.toml                       | 10 ++-
 red/__init__.py                      |  2 +-
 red/_validation.py                   | 12 ++--
 red/confidence_intervals.py          |  7 +-
 red/equilibration.py                 | 56 ++++++++--------
 red/ess.py                           | 70 +++++++++++---------
 red/gelman_rubin.py                  | 18 ++---
 red/plot.py                          | 28 ++++----
 red/sse.py                           | 20 +++---
 red/variance.py                      | 98 ++++++++++++++++------------
 11 files changed, 180 insertions(+), 143 deletions(-)

diff --git a/devtools/scripts/create_conda_env.py b/devtools/scripts/create_conda_env.py
index e3087e5..229f893 100644
--- a/devtools/scripts/create_conda_env.py
+++ b/devtools/scripts/create_conda_env.py
@@ -49,7 +49,7 @@
 
 @contextmanager
 def temp_cd():
-    """Temporary CD Helper"""
+    """Temporarily change the working directory (context manager helper)."""
     cwd = os.getcwd()
     with TemporaryDirectory() as td:
         try:
diff --git a/pyproject.toml b/pyproject.toml
index c073b5b..e395425 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -74,10 +74,18 @@ line-length = 100
 [tool.ruff.lint]
 ignore = ["PLR", "PLW", "C901"]
 select = ["B","C","E","F","W","B9"]
+# Enforce PEP257 docstring conventions (pydocstyle). The numpy convention matches the
+# numpy-style docstrings used throughout the package; it disables D212/D413/D415/D417 by
+# default, so we re-enable those explicitly. D213 is intentionally left disabled by the
+# numpy convention (D212 and D213 are mutually exclusive).
+extend-select = ["D", "D212", "D413", "D415", "D417"]
+
+[tool.ruff.lint.pydocstyle]
+convention = "numpy"
 
 [tool.ruff.lint.per-file-ignores]
 "__init__.py" = ["F401"]
-"red/tests/*.py" = ["F401", "F811"]
+"red/tests/*.py" = ["F401", "F811", "D"]
 
 [tool.setuptools]
 # This subkey is a beta stage development and keys may change in the future, see https://setuptools.pypa.io/en/latest/userguide/pyproject_config.html for more details
diff --git a/red/__init__.py b/red/__init__.py
index 3852470..6b83bf8 100644
--- a/red/__init__.py
+++ b/red/__init__.py
@@ -1,4 +1,4 @@
-"""Robust Equilibration Detection"""
+"""Robust Equilibration Detection."""
 
 from ._version import __version__
 from .confidence_intervals import get_conf_int_init_seq
diff --git a/red/_validation.py b/red/_validation.py
index 9381ffd..6e5a0e4 100644
--- a/red/_validation.py
+++ b/red/_validation.py
@@ -11,12 +11,11 @@
 def check_data(
     data: _npt.NDArray[_np.float64], one_dim_allowed: bool = False
 ) -> _npt.NDArray[_np.float64]:
-    """
-    Assert that data passed is a numpy array where
-    the first dimension is the number of chains and
-    the second dimension is the number of samples.
-    If the array is one dimensional, add a second
-    dimension with length 1.
+    """Validate and reshape input data to 2D ``(n_chains, n_samples)``.
+
+    Asserts that the data passed is a numpy array where the first dimension is the number
+    of chains and the second dimension is the number of samples. If the array is one
+    dimensional, a second dimension with length 1 is added.
 
     Parameters
     ----------
@@ -30,6 +29,7 @@ def check_data(
     -------
     np.ndarray
         Data with shape (n_chains, n_samples).
+
     """
     # Check that data is a numpy array.
     if not isinstance(data, _np.ndarray):
diff --git a/red/confidence_intervals.py b/red/confidence_intervals.py
index b72c4c6..c578c51 100644
--- a/red/confidence_intervals.py
+++ b/red/confidence_intervals.py
@@ -17,9 +17,9 @@ def get_conf_int_init_seq(
     min_max_lag_time: int = 3,
     max_max_lag_time: _Optional[int] = None,
 ) -> float:
-    """
-    Calculate the confidence interval for the mean of a time
-    series using initial sequence methods. See Geyer, 1992:
+    """Calculate the confidence interval for the mean of a time series.
+
+    Uses initial sequence methods. See Geyer, 1992:
     https://www.jstor.org/stable/2246094.
 
     Parameters
@@ -48,6 +48,7 @@ def get_conf_int_init_seq(
     -------
     float
         The standard error of the mean.
+
     """
     # Get the correlated estimate of the variance.
     var_cor, max_lag, acovf = _get_variance_initial_sequence(
diff --git a/red/equilibration.py b/red/equilibration.py
index f71f7c3..cb83508 100644
--- a/red/equilibration.py
+++ b/red/equilibration.py
@@ -34,13 +34,12 @@ def detect_equilibration_init_seq(
     data_y_label: str = r"$\Delta G$ / kcal mol$^{-1}$",
     plot_max_lags: bool = True,
 ) -> _Tuple[_Union[float, int], float, float]:
-    r"""
-    Detect the equilibration time of a time series by finding the minimum
-    squared standard error (SSE), or maximum effective sample size (ESS)
-    of the time series, using initial sequence estimators of the variance.
-    This is done by computing the SSE at each time point, discarding all
-    samples before the time point. The index of the time point with
-    the minimum SSE or maximum ESS is taken to be the point of equilibration.
+    r"""Detect the equilibration time of a time series by finding the minimum SSE or maximum ESS.
+
+    The variance is estimated using initial sequence estimators. This is done by computing
+    the squared standard error (SSE) at each time point, discarding all samples before the
+    time point. The index of the time point with the minimum SSE, or maximum effective
+    sample size (ESS), is taken to be the point of equilibration.
 
     Parameters
     ----------
@@ -108,6 +107,7 @@ def detect_equilibration_init_seq(
 
     equil_ess: float
         The effective sample size at the equilibration point.
+
     """
     # Check that data is valid.
     data = check_data(data, one_dim_allowed=True)
@@ -201,13 +201,12 @@ def detect_equilibration_window(
     data_y_label: str = r"$\Delta G$ / kcal mol$^{-1}$",
     plot_window_size: bool = True,
 ) -> _Tuple[_Union[float, int], float, float]:
-    r"""
-    Detect the equilibration time of a time series by finding the minimum
-    squared standard error (SSE) or maximum effective sample size (ESS)
-    of the time series, using window estimators of the variance. This is
-    done by computing the SSE at each time point, discarding all samples
-    before the time point. The index of the time point with the minimum
-    SSE is taken to be the point of equilibration.
+    r"""Detect the equilibration time of a time series by finding the minimum SSE or maximum ESS.
+
+    The variance is estimated using window estimators. This is done by computing the
+    squared standard error (SSE) at each time point, discarding all samples before the time
+    point. The index of the time point with the minimum SSE, or maximum effective sample
+    size (ESS), is taken to be the point of equilibration.
 
     Parameters
     ----------
@@ -269,6 +268,7 @@ def detect_equilibration_window(
 
     equil_ess: float
         The effective sample size at the equilibration point.
+
     """
     # Check that data is valid.
     data = check_data(data, one_dim_allowed=True)
@@ -350,11 +350,11 @@ def get_paired_t_p_timeseries(
     final_block_size: float = 0.5,
     t_test_sidedness: str = "two-sided",
 ) -> _Tuple[_npt.NDArray[_np.float64], _npt.NDArray[_Union[_np.int64, _np.float64]]]:
-    """
-    Get a timeseries of the p-values from a paired t-test on the differences
-    between sample means between intial and final portions of the data. The timeseries
-    is obtained by repeatedly discarding more data from the time series between
-    calculations of the p-value.
+    """Get a timeseries of the p-values from a paired t-test.
+
+    The p-values come from a paired t-test on the differences between sample means between
+    initial and final portions of the data. The timeseries is obtained by repeatedly
+    discarding more data from the time series between calculations of the p-value.
 
     Parameters
     ----------
@@ -392,6 +392,7 @@ def get_paired_t_p_timeseries(
 
     np.ndarray
         The times at which the p-values were calculated.
+
     """
     # Check that the data is valid.
     data = check_data(data, one_dim_allowed=False)
@@ -481,14 +482,14 @@ def detect_equilibration_paired_t_test(
     time_units: str = "ns",
     data_y_label: str = r"$\Delta G$ / kcal mol$^{-1}$",
 ) -> _Union[_np.int64, _np.float64]:
-    r"""
-    Detect the equilibration time of a time series by performing a paired
-    t-test between initial and final portions of the time series. This is repeated
-    , discarding more data from the time series between repeats. If the p-value
-    is greater than the threshold, there is no significant evidence that the data is
-    no equilibrated and the timeseries is taken to be equilibrated at this time
-    point. This test may be useful when we care only about systematic bias in the
-    data, and do not care about detecting inter-run differences.
+    r"""Detect the equilibration time of a time series by performing a paired t-test.
+
+    A paired t-test is performed between initial and final portions of the time series.
+    This is repeated, discarding more data from the time series between repeats. If the
+    p-value is greater than the threshold, there is no significant evidence that the data
+    is not equilibrated and the timeseries is taken to be equilibrated at this time point.
+    This test may be useful when we care only about systematic bias in the data, and do
+    not care about detecting inter-run differences.
 
     Parameters
     ----------
@@ -540,6 +541,7 @@ def detect_equilibration_paired_t_test(
     np.float64 | np.int64
         The time (or index, if no times are supplied) at which
         the time series is equilibrated.
+
     """
     # Validate data.
     data = check_data(data, one_dim_allowed=False)
diff --git a/red/ess.py b/red/ess.py
index 35150fb..324153c 100644
--- a/red/ess.py
+++ b/red/ess.py
@@ -16,21 +16,21 @@
 def convert_sse_series_to_ess_series(
     data: _npt.NDArray[_np.float64], sse_series: _npt.NDArray[_np.float64]
 ) -> _npt.NDArray[_np.float64]:
-    """
-    Convert a series of squared standard errors to a series of effective sample sizes.
+    """Convert a series of squared standard errors to a series of effective sample sizes.
 
     Parameters
     ----------
+    data : np.ndarray
+        The time series data with shape (n_chains, n_samples).
+
     sse_series : np.ndarray
         The squared standard error series.
 
-    uncor_vars : np.ndarray
-        The uncorrelated variances.
-
     Returns
     -------
     np.ndarray
         The effective sample size series.
+
     """
     # Validate the data.
     data = check_data(data, one_dim_allowed=True)
@@ -53,10 +53,10 @@ def get_ess_series_init_seq(
     smooth_lag_times: bool = False,
     frac_padding: float = 0.1,
 ) -> _Tuple[_npt.NDArray[_np.float64], _npt.NDArray[_np.float64]]:
-    """
-    Compute a series of effective sample sizes for a time series as data
-    is discarded from the beginning of the time series. The autocorrelation
-    is computed using the sequence estimator specified.
+    """Compute a series of effective sample sizes for a time series.
+
+    Effective sample sizes are computed as data is discarded from the beginning of the
+    time series. The autocorrelation is computed using the sequence estimator specified.
 
     Parameters
     ----------
@@ -95,6 +95,7 @@ def get_ess_series_init_seq(
 
     np.ndarray
         The maximum lag times used.
+
     """
     sse_series, max_lag_times = _get_sse_series_init_seq(
         data,
@@ -116,10 +117,11 @@ def get_ess_series_window(
     window_size_fn: _Optional[_Callable[[int], int]] = lambda x: round(x**0.5),
     window_size: _Optional[int] = None,
 ) -> _Tuple[_npt.NDArray[_np.float64], _npt.NDArray[_np.float64]]:
-    """
-    Compute a series of effective sample sizes for a time series as data
-    is discarded from the beginning of the time series. The squared standard
-    error is computed using the window size and kernel specified.
+    """Compute a series of effective sample sizes for a time series.
+
+    Effective sample sizes are computed as data is discarded from the beginning of the
+    time series. The squared standard error is computed using the window size and kernel
+    specified.
 
     Parameters
     ----------
@@ -144,6 +146,7 @@ def get_ess_series_window(
 
     np.ndarray
         The window sizes used.
+
     """
     sse_series, max_lag_times = _get_sse_series_window(
         data, kernel=kernel, window_size_fn=window_size_fn, window_size=window_size
@@ -155,10 +158,10 @@ def get_ess_series_window(
 
 
 def statistical_inefficiency_inter_variance(data: _npt.NDArray[_np.float64]) -> float:
-    """
-    Compute the statistical inefficiency of a time series by dividing
-    the inter-run variance estimate by the intra-run variance estimate.
-    More than one run is required.
+    """Compute the statistical inefficiency of a time series.
+
+    The statistical inefficiency is computed by dividing the inter-run variance estimate
+    by the intra-run variance estimate. More than one run is required.
 
     Parameters
     ----------
@@ -171,6 +174,7 @@ def statistical_inefficiency_inter_variance(data: _npt.NDArray[_np.float64]) ->
     -------
     float
         The statistical inefficiency.
+
     """
     g = inter_run_variance(data) / intra_run_variance(data)
     # Ensure that the statistical inefficiency is at least 1.
@@ -180,10 +184,11 @@ def statistical_inefficiency_inter_variance(data: _npt.NDArray[_np.float64]) ->
 def statistical_inefficiency_lugsail_variance(
     data: _npt.NDArray[_np.float64], n_pow: float = 1 / 3
 ) -> float:
-    """
-    Compute the statistical inefficiency of a time series by dividing
-    the lugsail replicated batch means variance estimate by the
-    intra-run variance estimate. This is applicable to a single run.
+    """Compute the statistical inefficiency of a time series.
+
+    The statistical inefficiency is computed by dividing the lugsail replicated batch
+    means variance estimate by the intra-run variance estimate. This is applicable to a
+    single run.
 
     Parameters
     ----------
@@ -200,6 +205,7 @@ def statistical_inefficiency_lugsail_variance(
     -------
     float
         The statistical inefficiency.
+
     """
     g = lugsail_variance(data, n_pow=n_pow) / intra_run_variance(data)
     # Ensure that the statistical inefficiency is at least 1.
@@ -207,11 +213,11 @@ def statistical_inefficiency_lugsail_variance(
 
 
 def ess_inter_variance(data: _npt.NDArray[_np.float64]) -> float:
-    """
-    Compute the effective sample size of a time series by dividing
-    the total number of samples by the statistical inefficiency, where
-    the statistical inefficiency is calculated using the ratio of the
-    inter-run and intra-run variance estimates.
+    """Compute the effective sample size of a time series.
+
+    The effective sample size is computed by dividing the total number of samples by the
+    statistical inefficiency, where the statistical inefficiency is calculated using the
+    ratio of the inter-run and intra-run variance estimates.
 
     Parameters
     ----------
@@ -224,6 +230,7 @@ def ess_inter_variance(data: _npt.NDArray[_np.float64]) -> float:
     -------
     float
         The effective sample size.
+
     """
     data = check_data(data, one_dim_allowed=False)
     n_runs: int = data.shape[0]
@@ -233,11 +240,11 @@ def ess_inter_variance(data: _npt.NDArray[_np.float64]) -> float:
 
 
 def ess_lugsail_variance(data: _npt.NDArray[_np.float64], n_pow: float = 1 / 3) -> float:
-    """
-    Compute the effective sample size of a time series by dividing
-    the total number of samples by the statistical inefficiency, where
-    the statistical inefficiency is calculated using the ratio of the
-    lugsail replicated batch means and intra-run variance estimates.
+    """Compute the effective sample size of a time series.
+
+    The effective sample size is computed by dividing the total number of samples by the
+    statistical inefficiency, where the statistical inefficiency is calculated using the
+    ratio of the lugsail replicated batch means and intra-run variance estimates.
 
     Parameters
     ----------
@@ -254,6 +261,7 @@ def ess_lugsail_variance(data: _npt.NDArray[_np.float64], n_pow: float = 1 / 3)
     -------
     float
         The effective sample size.
+
     """
     data = check_data(data, one_dim_allowed=True)
     n_runs: int = data.shape[0]
diff --git a/red/gelman_rubin.py b/red/gelman_rubin.py
index 03f4480..6bae94b 100644
--- a/red/gelman_rubin.py
+++ b/red/gelman_rubin.py
@@ -16,10 +16,10 @@
 
 
 def gelman_rubin(data: _npt.NDArray[_np.float64]) -> float:
-    """
-    Compute the Gelman-Rubin diagnostic according to
-    equation 4 in  Statist. Sci. 36(4): 518-529
-    (November 2021). DOI: 10.1214/20-STS812
+    """Compute the Gelman-Rubin diagnostic.
+
+    Uses equation 4 in Statist. Sci. 36(4): 518-529 (November 2021).
+    DOI: 10.1214/20-STS812.
 
     Parameters
     ----------
@@ -32,6 +32,7 @@ def gelman_rubin(data: _npt.NDArray[_np.float64]) -> float:
     -------
     float
         The Gelman-Rubin diagnostic.
+
     """
     # Check that the data is valid.
     data = _check_data(data, one_dim_allowed=False)
@@ -53,11 +54,10 @@ def gelman_rubin(data: _npt.NDArray[_np.float64]) -> float:
 
 
 def stable_gelman_rubin(data: _npt.NDArray[_np.float64], n_pow: float = 1 / 3) -> float:
-    """
-    Compute the stable Gelman-Rubin diagnostic according to
-    equation 7 in  Statist. Sci. 36(4): 518-529
-    (November 2021). DOI: 10.1214/20-STS812. This is applicable to
-    a single run.
+    """Compute the stable Gelman-Rubin diagnostic.
+
+    Uses equation 7 in Statist. Sci. 36(4): 518-529 (November 2021).
+    DOI: 10.1214/20-STS812. This is applicable to a single run.
     """
     # Validate the data.
     data = _check_data(data, one_dim_allowed=True)
diff --git a/red/plot.py b/red/plot.py
index feb68cf..38d0622 100644
--- a/red/plot.py
+++ b/red/plot.py
@@ -28,8 +28,7 @@ def plot_timeseries(
     time_units: str = "ns",
     y_label: str = r"$\Delta G$ / kcal mol$^{-1}$",
 ) -> None:
-    r"""
-    Plot the (multi-run) time series data.
+    r"""Plot the (multi-run) time series data.
 
     Parameters
     ----------
@@ -55,6 +54,7 @@ def plot_timeseries(
 
     y_label : str, optional
         The y-axis label. The default is "$\Delta G$ / kcal mol$^{-1}$".
+
     """
     # Check that data is valid.
     data = check_data(data, one_dim_allowed=True)
@@ -122,8 +122,7 @@ def plot_p_values(
     time_units: str = "ns",
     threshold_times: _Optional[_npt.NDArray[_np.float64]] = None,
 ) -> None:
-    """
-    Plot the p-values of the paired t-test.
+    """Plot the p-values of the paired t-test.
 
     Parameters
     ----------
@@ -146,6 +145,7 @@ def plot_p_values(
         The times to plot the p-value threshold at. If None, this is
         set to the times at which the data was sampled. Useful when
         using this plot underneath a time series plot.
+
     """
     # Check that p_values is valid.
     if not isinstance(p_values, _np.ndarray) or not isinstance(times, _np.ndarray):
@@ -208,8 +208,7 @@ def plot_sse(
     variance_y_label: str = r"$\frac{1}{\sigma^2(\Delta G)}$ / kcal$^{-2}$ mol$^2$",
     reciprocal: bool = True,
 ) -> _Tuple[_List[_Artist], _List[_Any]]:
-    r"""
-    Plot the squared standard error (SSE) estimate against time.
+    r"""Plot the squared standard error (SSE) estimate against time.
 
     Parameters
     ----------
@@ -245,6 +244,7 @@ def plot_sse(
 
     labels : List[str]
         The labels for the legend.
+
     """
     # Check that sse is valid.
     if not isinstance(sse, _np.ndarray) or not isinstance(times, _np.ndarray):
@@ -318,16 +318,16 @@ def plot_equilibration_paired_t_test(
     time_units: str = "ns",
     data_y_label: str = r"$\Delta G$ / kcal mol$^{-1}$",
 ) -> _Tuple[_Axes, _Axes]:
-    r"""
-    Plot the p-values of the paired t-test against time, underneath the
-    time series data.
+    r"""Plot the p-values of the paired t-test against time.
+
+    The p-values are plotted underneath the time series data.
 
     Parameters
     ----------
     fig : plt.Figure
         The figure to plot on.
 
-    gridspec_obj : plt.GridSpec
+    subplot_spec : plt.GridSpec
         The gridspec to use for the plot.
 
     data : np.ndarray
@@ -361,6 +361,7 @@ def plot_equilibration_paired_t_test(
 
     ax_bottom : Axes
         The axes for the p-values.
+
     """
     with _plt.style.context(PLT_STYLE):
         # We need to split the gridspec into two subplots, one for the time series data (above)
@@ -418,9 +419,9 @@ def plot_equilibration_min_sse(
     variance_y_label: str = r"$\frac{1}{\sigma^2(\Delta G)}$ / kcal$^{-2}$ mol$^2$",
     reciprocal: bool = True,
 ) -> _Tuple[_Axes, _Axes]:
-    r"""
-    Plot the (reciprocal of the) squared standard error (SSE)
-    estimates against time, underneath the time series data.
+    r"""Plot the (reciprocal of the) squared standard error (SSE) estimates against time.
+
+    The estimates are plotted underneath the time series data.
 
     Parameters
     ----------
@@ -475,6 +476,7 @@ def plot_equilibration_min_sse(
 
     ax_bottom : Axes
         The axes for the p-values.
+
     """
     with _plt.style.context(PLT_STYLE):
         data = check_data(data, one_dim_allowed=True)
diff --git a/red/sse.py b/red/sse.py
index 19a1eda..a36c756 100644
--- a/red/sse.py
+++ b/red/sse.py
@@ -19,10 +19,11 @@ def get_sse_series_init_seq(
     smooth_lag_times: bool = False,
     frac_padding: float = 0.1,
 ) -> _Tuple[_npt.NDArray[_np.float64], _npt.NDArray[_np.float64]]:
-    """
-    Compute a series of squared standard errors for a time series as data
-    is discarded from the beginning of the time series. The squared standard
-    error is computed using the sequence estimator specified.
+    """Compute a series of squared standard errors for a time series.
+
+    Squared standard errors are computed as data is discarded from the beginning of the
+    time series. The squared standard error is computed using the sequence estimator
+    specified.
 
     Parameters
     ----------
@@ -61,6 +62,7 @@ def get_sse_series_init_seq(
 
     np.ndarray
         The maximum lag times used.
+
     """
     # Validate the data.
     data = check_data(data, one_dim_allowed=True)
@@ -91,10 +93,11 @@ def get_sse_series_window(
     window_size: _Optional[int] = None,
     frac_padding: float = 0.1,
 ) -> _Tuple[_npt.NDArray[_np.float64], _npt.NDArray[_np.float64]]:
-    """
-    Compute a series of squared standard errors for a time series as data
-    is discarded from the beginning of the time series. The squared standard
-    error is computed using the window size and kernel specified.
+    """Compute a series of squared standard errors for a time series.
+
+    Squared standard errors are computed as data is discarded from the beginning of the
+    time series. The squared standard error is computed using the window size and kernel
+    specified.
 
     Parameters
     ----------
@@ -125,6 +128,7 @@ def get_sse_series_window(
 
     np.ndarray
         The window sizes used.
+
     """
     # Validate the data.
     data = check_data(data, one_dim_allowed=True)
diff --git a/red/variance.py b/red/variance.py
index b1ab069..6116409 100644
--- a/red/variance.py
+++ b/red/variance.py
@@ -1,5 +1,4 @@
-"""
-Functions to calculate the variance of a time series, accounting for autocorrelation.
+"""Functions to calculate the variance of a time series, accounting for autocorrelation.
 
 Methods implemented:
 
@@ -35,8 +34,8 @@
 def _compute_autocovariance_no_fft(
     data: _npt.NDArray[_np.float64], max_lag: int
 ) -> _npt.NDArray[_np.float64]:
-    """
-    Calculate the auto-covariance as a function of lag time for a time series.
+    """Calculate the auto-covariance as a function of lag time for a time series.
+
     Avoids using statsmodel's acovf function as using numpy's dot function and jit
     gives a substantial speedup.
 
@@ -54,6 +53,7 @@ def _compute_autocovariance_no_fft(
     -------
     numpy.ndarray
         The auto-correlation function of the time series.
+
     """
     # Don't use statsmodel's acovf as we can get a substantial speedup by using
     # numpy's dot function and jit.
@@ -74,10 +74,11 @@ def _compute_autocovariance_no_fft(
 def _compute_autocovariance_fft(
     data: _npt.NDArray[_np.float64], max_lag: int
 ) -> _npt.NDArray[_np.float64]:
-    """
-    Calculate the autocovariance using the FFT method, as implemented in statsmodels.
+    """Calculate the autocovariance using the FFT method, as implemented in statsmodels.
+
     Note that we can speed this up for large arrays by rewriting to directly use numpy's fft
     function and using jit with rocket-fft https://github.com/styfenschaer/rocket-fft.
+
     Parameters
     ----------
     data : numpy.ndarray
@@ -92,6 +93,7 @@ def _compute_autocovariance_fft(
     -------
     numpy.ndarray
     The auto-correlation function of the time series.
+
     """
     autocov_fn = _acovf(data, adjusted=False, nlag=max_lag, fft=True, demean=False)
     autocov_fn = autocov_fn.astype(_np.float64, copy=False)
@@ -104,8 +106,7 @@ def _get_autocovariance(
     mean: _Union[None, float] = None,
     fft: bool = False,
 ) -> _npt.NDArray[_np.float64]:
-    """
-    Calculate the auto-covariance as a function of lag time for a time series.
+    """Calculate the auto-covariance as a function of lag time for a time series.
 
     Parameters
     ----------
@@ -130,6 +131,7 @@ def _get_autocovariance(
     -------
     numpy.ndarray
         The auto-correlation function of the time series.
+
     """
     # Copy the data so we don't modify the original.
     data = data.copy()
@@ -157,8 +159,7 @@ def _get_autocovariance(
 def _get_gamma_cap(
     autocov_series: _npt.NDArray[_np.float64],
 ) -> _npt.NDArray[_np.float64]:
-    """
-    Compute the capitial gamma function from the auto-covariance function.
+    """Compute the capitial gamma function from the auto-covariance function.
 
     Parameters
     ----------
@@ -169,6 +170,7 @@ def _get_gamma_cap(
     -------
     numpy.ndarray
         The capital gamma function of the time series.
+
     """
     # Get the length of the time series.
     n_samples = autocov_series.shape[0]
@@ -195,8 +197,8 @@ def _get_initial_positive_sequence(
     gamma_cap: _npt.NDArray[_np.float64],
     min_max_lag_time: int = 3,
 ) -> _npt.NDArray[_np.float64]:
-    """ "
-    Get the initial positive sequence from the capital gamma function of a time series.
+    """Get the initial positive sequence from the capital gamma function of a time series.
+
     See Geyer, 1992: https://www.jstor.org/stable/2246094.
 
     Parameters
@@ -212,6 +214,7 @@ def _get_initial_positive_sequence(
     -------
     numpy.ndarray
         The initial positive sequence.
+
     """
     # Make a copy of gamma_cap so we don't modify the original.
     # gamma_cap = _deepcopy(gamma_cap)
@@ -231,8 +234,8 @@ def _get_initial_monotone_sequence(
     gamma_cap: _npt.NDArray[_np.float64],
     min_max_lag_time: int = 3,
 ) -> _npt.NDArray[_np.float64]:
-    """
-    Get the initial monotone sequence from the capital gamma function of a time series.
+    """Get the initial monotone sequence from the capital gamma function of a time series.
+
     See Geyer, 1992: https://www.jstor.org/stable/2246094.
 
     Parameters
@@ -247,6 +250,7 @@ def _get_initial_monotone_sequence(
     -------
     numpy.ndarray
         The initial monotone sequence.
+
     """
     # Make a copy of gamma_cap so we don't modify the original.
     gamma_cap = gamma_cap.copy()
@@ -267,8 +271,8 @@ def _get_initial_convex_sequence(
     gamma_cap: _npt.NDArray[_np.float64],
     min_max_lag_time: int = 3,
 ) -> _npt.NDArray[_np.float64]:
-    """
-    Get the initial convex sequence from the capital gamma function of a time series.
+    """Get the initial convex sequence from the capital gamma function of a time series.
+
     See Geyer, 1992: https://www.jstor.org/stable/2246094.
 
     Parameters
@@ -290,6 +294,7 @@ def _get_initial_convex_sequence(
     MIT License.
     YEAR: 2005, 2009, 2010, 2012
     COPYRIGHT HOLDER: Charles J. Geyer and Leif T. Johnson
+
     """
     # Make a copy of gamma_cap so we don't modify the original.
     gamma_con = gamma_cap.copy()
@@ -346,8 +351,7 @@ def _get_autocovariance_window(
     kernel: _Callable[[int], _npt.NDArray[_np.float64]] = _np.bartlett,  # type: ignore
     window_size: int = 10,
 ) -> _npt.NDArray[_np.float64]:
-    """
-    Calculate the autocovariance of a time series using window estimators.
+    """Calculate the autocovariance of a time series using window estimators.
 
     Parameters
     ----------
@@ -366,6 +370,7 @@ def _get_autocovariance_window(
     -------
     np.ndarray
         The autocovariance of the time series as a function of lag time.
+
     """
     n_runs, n_samples = data.shape
     if n_samples < window_size:
@@ -396,10 +401,11 @@ def _get_autocovariance_window(
 def _smoothen_max_lag_times(
     max_lag_times: _npt.NDArray[_np.float64],
 ) -> _npt.NDArray[_np.int64]:
-    """
-    Smoothen a list of maximum lag times by a) converting them to a monotinically
-    decreasing sequence and b) linearly interpolating between points where the sequence
-    changes. This may be useful when the max lag times are noisy.
+    """Smoothen a list of maximum lag times.
+
+    This is done by a) converting them to a monotonically decreasing sequence and b)
+    linearly interpolating between points where the sequence changes. This may be useful
+    when the max lag times are noisy.
 
     Parameters
     ----------
@@ -410,6 +416,7 @@ def _smoothen_max_lag_times(
     -------
     numpy.ndarray
         The smoothened maximum lag times.
+
     """
     # Get a monotinically decreasing sequence.
     max_lag_times_monotonic = _get_initial_monotone_sequence(max_lag_times, min_max_lag_time=0)
@@ -445,8 +452,8 @@ def get_variance_initial_sequence(
     max_max_lag_time: _Optional[int] = None,
     autocov: _Optional[_npt.NDArray[_np.float64]] = None,
 ) -> _Tuple[float, int, _npt.NDArray[_np.float64]]:
-    """
-    Calculate the variance of a time series using initial sequence methods.
+    """Calculate the variance of a time series using initial sequence methods.
+
     See Geyer, 1992: https://www.jstor.org/stable/2246094.
 
     Parameters
@@ -482,6 +489,7 @@ def get_variance_initial_sequence(
 
     numpy.ndarray
         The auto-covariance function of the time series.
+
     """
     # Validate the data.
     data = _check_data(data, one_dim_allowed=True)
@@ -611,10 +619,11 @@ def get_variance_series_initial_sequence(
     smooth_lag_times: bool = False,
     frac_padding: float = 0.1,
 ) -> _Tuple[_npt.NDArray[_np.float64], _npt.NDArray[_np.float64]]:
-    """
-    Repeatedly calculate the variance of a time series while discarding increasing
-    numbers of samples from the start of the time series. The variance is calculated
-    using initial sequence methods. See Geyer, 1992: https://www.jstor.org/stable/2246094.
+    """Repeatedly calculate the variance of a time series while discarding samples.
+
+    Increasing numbers of samples are discarded from the start of the time series. The
+    variance is calculated using initial sequence methods. See Geyer, 1992:
+    https://www.jstor.org/stable/2246094.
 
     Parameters
     ----------
@@ -653,6 +662,7 @@ def get_variance_series_initial_sequence(
 
     numpy.ndarray
         The maximum lag time used when calculating the auto-correlated variance.
+
     """
     # Check that the data is valid.
     data = _check_data(data, one_dim_allowed=True)
@@ -740,8 +750,7 @@ def get_variance_window(
     kernel: _Callable[[int], _npt.NDArray[_np.float64]] = _np.bartlett,  # type: ignore
     window_size: int = 10,
 ) -> float:
-    """
-    Calculate the variance of a time series using window estimators.
+    """Calculate the variance of a time series using window estimators.
 
     Parameters
     ----------
@@ -759,6 +768,7 @@ def get_variance_window(
     -------
     float
         The estimated variance of the time series.
+
     """
     # Check that the data is valid.
     data = _check_data(data, one_dim_allowed=True)
@@ -802,10 +812,10 @@ def get_variance_series_window(
     window_size: _Optional[int] = None,
     frac_padding: float = 0.1,
 ) -> _Tuple[_npt.NDArray[_np.float64], _npt.NDArray[_np.float64]]:
-    """
-    Repeatedly calculate the variance of a time series while discarding increasing
-    numbers of samples from the start of the time series. The variance is calculated
-    using window estimators.
+    """Repeatedly calculate the variance of a time series while discarding samples.
+
+    Increasing numbers of samples are discarded from the start of the time series. The
+    variance is calculated using window estimators.
 
     Parameters
     ----------
@@ -836,6 +846,7 @@ def get_variance_series_window(
 
     numpy.ndarray
         The window size used at each index.
+
     """
     # Check that the data is valid.
     data = _check_data(data, one_dim_allowed=True)
@@ -888,8 +899,8 @@ def get_variance_series_window(
 
 
 def replicated_batch_means_variance(data: _npt.NDArray[_np.float64], batch_size: int) -> float:
-    """
-    Estimate the variance of a time series using the replicated batch means method.
+    """Estimate the variance of a time series using the replicated batch means method.
+
     See section 3.1 in Statist. Sci. 36(4): 518-529 (November 2021).
     DOI: 10.1214/20-STS812 .
 
@@ -905,6 +916,7 @@ def replicated_batch_means_variance(data: _npt.NDArray[_np.float64], batch_size:
     -------
     float
         The estimated variance.
+
     """
     data = _check_data(data, one_dim_allowed=True)
 
@@ -935,8 +947,8 @@ def replicated_batch_means_variance(data: _npt.NDArray[_np.float64], batch_size:
 
 
 def lugsail_variance(data: _npt.NDArray[_np.float64], n_pow: float = 1 / 3) -> float:
-    """
-    Estimate the variance of a time series using the lugsail method.
+    """Estimate the variance of a time series using the lugsail method.
+
     See section 3.2 in Statist. Sci. 36(4): 518-529 (November 2021).
     DOI: 10.1214/20-STS812 .
 
@@ -953,6 +965,7 @@ def lugsail_variance(data: _npt.NDArray[_np.float64], n_pow: float = 1 / 3) -> f
     -------
     float
         The estimated variance.
+
     """
     # Check that the data is valid.
     data = _check_data(data, one_dim_allowed=True)
@@ -983,9 +996,7 @@ def lugsail_variance(data: _npt.NDArray[_np.float64], n_pow: float = 1 / 3) -> f
 
 
 def inter_run_variance(data: _npt.NDArray[_np.float64]) -> float:
-    """
-    Compute the variance based on the inter-run differences
-    between means.
+    """Compute the variance based on the inter-run differences between means.
 
     Parameters
     ----------
@@ -996,6 +1007,7 @@ def inter_run_variance(data: _npt.NDArray[_np.float64]) -> float:
     -------
     float
         The estimated variance.
+
     """
     # Check that the data is valid.
     data = _check_data(data, one_dim_allowed=False)
@@ -1011,8 +1023,7 @@ def inter_run_variance(data: _npt.NDArray[_np.float64]) -> float:
 
 
 def intra_run_variance(data: _npt.NDArray[_np.float64]) -> float:
-    """
-    Compute the average intra-run variance estimate.
+    """Compute the average intra-run variance estimate.
 
     Parameters
     ----------
@@ -1023,6 +1034,7 @@ def intra_run_variance(data: _npt.NDArray[_np.float64]) -> float:
     -------
     float
         The mean intra-run variance estimate.
+
     """
     # Check that the data is valid.
     data = _check_data(data, one_dim_allowed=True)

From 7888c80cdfc5a9aff6047356efbfd167d3d448b9 Mon Sep 17 00:00:00 2001
From: Finlay Clark <finlay.clark@newcastle.ac.uk>
Date: Thu, 6 Aug 2026 12:27:57 +0100
Subject: [PATCH 2/3] Fix typing messed up in merge

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

diff --git a/red/ess.py b/red/ess.py
index ad4185e..b9dfe94 100644
--- a/red/ess.py
+++ b/red/ess.py
@@ -111,9 +111,9 @@ def get_ess_series_init_seq(
 
 def get_ess_series_window(
     data: _npt.NDArray[_np.float64],
-    kernel: _Callable[[int], _npt.NDArray[_np.float64]] = _np.bartlett,  # type: ignore
-    window_size_fn: _Optional[_Callable[[int], int]] = lambda x: round(x**0.5),
-    window_size: _Optional[int] = None,
+    kernel: _Callable[[int], _npt.NDArray[_np.float64]] = _np.bartlett,
+    window_size_fn: _Callable[[int], int] | None = lambda x: round(x**0.5),
+    window_size: int | None = None,
 ) -> tuple[_npt.NDArray[_np.float64], _npt.NDArray[_np.float64]]:
     """Compute a series of effective sample sizes for a time series.
 

From 5534ef1a13bd816d13005511ff7077837e5da2d7 Mon Sep 17 00:00:00 2001
From: Finlay Clark <finlay.clark@newcastle.ac.uk>
Date: Thu, 6 Aug 2026 12:39:21 +0100
Subject: [PATCH 3/3] fix: ruff issues

---
 README.md         | 15 +++++----------
 docs/changelog.md |  1 +
 docs/examples.md  | 28 +++++++++++++++++++++-------
 red/_version.py   |  2 +-
 red/plot.py       |  7 +++----
 5 files changed, 31 insertions(+), 22 deletions(-)

diff --git a/README.md b/README.md
index 0758dd0..dd85a79 100644
--- a/README.md
+++ b/README.md
@@ -76,21 +76,17 @@ my_timeseries = ...
 # of the first sample after equilibration, g is the statistical
 # inefficiency of the equilibrated sample, and ess is the effective sample
 # size of the equilibrated sample.
-idx, g, ess = red.detect_equilibration_window(my_timeseries,
-                                              method="min_sse",
-                                              plot=True)
+idx, g, ess = red.detect_equilibration_window(my_timeseries, method="min_sse", plot=True)
 
 # Alternatively, use Geyer's initial convex sequence method to account
 # for autocorrelation.
-idx, g, ess = red.detect_equilibration_init_seq(my_timeseries,
-                                                method="min_sse",
-                                                plot=True)
+idx, g, ess = red.detect_equilibration_init_seq(my_timeseries, method="min_sse", plot=True)
 
 # We can also determine equilibration in the same way as in
 # pymbar.timeseries.detect_equilibration(my_timeseries, fast=False)
-idx, g, ess = red.detect_equilibration_init_seq(my_timeseries,
-                                                method="max_ess",
-                                                sequence_estimator="positive")
+idx, g, ess = red.detect_equilibration_init_seq(
+    my_timeseries, method="max_ess", sequence_estimator="positive"
+)
 ```
 
 #### Uncertainty Quantification
@@ -99,7 +95,6 @@ idx, g, ess = red.detect_equilibration_init_seq(my_timeseries,
 # Estimate the 95 % confidence interval, accounting for autocorrelation using Geyer's initial
 # convex sequence method.
 ci_95 = red.get_conf_int_init_seq(my_timeseries, alpha_two_tailed=0.05)
-
 ```
 
 For more examples, see the [documentation](https://fjclark.github.io/red/latest/examples/).
diff --git a/docs/changelog.md b/docs/changelog.md
index b3488fa..cc10d7c 100644
--- a/docs/changelog.md
+++ b/docs/changelog.md
@@ -2,6 +2,7 @@
 
 ## 0.2.0 - 2026-08-06
 
+- Make ruff stricter about docstrings and fix them up.
 - Switch to from conda, make, and mypy to uv, just, and ty. Modernise type hinting and switch to Python > 3.11.
 - Drop polyfill.io due to security issues.
 
diff --git a/docs/examples.md b/docs/examples.md
index 1a9577a..5470f39 100644
--- a/docs/examples.md
+++ b/docs/examples.md
@@ -26,7 +26,9 @@ For all examples, `my_timeseries` should be a numpy array with shape `(n_samples
 To use any of Geyer's initial sequence methods ([Geyer, 1992](https://www.jstor.org/stable/2246094)), you can specify the "sequence_estimator" to be "initial_positive" (the least strict), "initial_monotone", or "initial_convex" (the strictest):
 
 ```python
-idx, g, ess = red.detect_equilibration_init_seq(my_timeseries, sequence_estimator="initial_convex", plot=True)
+idx, g, ess = red.detect_equilibration_init_seq(
+    my_timeseries, sequence_estimator="initial_convex", plot=True
+)
 my_truncated_timeseries = my_timeseries[idx:]
 ```
 To use Chodera's method of simply truncating the autocovariance series at the first negative value ([Chodera, 2016](https://pubs.acs.org/doi/full/10.1021/acs.jctc.5b00784)), you can specify the "sequence estimator" to be "positive".
@@ -36,7 +38,9 @@ To use Chodera's method of simply truncating the autocovariance series at the fi
 When using window methods, you can either specify a fixed window size, or a window size function which computes the window size as a function of the number of data points (which decreases as the truncation point increases). These are specified via `window_size` and `window_size_fn`, respectively (one must be specified and the other must be `None`). The default window size function is `lambda x: round(x**0.5)` - explicitly:
 
 ```python
-idx, g, ess = red.detect_equilibration_window(my_timeseries, window_size=None, window_size_fn=lambda x: round(x**0.5), plot=True)
+idx, g, ess = red.detect_equilibration_window(
+    my_timeseries, window_size=None, window_size_fn=lambda x: round(x**0.5), plot=True
+)
 # This is equivalent to:
 idx, g, ess = red.detect_equilibration_window(my_timeseries, plot=True)
 ```
@@ -44,7 +48,9 @@ idx, g, ess = red.detect_equilibration_window(my_timeseries, plot=True)
 To use a window size of 10:
 
 ```python
-idx, g, ess = red.detect_equilibration_window(my_timeseries, window_size=10, window_size_fn = None, plot=True)
+idx, g, ess = red.detect_equilibration_window(
+    my_timeseries, window_size=10, window_size_fn=None, plot=True
+)
 ```
 
 You can also play with the kernel function used in the window method by specifying the `kernel` argument. You should supply the function directly - the default is `np.bartlett`.
@@ -54,7 +60,9 @@ You can also play with the kernel function used in the window method by specifyi
 To use White's original Marginal Standard Error Rule ([White, 1997](https://journals.sagepub.com/doi/abs/10.1177/003754979706900601)), you can use the window method with a window size of 1:
 
 ```python
-idx, g, ess = red.detect_equilibration_window(my_timeseries, window_size=1, window_size_fn=None, plot=True)
+idx, g, ess = red.detect_equilibration_window(
+    my_timeseries, window_size=1, window_size_fn=None, plot=True
+)
 ```
 
 ### Maximum Effective Sample Size and Chodera's Method
@@ -62,7 +70,9 @@ idx, g, ess = red.detect_equilibration_window(my_timeseries, window_size=1, wind
 To select the truncation point according to the maximum effective sample size (instead of the minimum squared standard error), you can specify the `method` argument to be "max_ess". To use Chodera's method ([Chodera, 2016](https://pubs.acs.org/doi/full/10.1021/acs.jctc.5b00784)) as implemented in `pymbar.timeseries`, you can specify the `sequence_estimator` to be "positive":
 
 ```python
-idx, g, ess = red.detect_equilibration_init_seq(my_timeseries, method="max_ess", sequence_estimator="positive", plot=True)
+idx, g, ess = red.detect_equilibration_init_seq(
+    my_timeseries, method="max_ess", sequence_estimator="positive", plot=True
+)
 # Equivalent to pymbar.timeseries.detect_equilibration(my_timeseries, fast=False)
 ```
 
@@ -71,7 +81,9 @@ idx, g, ess = red.detect_equilibration_init_seq(my_timeseries, method="max_ess",
 To save a plot showing the (block-averaged) time series and variance of the mean/ effective sample size against truncation time, simply specify `plot=True` and, optionally, specify a name for the plot with `plot_name`. This works for either of the equilbration detection functions.
 
 ```python
-idx, g, ess = red.detect_equilibration_window(my_timeseries, plot=True, plot_name="my_equilibration_plot.png")
+idx, g, ess = red.detect_equilibration_window(
+    my_timeseries, plot=True, plot_name="my_equilibration_plot.png"
+)
 ```
 
 ## Estimating Uncertainty
@@ -79,7 +91,9 @@ idx, g, ess = red.detect_equilibration_window(my_timeseries, plot=True, plot_nam
 To calculate uncertainty, we recommend using Geyer's initial convex sequence method ([Geyer, 1992](https://www.jstor.org/stable/2246094)), which is the default for [`get_conf_int_init_seq`][red.confidence_intervals.get_conf_int_init_seq]. For example, to estimate a 95 % confidence interval:
 
 ```python
-ci_95 = red.get_conf_int_init_seq(my_timeseries, sequence_estimator="initial_convex", alpha_two_tailed=0.05)
+ci_95 = red.get_conf_int_init_seq(
+    my_timeseries, sequence_estimator="initial_convex", alpha_two_tailed=0.05
+)
 ```
 
 This function has a similar interface to [`detect_equilibration_init_seq`][red.equilibration.detect_equilibration_init_seq], and you can specify the "sequence_estimator" in the same way. Note that this assumes we have a reasonable effective sample size, and hence that the means are approximately normally distributed by the central limit theorem.
diff --git a/red/_version.py b/red/_version.py
index 3c07d0f..207b274 100644
--- a/red/_version.py
+++ b/red/_version.py
@@ -1 +1 @@
-__version__ = "0.1.4+6.gf1edc4b.dirty"
+__version__ = "0.1.4+5.g7888c80.dirty"
diff --git a/red/plot.py b/red/plot.py
index dc4b250..5f34097 100644
--- a/red/plot.py
+++ b/red/plot.py
@@ -204,8 +204,7 @@ def plot_sse(
     variance_y_label: str = r"$\frac{1}{\sigma^2(\Delta G)}$ / kcal$^{-2}$ mol$^2$",
     reciprocal: bool = True,
 ) -> tuple[list[_Artist], list[_Any]]:
-    """
-    Plot the squared standard error (SSE) estimate against time.
+    r"""Plot the squared standard error (SSE) estimate against time.
 
     Parameters
     ----------
@@ -229,7 +228,7 @@ def plot_sse(
 
     variance_y_label : str, optional
         The y-axis label for the variance. The default is
-        "$\frac{1}{\sigma^2(\Delta G)}$ / kcal$^{-2}$ mol$^2$".
+        "$\frac{1}{\\sigma^2(\\Delta G)}$ / kcal$^{-2}$ mol$^2$".
 
     reciprocal : bool, optional, default=True
         Whether to plot the reciprocal of the SSE.
@@ -316,7 +315,7 @@ def plot_equilibration_paired_t_test(
     time_units: str = "ns",
     data_y_label: str = r"$\Delta G$ / kcal mol$^{-1}$",
 ) -> tuple[_Axes, _Axes]:
-    """Plot the p-values of the paired t-test against time.
+    r"""Plot the p-values of the paired t-test against time.
 
     The p-values are plotted underneath the time series data.