alnn_rewrite/analysis/plots/bootstrap_ci.py

"""Confidence intervals for full-ensemble accuracy.

Patient-clustered percentile bootstrap: whole PATIENTS are resampled with
replacement, not individual scans, because scans from one patient are correlated
and a per-scan bootstrap would understate the interval. Holds the trained models
fixed, so it captures test-set sampling noise only -- an honest lower bound on
total uncertainty (it excludes training/initialisation variability).
"""

from typing import Any, Dict

import numpy as np

from analysis.plotting import panel_grid, save_figure, save_json, style_axes
from analysis.sources import FAMILY_LABELS, ordered_families
from analysis.context import AnalysisContext, FAMILY_COLORS
from analysis.registry import register

BOOTSTRAP_ITERS = 2000
CI_LEVEL = 0.95


def _bootstrap_accuracy_ci(
    correct: np.ndarray,
    ptids: np.ndarray,
    iters: int,
    ci_level: float,
    rng: np.random.Generator,
) -> Dict[str, Any]:
    """Patient-clustered percentile bootstrap CI for a mean-accuracy statistic.

    Reduces to an ordinary bootstrap when every patient has exactly one scan.
    """
    point = float(np.mean(correct)) if correct.size else 0.0
    _, inverse = np.unique(ptids, return_inverse=True)
    n_patients = int(inverse.max()) + 1 if inverse.size else 0
    groups = [np.where(inverse == p)[0] for p in range(n_patients)]

    boot = np.empty(iters, dtype=float)
    for b in range(iters):
        chosen = rng.integers(0, n_patients, size=n_patients)
        boot[b] = correct[np.concatenate([groups[c] for c in chosen])].mean()

    tail = (1.0 - ci_level) / 2.0
    return {
        "point": point,
        "ci_low": float(np.quantile(boot, tail)),
        "ci_high": float(np.quantile(boot, 1.0 - tail)),
        "n_samples": int(correct.size),
        "n_patients": n_patients,
    }


@register(
    "bootstrap_ensemble_ci",
    title="Bootstrap confidence intervals for ensemble accuracy",
)
def bootstrap_ensemble_ci(ctx: AnalysisContext) -> None:
    base_seed = ctx.base_seed()

    results: Dict[str, Dict[str, Any]] = {}
    for index, family in enumerate(ctx.families):
        selected = ctx.clean_source(family)
        if selected is None:
            continue
        stem, ds = selected
        probs = ctx.probs(family, ds)
        if probs.shape[0] < 1:
            continue
        # Full-ensemble prediction: argmax of the mean over ALL members.
        correct = (
            probs.mean(axis=0).argmax(axis=1) == ctx.true_index(ds)
        ).astype(float)
        rng = np.random.default_rng([base_seed, 1000 + index])
        ci = _bootstrap_accuracy_ci(
            correct, ctx.ptids(ds), BOOTSTRAP_ITERS, CI_LEVEL, rng
        )
        ci["source"] = stem
        results[family] = ci
        ctx.log.info(
            f"bootstrap_ensemble_ci: {family} full-ensemble acc = "
            f"{ci['point'] * 100:.1f}% (95% CI {ci['ci_low'] * 100:.1f}-"
            f"{ci['ci_high'] * 100:.1f}%, n={ci['n_samples']} scans / "
            f"{ci['n_patients']} patients)."
        )

    if not results:
        ctx.log.error("bootstrap_ensemble_ci: nothing to compute.")
        return

    ordered = ordered_families(results.keys())
    fig, axes = panel_grid(
        1, n_cols=1, panel_size=(7.0, 1.4 + 0.8 * len(ordered))
    )
    ax = axes[0]
    for row, family in enumerate(ordered):
        ci = results[family]
        color = FAMILY_COLORS.get(family, "#666666")
        ax.plot(
            [ci["ci_low"], ci["ci_high"]],
            [row, row],
            color=color,
            linewidth=2.5,
            solid_capstyle="round",
        )
        ax.plot([ci["point"]], [row], "o", color=color, markersize=8)
        ax.annotate(
            f"{ci['point'] * 100:.1f}%  "
            f"[{ci['ci_low'] * 100:.1f}, {ci['ci_high'] * 100:.1f}]",
            (ci["point"], row),
            xytext=(0, 10),
            textcoords="offset points",
            ha="center",
            va="bottom",
            fontsize=9,
            color="#222222",
        )

    ax.set_yticks(range(len(ordered)))
    ax.set_yticklabels([FAMILY_LABELS.get(k, k) for k in ordered])
    ax.invert_yaxis()  # first family on top
    ax.set_xlabel("Full-ensemble accuracy (95% patient-clustered bootstrap CI)")
    ax.set_title("Ensemble accuracy with bootstrap confidence intervals")
    ax.margins(x=0.18, y=0.35)
    style_axes(ax, grid_axis="x")
    ax.spines["left"].set_visible(False)
    ax.tick_params(left=False)

    png = save_figure(fig, ctx.out_dir, "bootstrap_ci")
    save_json(
        {"iters": BOOTSTRAP_ITERS, "ci_level": CI_LEVEL, "families": results},
        ctx.out_dir,
        "bootstrap_ci",
    )
    ctx.log.info(
        f"bootstrap_ensemble_ci: wrote {png.name} ({len(results)} family/families)."
    )