Breakdown grid v2: paired six-arm estimator, honest bootstrap errors, degeneracy guards

scripts/python/plot_breakdown_grid.py

@@ -2,42 +2,70 @@
 # requires-python = ">=3.10"
 # dependencies = ["numpy", "matplotlib", "seaborn"]
 # ///
-"""Column figure for the resolution x noise breakdown grid.
-
-Reads the JSON emitted by ``run_breakdown_grid.py`` and draws three vertically
-stacked panels sharing the x-axis -- multiplicative bias ``m``, metacal response
-``R11``, and additive bias ``|c1|`` vs the resolution ratio ``r`` -- with the
-three S/N bands as hue. The ``m`` panel is linear with a widened y-range so the
-moderate/low-S/N climb stays on-plot while the +/-5e-3 tolerance band is still
-drawn. Seaborn is not in the ShapePipe sif, so this runs on the host:
+"""Column figure for the resolution x noise breakdown grid (v2).
+
+Reads the JSON emitted by ``run_breakdown_grid.py`` (v2 schema) and draws three
+vertically stacked panels sharing the x-axis -- multiplicative bias ``m1`` and
+``m2``, and additive bias ``|c1|`` vs the resolution ratio ``r`` -- with the
+three S/N bands as hue. **Every point carries its bootstrap error bar**; that is
+the whole point of v2. Cells whose response was consistent with zero (the
+degeneracy guard fired, ``m`` null) are drawn as open markers at 0 with a hatch
+so a reader never mistakes "we refused to divide" for "m = 0". The ``m`` panels
+keep the ``+/-5e-3`` tolerance band drawn as a sliver at 0. Seaborn is not in the
+ShapePipe sif, so this runs on the host:
 
     uv run scripts/python/plot_breakdown_grid.py \\
-        --input <run>/results/baseline/breakdown_grid/breakdown_grid.json \\
+        --input <run>/results/baseline/breakdown_grid/breakdown_grid_v2.json \\
         --output <figures>/breakdown_grid.png
 """
 import argparse
 import json
 
 import numpy as np
 import matplotlib.pyplot as plt
+from matplotlib.lines import Line2D
 import seaborn as sns
 
 BANDS = [("high", r"S/N $\approx$ 1000"), ("moderate", r"S/N $\approx$ 50"),
          ("low", r"S/N $\approx$ 15")]
 
 
-def series(cells, band):
+def series(cells, band, key):
+    """(r, value, err, degenerate-mask) for one band and one estimator key.
+
+    ``value``/``err`` are NaN where the cell is degenerate (m null); the mask
+    flags those so the plotter can render them as open markers at 0."""
     rows = sorted((c for c in cells if c["s2n_band"] == band),
                   key=lambda c: c["res_ratio"])
-    return (np.array([c["res_ratio"] for c in rows]),
-            np.array([c["m1"] for c in rows]),
-            np.array([c["R11"] for c in rows]),
-            np.array([abs(c["c1"]) for c in rows]))
+    r = np.array([c["res_ratio"] for c in rows])
+    deg_key = "c_degenerate" if key.startswith("c") else f"{key}_degenerate"
+    deg = np.array([bool(c.get(deg_key, False)) or c[key] is None for c in rows],
+                   dtype=bool)
+    val = np.array([np.nan if c[key] is None else c[key] for c in rows], dtype=float)
+    err = np.array([np.nan if c.get(f"{key}_err") is None else c[f"{key}_err"]
+                    for c in rows], dtype=float)
+    return r, val, err, deg
+
+
+def draw_panel(ax, cells, key, pal, transform=lambda v: v):
+    """Plot one estimator across all bands with error bars; degenerate cells as
+    open hatched markers at 0."""
+    for band, _ in BANDS:
+        r, val, err, deg = series(cells, band, key)
+        good = ~deg
+        ax.errorbar(r[good], transform(val[good]),
+                    yerr=None if err is None else err[good],
+                    color=pal[band], marker="o", ms=7, lw=2, capsize=3,
+                    elinewidth=1.4, label=dict(BANDS)[band])
+        if deg.any():  # open markers at 0 mark "response ~ 0, refused to divide"
+            ax.scatter(r[deg], np.zeros(deg.sum()), s=90, facecolors="none",
+                       edgecolors=pal[band], linewidths=1.8, hatch="////",
+                       zorder=5)
 
 
 def main():
     p = argparse.ArgumentParser(description=__doc__.split("\n\n")[0])
-    p.add_argument("--input", required=True, help="breakdown_grid.json path.")
+    p.add_argument("--input", required=True, help="breakdown_grid JSON path (v2).")
     p.add_argument("--output", required=True, help="output PNG path.")
     a = p.parse_args()
 
@@ -48,38 +76,35 @@ def main():
     sns.set_theme(style="ticks", context="talk")
     pal = dict(zip([b for b, _ in BANDS], sns.color_palette("husl", len(BANDS))))
 
-    fig, (ax_m, ax_R, ax_c) = plt.subplots(3, 1, figsize=(7.6, 12.5), sharex=True)
-
-    for band, label in BANDS:
-        r, m, R11, ac1 = series(cells, band)
-        kw = dict(color=pal[band], marker="o", ms=7, lw=2, label=label)
-        ax_m.plot(r, m, **kw)
-        ax_R.plot(r, R11, **kw)
-        ax_c.plot(r, ac1, **kw)
-
-    # m panel: linear, y-range widened to hold the moderate/low-S/N climb while
-    # the +/-tol band stays drawn as a sliver at 0.
-    ax_m.axhspan(-m_tol, m_tol, color="0.75", alpha=0.4, zorder=0,
-                 label=r"$|m| < 5\times10^{-3}$ tol.")
-    ax_m.axhline(0, color="0.5", lw=0.8, zorder=0)
-    ax_m.set_ylim(-0.05, 1.5)
-    ax_m.set_ylabel("multiplicative bias  $m$")
-    ax_m.set_title("Multiplicative bias")
-
-    ax_R.axhline(1.0, color="0.5", lw=0.8, ls="--", zorder=0)
-    ax_R.set_ylim(0, 1.15)
-    ax_R.set_ylabel("metacal response  $R_{11}$")
-    ax_R.set_title("Response")
-
-    ax_c.set_yscale("log")
-    ax_c.set_ylabel("additive bias  $|c_1|$")
+    fig, (ax_m1, ax_m2, ax_c) = plt.subplots(3, 1, figsize=(7.6, 12.5), sharex=True)
+
+    for ax, key, title, ylab in (
+        (ax_m1, "m1", "Multiplicative bias $m_1$", r"multiplicative bias  $m_1$"),
+        (ax_m2, "m2", "Multiplicative bias $m_2$", r"multiplicative bias  $m_2$"),
+    ):
+        draw_panel(ax, cells, key, pal)
+        ax.axhspan(-m_tol, m_tol, color="0.75", alpha=0.4, zorder=0,
+                   label=r"$|m| < 5\times10^{-3}$ tol.")
+        ax.axhline(0, color="0.5", lw=0.8, zorder=0)
+        ax.set_ylabel(ylab)
+        ax.set_title(title)
+
+    # |c1| on a symlog axis so the tolerance-scale values and any blow-up both
+    # stay legible while error bars that cross zero still render.
+    draw_panel(ax_c, cells, "c1", pal, transform=np.abs)
+    ax_c.set_yscale("symlog", linthresh=1e-4)
+    ax_c.set_ylabel(r"additive bias  $|c_1|$")
     ax_c.set_title("Additive bias")
     ax_c.set_xlabel("resolution ratio  $r = \\mathrm{gal\\ hlr}\\,/\\,\\mathrm{psf\\ fwhm}$")
 
-    for ax in (ax_m, ax_R, ax_c):
+    for ax in (ax_m1, ax_m2, ax_c):
         ax.margins(x=0.04)
 
-    ax_m.legend(frameon=False, fontsize=12, loc="upper right")
+    handles = [Line2D([0], [0], color=pal[b], marker="o", lw=2, label=lab)
+               for b, lab in BANDS]
+    handles.append(Line2D([0], [0], marker="o", ls="none", mfc="none",
+                          mec="0.3", mew=1.8, label="degenerate ($R\\approx0$)"))
+    ax_m1.legend(handles=handles, frameon=False, fontsize=12, loc="upper left")
     sns.despine(fig)
     fig.tight_layout()
     fig.savefig(a.output, dpi=140, bbox_inches="tight")

scripts/python/plot_s4_ablations.py

@@ -0,0 +1,218 @@
+# /// script
+# requires-python = ">=3.10"
+# dependencies = ["numpy", "matplotlib", "seaborn"]
+# ///
+"""S4 comparison figure: metacal breakdown-grid ablations vs baseline.
+
+Reads the v2 breakdown_grid JSON (see ``run_breakdown_grid.py`` / ``plot_breakdown_grid.py``)
+for a baseline run and one or more ablation runs, and draws a 3 (S/N band) x 3 (m1, m2, c1)
+grid of panels vs resolution ratio, baseline in dark gray and ablations in husl hues, so a
+reader can see at a glance whether any ablation moves the bias trends outside the tolerance
+bands. Cells with a null estimator (the degeneracy guard fired) are skipped -- there is no
+value to plot. Grids may be partial (a single cell, as in the S4 smoke tests); the script
+draws whatever is present. Seaborn is not in the ShapePipe sif, so this runs on the host:
+
+    uv run scripts/python/plot_s4_ablations.py \\
+        --baseline <run>/results/baseline/breakdown_grid \\
+        --ablations <run>/results/wcscen/breakdown_grid <run>/results/jac/breakdown_grid \\
+        --labels "WCS centroid" "Jacobian" \\
+        --output <figures>
+
+Add ``--paired <ablation-dir>`` when an ablation's images are bit-identical to baseline
+(e.g. the WCS-centroid arm) to also recompute a paired Delta-m1 per cell straight from the
+raw metacal records in ``<dir>/cells/*.npz``, which cancels shape noise and gives an error
+roughly two orders of magnitude tighter than the independent-run comparison.
+"""
+import argparse
+import json
+from pathlib import Path
+
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.lines import Line2D
+from matplotlib.ticker import NullFormatter, ScalarFormatter
+import seaborn as sns
+
+BANDS = [("high", r"S/N $\approx$ 1000"), ("moderate", r"S/N $\approx$ 50"),
+         ("low", r"S/N $\approx$ 15")]
+COLUMNS = [("m1", r"$m_1$", 5e-3), ("m2", r"$m_2$", 5e-3), ("c1", r"$c_1$", 1e-3)]
+GAMMA = 0.02
+
+
+def load_grid(grid_dir):
+    """{band: {res_ratio: cell}} for a breakdown_grid_v2.json, skipping cells with any
+    null estimator -- a null means the degeneracy guard fired and there is nothing to plot."""
+    d = json.load(open(Path(grid_dir) / "breakdown_grid_v2.json"))
+    by_band = {}
+    for c in d["cells"]:
+        if c["m1"] is None or c["m2"] is None or c["c1"] is None:
+            continue
+        by_band.setdefault(c["s2n_band"], {})[c["res_ratio"]] = c
+    return by_band, d["meta"]
+
+
+def series(by_band, band, key):
+    """(res_ratio, value, err) for one band and one estimator key, sorted by res_ratio."""
+    cells = by_band.get(band, {})
+    rows = sorted(cells.values(), key=lambda c: c["res_ratio"])
+    r = np.array([c["res_ratio"] for c in rows])
+    val = np.array([c[key] for c in rows], dtype=float)
+    err = np.array([c[f"{key}_err"] for c in rows], dtype=float)
+    return r, val, err
+
+
+def draw_panel(ax, runs, band, key, tol, pal, dodge_frac=0.015):
+    """One (band, key) panel: baseline dark gray, ablations husl-hued, slight x-dodge so
+    error bars at the same res_ratio don't overlap."""
+    n = len(runs)
+    for i, (label, by_band, color) in enumerate(runs):
+        r, val, err = series(by_band, band, key)
+        if r.size == 0:
+            continue
+        dodge = (i - (n - 1) / 2) * dodge_frac
+        ax.errorbar(r * (1 + dodge), val, yerr=err, color=color, marker="o", ms=6,
+                    lw=1.6, capsize=3, elinewidth=1.2, label=label, zorder=3 if i == 0 else 2)
+    ax.axhspan(-tol, tol, color="0.75", alpha=0.35, zorder=0)
+    ax.axhline(0, color="0.5", lw=0.8, zorder=0)
+    ax.set_xscale("log")
+    ax.xaxis.set_minor_formatter(NullFormatter())  # dense log-minor labels collide otherwise
+    ax.xaxis.set_major_formatter(ScalarFormatter())
+
+
+def make_figure(baseline, ablations, labels, out_png):
+    base_by_band, _ = load_grid(baseline)
+    abl_grids = [load_grid(d)[0] for d in ablations]
+    pal = sns.color_palette("husl", len(ablations))
+    runs = [("baseline", base_by_band, "0.25")] + list(zip(labels, abl_grids, pal))
+
+    sns.set_theme(style="ticks", context="talk")
+    fig, axes = plt.subplots(len(BANDS), len(COLUMNS), figsize=(4.4 * len(COLUMNS), 3.6 * len(BANDS)),
+                              sharex="col", squeeze=False)
+
+    for row, (band, band_lab) in enumerate(BANDS):
+        for col, (key, key_lab, tol) in enumerate(COLUMNS):
+            ax = axes[row][col]
+            draw_panel(ax, runs, band, key, tol, pal)
+            if row == 0:
+                ax.set_title(key_lab)
+            if col == 0:
+                ax.set_ylabel(f"{band_lab}\nvalue")
+            if row == len(BANDS) - 1:
+                ax.set_xlabel(r"resolution ratio  $r$")
+
+    handles = [Line2D([0], [0], color=color, marker="o", lw=1.6, label=label)
+               for label, _, color in runs]
+    fig.legend(handles=handles, frameon=False, loc="upper center",
+               ncol=min(len(runs), 5), bbox_to_anchor=(0.5, 1.02))
+    sns.despine(fig)
+    fig.tight_layout()
+    fig.savefig(out_png, dpi=150, bbox_inches="tight")
+    print("wrote", out_png)
+
+
+def print_ablation_table(baseline, ablations, labels):
+    """Per (ablation, band): max |m1|/|m2| over resolved cells (res_ratio >= 0.35), its
+    error, and the max |Delta| vs baseline with combined error and significance."""
+    base_by_band, _ = load_grid(baseline)
+    for label, abl_dir in zip(labels, ablations):
+        abl_by_band, _ = load_grid(abl_dir)
+        print(f"\n=== {label} ===")
+        header = f"{'band':<10}{'stat':<6}{'max|val|':>10}{'err':>10}{'max|Δ|':>10}{'σ_Δ':>10}{'|Δ|/σ':>8}"
+        print(header)
+        for band, _ in BANDS:
+            for key in ("m1", "m2"):
+                base_r, base_val, base_err = series(base_by_band, band, key)
+                abl_r, abl_val, abl_err = series(abl_by_band, band, key)
+                common = sorted(set(base_r[base_r >= 0.35]) & set(abl_r[abl_r >= 0.35]))
+                if not common:
+                    continue
+                base_idx = {r: i for i, r in enumerate(base_r)}
+                abl_idx = {r: i for i, r in enumerate(abl_r)}
+                i_max = max(range(len(abl_r)), key=lambda i: abs(abl_val[i]) if abl_r[i] >= 0.35 else -1)
+                max_val, max_err = abs(abl_val[i_max]), abl_err[i_max]
+
+                deltas = np.array([abl_val[abl_idx[r]] - base_val[base_idx[r]] for r in common])
+                sigmas = np.array([np.sqrt(abl_err[abl_idx[r]]**2 + base_err[base_idx[r]]**2)
+                                   for r in common])
+                j = int(np.argmax(np.abs(deltas)))
+                max_delta, sigma_delta = abs(deltas[j]), sigmas[j]
+                sig = max_delta / sigma_delta if sigma_delta > 0 else np.inf
+                flag = " ***" if sig > 3 else ""
+                print(f"{band:<10}{key:<6}{max_val:>10.4f}{max_err:>10.4f}"
+                      f"{max_delta:>10.4f}{sigma_delta:>10.4f}{sig:>8.2f}{flag}")
+
+
+def paired_delta_m1(baseline_dir, ablation_dir, B=1000, seed=0):
+    """Paired Delta-m1 per cell: recompute m1 directly from raw metacal records so shape
+    noise -- common to both runs, since the images are bit-identical -- cancels. Bootstraps
+    over seeds with one shared index resample per replicate across all arrays."""
+    base_cells = sorted(Path(baseline_dir, "cells").glob("*.npz"))
+    abl_cells = sorted(Path(ablation_dir, "cells").glob("*.npz"))
+    base_by_name = {p.name: p for p in base_cells}
+    abl_by_name = {p.name: p for p in abl_cells}
+    common_names = sorted(set(base_by_name) & set(abl_by_name))
+    if not common_names:
+        raise ValueError(f"no matching cell files between {baseline_dir} and {ablation_dir}")
+
+    rng = np.random.default_rng(seed)
+    results = []
+    for name in common_names:
+        base = np.load(base_by_name[name])
+        abl = np.load(abl_by_name[name])
+        ok = (base["g1p__flags"] == 0) & (base["g1p__finite"]) \
+            & (base["g1m__flags"] == 0) & (base["g1m__finite"]) \
+            & (abl["g1p__flags"] == 0) & (abl["g1p__finite"]) \
+            & (abl["g1m__flags"] == 0) & (abl["g1m__finite"])
+        if not ok.any():
+            raise ValueError(f"no jointly-valid seeds in {name}")
+
+        base_g1p, base_g1m = base["g1p__g1"][ok], base["g1m__g1"][ok]
+        base_R11p, base_R11m = base["g1p__R11"][ok], base["g1m__R11"][ok]
+        abl_g1p, abl_g1m = abl["g1p__g1"][ok], abl["g1m__g1"][ok]
+
+        Rbar_base = 0.5 * (base_R11p.mean() + base_R11m.mean())
+        delta_m1 = ((abl_g1p - base_g1p).mean() - (abl_g1m - base_g1m).mean()) / (2 * GAMMA * Rbar_base)
+
+        n = ok.sum()
+        boot = np.empty(B)
+        for b in range(B):
+            idx = rng.integers(0, n, size=n)
+            Rbar_b = 0.5 * (base_R11p[idx].mean() + base_R11m[idx].mean())
+            boot[b] = ((abl_g1p[idx] - base_g1p[idx]).mean()
+                       - (abl_g1m[idx] - base_g1m[idx]).mean()) / (2 * GAMMA * Rbar_b)
+        results.append((name, float(delta_m1), float(boot.std())))
+    return results
+
+
+def print_paired_table(baseline_dir, ablation_dir):
+    print(f"\n=== paired Δm1: {ablation_dir} vs {baseline_dir} ===")
+    for name, delta_m1, err in paired_delta_m1(baseline_dir, ablation_dir):
+        print(f"{name:<24}Δm1 = {delta_m1:+.6f} ± {err:.6f}")
+
+
+def main():
+    p = argparse.ArgumentParser(description=__doc__.split("\n\n")[0])
+    p.add_argument("--baseline", required=True, help="baseline breakdown_grid dir.")
+    p.add_argument("--ablations", required=True, nargs="+", help="ablation breakdown_grid dirs.")
+    p.add_argument("--labels", required=True, nargs="+", help="one label per ablation dir.")
+    p.add_argument("--output", required=True, help="output figure directory.")
+    p.add_argument("--paired", default=None,
+                    help="ablation dir (must also be in --ablations) to additionally run "
+                         "the paired same-noise Δm1 estimator against --baseline.")
+    a = p.parse_args()
+
+    if len(a.ablations) != len(a.labels):
+        raise ValueError(f"--ablations has {len(a.ablations)} entries but --labels has "
+                          f"{len(a.labels)}; need one label per ablation dir.")
+
+    out_dir = Path(a.output)
+    out_dir.mkdir(parents=True, exist_ok=True)
+    make_figure(a.baseline, a.ablations, a.labels, out_dir / "s4_ablations.png")
+    print_ablation_table(a.baseline, a.ablations, a.labels)
+
+    if a.paired is not None:
+        print_paired_table(a.baseline, a.paired)
+
+
+if __name__ == "__main__":
+    main()