diff --git a/scripts/python/plot_breakdown_grid.py b/scripts/python/plot_breakdown_grid.py index 671dc6cf..c24ef1bc 100644 --- a/scripts/python/plot_breakdown_grid.py +++ b/scripts/python/plot_breakdown_grid.py @@ -2,17 +2,20 @@ # 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 /results/baseline/breakdown_grid/breakdown_grid.json \\ + --input /results/baseline/breakdown_grid/breakdown_grid_v2.json \\ --output /breakdown_grid.png """ import argparse @@ -20,24 +23,49 @@ 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") diff --git a/scripts/python/plot_s4_ablations.py b/scripts/python/plot_s4_ablations.py new file mode 100644 index 00000000..6c074385 --- /dev/null +++ b/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 /results/baseline/breakdown_grid \\ + --ablations /results/wcscen/breakdown_grid /results/jac/breakdown_grid \\ + --labels "WCS centroid" "Jacobian" \\ + --output + +Add ``--paired `` 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 ``/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() diff --git a/scripts/python/run_breakdown_grid.py b/scripts/python/run_breakdown_grid.py index 1bbe366b..7ec71ad9 100644 --- a/scripts/python/run_breakdown_grid.py +++ b/scripts/python/run_breakdown_grid.py @@ -1,5 +1,5 @@ #!/usr/bin/env python -"""Runnable resolution x noise breakdown grid for the ShapePipe digital twin. +"""Statistically rigorous resolution x noise breakdown grid (v2). The Tier-1 companion of ``run_mbias.py``. Where that recipe reports the multiplicative bias at a single well-resolved, high-S/N operating point, this @@ -7,20 +7,48 @@ tolerance and where it breaks down: galaxy size relative to the PSF (the resolution ratio ``r = gal_hlr / psf_fwhm``) and image noise (three S/N bands). -Each grid cell is one self-contained in-memory GalSim simulation: an exponential -galaxy convolved with a Moffat (beta = 2.5) PSF, sampled at 0.1857" pixels, with -a known shear injected. The stamp is passed through the *live* pipeline shear -estimator (``do_ngmix_metacal``) -- the same code the real survey uses -- and the -recovered shear, the metacal response ``R``, and the additive bias ``c`` are read -back. Averaging over ``n_seeds`` per cell gives ``m``, ``R11`` and ``c1`` for that -cell. Estimator: *ratio of averages* -- ``m = / / gamma - 1``, -``c1 = / `` -- the same shape estimator the twin grounds against. +**What changed from v1.** v1 ran a single-signed arm per axis and reported bare +``m = //gamma - 1`` with no uncertainty; at moderate/low S/N the per-seed +shear noise gives sigma_m ~ 0.5, so v1's scary S/N trends were *unmeasured*. v2 +makes every number carry an honest bootstrap error and cancels measurement noise +by pairing signed arms that share an identical noise realisation: + +* **Paired shear response.** Six arms per seed -- (g1+, g1-), (g2+, g2-), + (psf+, psf-) -- driven so the noise draw and sub-pixel shift are BIT-IDENTICAL + across arms of a seed (same ``RandomState(seed+1000)`` into ``make_data``; + galsim ``.shear()`` adds no rng draw). The +/- difference then cancels the + shared shape noise that dominated v1. ``m1 = ( - )/(2 gamma + Rbar11) - 1`` with ``Rbar11`` the mean R11 over both g1 arms. +* **Additive bias from a PSF-shear pair.** ``c1 = ( - )/(2 + Rbar11_psf)``, leakage ``alpha = c1 / 0.05``. c2 (expected 0) is a null check. +* **Bootstrap-over-seeds errors** (B=1000) on every cell estimator. +* **Degeneracy guard.** When ``|Rbar| < 5 se(Rbar)`` the m is set null with + ``degenerate: true`` rather than dividing by a response consistent with zero. +* **Legacy estimator** (v1's ratio-of-averages) carried alongside for continuity. + +**Production-parity knobs (S4).** Three CLI knobs turn the clean-room grid into +production-settings ablations, each defaulting to the clean baseline: + +* ``--centroid-source {hsm,wcs}`` — "wcs" is the production default since + ``31ae736c``: the ngmix jacobian centre comes from the catalogue sky position + projected through the epoch WCS instead of HSM re-centering. The harness + fabricates the astrometry truth (a TAN WCS whose CD is the drawing jacobian; + ra/dec evaluated through that same WCS at the object's true pixel position) — + the toy analogue of coadd (x,y) -> (ra,dec) -> epoch pixel. +* ``--wcs-g1/--wcs-g2/--wcs-theta-deg`` — shear/rotate the drawing WCS jacobian + away from a pure pixel scale (the esheldon/ngmix#72 sensitivity axis, where a + mishandled jacobian produced m ~ -0.2 at wcs_g1=0.1). +* ``--n-epochs`` — pre-existing; n_epochs=2 matches the old CI operating point. + +Estimator math is exact; the finite-difference response step is hardcoded to +``do_ngmix_metacal``'s mainline value (Spec-Code Invariant: if that line moves, +``STEP`` moves with it). Runs ONLY inside the ShapePipe apptainer container with the live source tree bound over the baked copy -- same invocation as ``run_mbias.py`` (two binds: the live shapepipe src AND the live cs_util, since mainline ngmix.py imports -``cs_util.size``). ~2 min on a candide login node for the default 8x3 grid, -16 seeds/cell: +``cs_util.size``). Cells are embarrassingly parallel (``--jobs``); the full +24-cell production grid is driven by the orchestrator on SLURM, not here. apptainer exec \\ --bind /home,/scratch,/automnt,/n17data,/n23data1,/n09data \\ @@ -30,17 +58,25 @@ python scripts/python/run_breakdown_grid.py \\ --output /results/baseline/breakdown_grid -Emits ``/breakdown_grid.json``; the column figure is drawn from that JSON -by the host-side ``plot_breakdown_grid.py`` (seaborn, not in the sif). +Emits ``/breakdown_grid_v2.json`` (meta + provenance + per-cell summaries +with all errors) plus a mirror ``breakdown_grid.json`` for the status-page reader, +and per-cell ``/cells/_r.npz`` holding the raw per-(seed,arm) +arrays so a cell can be re-analysed without refitting. The column figure is drawn +from the JSON by the host-side ``plot_breakdown_grid.py`` (seaborn, not in the sif). """ import argparse import json import os +import subprocess import time +from concurrent.futures import ProcessPoolExecutor +import astropy.wcs import galsim import ngmix import numpy as np +from modopt.math.stats import sigma_mad +from ngmix.gexceptions import BootGalFailure, BootPSFFailure from shapepipe.modules.ngmix_package.ngmix import ( Postage_stamp, @@ -53,75 +89,276 @@ # do_ngmix_metacal's mainline value (ngmix.py); the Spec-Code Invariant says if # that line moves, this must move with it. STEP = 0.01 +# PSF-shear amplitude for the additive-bias (leakage) arms; c1 = dgn / (2*this). +PSF_SHEAR = 0.05 +# MegaCam pixel scale (arcsec/pixel): make_data's default and the prior width. +PIXEL_SCALE = 0.1857 # The three noise levels map to nominal S/N ~ 1000 / 50 / 15, calibrated at the # grid midpoint (r = 0.5). S/N drifts with galaxy size, so mean_s2n is recorded # per cell rather than trusted from the label. NOISE = {"high": 0.08, "moderate": 1.6, "low": 5.0} +# Default seeds per band: moderate/low need many more to beat their shear noise +# down even after pairing. CLI-overridable. +N_SEEDS_DEFAULT = {"high": 100, "moderate": 400, "low": 1600} +# The six signed arms: (shear, psf_shear). g1/g2 probe multiplicative response; +# psf probes additive leakage. Arms of one seed share a noise realisation. +ARMS = { + "g1p": ((0.02, 0.0), (0.0, 0.0)), "g1m": ((-0.02, 0.0), (0.0, 0.0)), + "g2p": ((0.0, 0.02), (0.0, 0.0)), "g2m": ((0.0, -0.02), (0.0, 0.0)), + "psfp": ((0.0, 0.0), (PSF_SHEAR, 0.0)), "psfm": ((0.0, 0.0), (-PSF_SHEAR, 0.0)), +} +# The per-(seed,arm) scalars persisted to the .npz for refit-free re-analysis. +RECORD_KEYS = ["g1", "g2", "R11", "R12", "R21", "R22", "s2n", "T", "T_psf_orig", + "flags", "finite", "sigma_mad_est", "noise_true"] + + +def build_wcs(g1, g2, theta_deg): + """Uniform drawing WCS from the S4 parity knobs. + + All-zero knobs return None -> make_data's bit-identical legacy PixelScale + path. Otherwise: shear the pixel scale (area-preserving, det = scale**2 -- + the esheldon/ngmix#72 construction), then rotate by theta_deg.""" + if g1 == g2 == theta_deg == 0.0: + return None + jac = galsim.ShearWCS(PIXEL_SCALE, galsim.Shear(g1=g1, g2=g2)).jacobian() + th = np.deg2rad(theta_deg) + rot = np.array([[np.cos(th), -np.sin(th)], [np.sin(th), np.cos(th)]]) + m = rot @ np.array([[jac.dudx, jac.dudy], [jac.dvdx, jac.dvdy]]) + return galsim.JacobianWCS(m[0, 0], m[0, 1], m[1, 0], m[1, 1]) + + +def fabricate_astrometry(stamp, centers, jacob, img_size): + """Populate stamp.wcs/ra/dec -- the truth the "wcs" centroid path consumes. + Production reads the catalogue sky position and projects it through the + epoch WCS; the toy analogue is a TAN WCS whose CD matrix IS the drawing + jacobian (deg/pixel), with ra/dec evaluated through that same WCS at the + object's true pixel position. The round trip toWorld -> toImage is then + self-cancelling (TAN nonlinearity over half a stamp ~ 5e-9 px), so the + centroid the estimator recovers is the truth -- modulo the same + round-to-nearest-pixel logic production applies. Epochs + get independent sub-pixel shifts, so per-epoch ra/dec differ; the estimator + reads each epoch independently, making this equivalent to one sky position + seen through slightly different epoch WCSs, as in production.""" + # Even stamps put the galsim centre on a half-integer, so np.round() in the + # consumer (ngmix.py "wcs" path) lands one pixel off for EVERY object -- a + # systematic half-pixel origin error that would read as a production + # pathology. Production VIGNETs are odd; hold the harness to it. + assert img_size % 2 == 1, "--centroid-source wcs requires odd --img-size" + w = astropy.wcs.WCS(naxis=2) + w.wcs.ctype = ["RA---TAN", "DEC--TAN"] + w.wcs.crpix = [(img_size + 1) / 2, (img_size + 1) / 2] + w.wcs.crval = [150.0, 0.0] + w.wcs.cd = np.array( + [[jacob.dudx, jacob.dudy], [jacob.dvdx, jacob.dvdy]] + ) / 3600.0 + g_wcs = galsim.fitswcs.AstropyWCS(wcs=w) + for cen in centers: + world = g_wcs.toWorld(cen) + stamp.wcs.append(w) + stamp.ra.append(world.ra / galsim.degrees) + stamp.dec.append(world.dec / galsim.degrees) -def one_stamp(noise, gal_hlr, psf_fwhm, img_size, n_epochs, shear, psf_shear, seed): - """One injected-truth realisation -> (metacal resdict, original-fit dict).""" + +def one_stamp(noise, gal_hlr, psf_fwhm, img_size, n_epochs, shear, psf_shear, + seed, true_noise, wcs=None, centroid_source="hsm"): + """One injected-truth realisation -> per-arm record dict. + + Uses ``RandomState(seed)`` for the fit rng and ``RandomState(seed+1000)`` for + make_data, so every arm of a given seed sees the IDENTICAL noise draw and + sub-pixel shift (galsim ``.shear()`` consumes no rng) -- the property the + paired estimator rests on. ``true_noise`` routes the per-pixel true-inverse- + variance path (see ``--true-noise``). ``wcs``/``centroid_source`` are the + S4 parity knobs (see module docstring).""" rng = np.random.RandomState(seed) - prior = get_prior(0.1857, rng) - gals, psfs, _, weights, flags, jacobs = make_data( + prior = get_prior(PIXEL_SCALE, rng) + gals, psfs, _, weights, flags, jacobs, centers = make_data( rng=np.random.RandomState(seed + 1000), shear=shear, noise=noise, n_epochs=n_epochs, gal_hlr=gal_hlr, psf_fwhm=psf_fwhm, img_size=img_size, - psf_shear=psf_shear, + psf_shear=psf_shear, wcs=wcs, return_centers=True, ) stamp = Postage_stamp(bkg_sub=False, megacam_flip=False) stamp.gals, stamp.psfs, stamp.weights, stamp.flags, stamp.jacobs = ( gals, psfs, weights, flags, jacobs, ) - res = do_ngmix_metacal(stamp, prior, 1.0, rng) - return res.resdict, res.orig + if centroid_source == "wcs": + fabricate_astrometry(stamp, centers, jacobs[0], img_size) + if true_noise: + # True per-pixel sigma into bkg_rms -> prepare_ngmix_weights takes the + # inverse-variance path and metacal fixnoise gets the true sigma. + stamp.bkg_rms = [np.full((img_size, img_size), noise) for _ in gals] + try: + res = do_ngmix_metacal(stamp, prior, 1.0, rng, + centroid_source=centroid_source) + except (BootPSFFailure, BootGalFailure): + # Production treats a bootstrap failure as a flagged object; the + # pair-drop policy drops the seed and fail_frac records the rate (a + # SYSTEMATIC failure mode still surfaces loudly there). Without this, + # one unlucky fit in ~1e5 kills the entire grid (job 808259). + return dict( + g1=np.nan, g2=np.nan, R11=np.nan, R12=np.nan, R21=np.nan, + R22=np.nan, s2n=np.nan, T=np.nan, T_psf_orig=np.nan, + flags=-1, finite=False, + sigma_mad_est=float(sigma_mad(gals[0])), noise_true=float(noise), + ) + d, orig = res.resdict, res.orig + ns = d["noshear"] + finite = bool(np.all(np.isfinite(ns["g"]))) + return dict( + g1=ns["g"][0] if finite else np.nan, g2=ns["g"][1] if finite else np.nan, + R11=response(d, 1, 1), R12=response(d, 1, 2), + R21=response(d, 2, 1), R22=response(d, 2, 2), + s2n=ns.get("s2n", np.nan), T=ns.get("T", np.nan), + T_psf_orig=orig["T_psf"], flags=int(ns["flags"]), finite=finite, + # sigma_mad on the epoch-0 galaxy image is the noise the estimator would + # infer from the data; its ratio to noise_true is the contamination + # diagnostic (>1 when object flux inflates the MAD, the moderate/low path). + sigma_mad_est=float(sigma_mad(gals[0])), noise_true=float(noise), + ) def response(d, i, j): """R_ij = (g_i[p] - g_i[m]) / (2*step).""" p, m = {1: "1p", 2: "2p"}[j], {1: "1m", 2: "2m"}[j] - return (d[p]["g"][i - 1] - d[m]["g"][i - 1]) / (2 * STEP) - - -def harvest_cell(noise, gal_hlr, psf_fwhm, img_size, n_epochs, gamma, seeds): - """Ratio-of-averages m, R11, c1 for one (resolution, noise) cell. - - Three arms per seed: a g1-shear arm (m1, R11), a g2-shear arm (m2, R22), and - a PSF-shear-only arm (the additive-bias null c1). Failed fits (nonzero flags - or non-finite shear) are dropped and counted in fail_frac.""" - g1ns, g2ns, R11, R22, c1ns, cR, s2n, trat = ([] for _ in range(8)) - nfail = ntot = 0 - arms = { - "g1": ((gamma, 0.0), (0.0, 0.0)), - "g2": ((0.0, gamma), (0.0, 0.0)), - "psf": ((0.0, 0.0), (0.05, 0.0)), - } + gp, gm = d[p]["g"], d[m]["g"] + if not (np.all(np.isfinite(gp)) and np.all(np.isfinite(gm))): + return np.nan + return (gp[i - 1] - gm[i - 1]) / (2 * STEP) + + +BOOTSTRAP_B = 1000 + + +def _boot_idx(n, B, rng): + """One seed-index array per bootstrap replicate. + + Every quantity inside a replicate must be computed from the SAME resampled + seed set: resampling each arm independently silently destroys the ±γ + pairing and inflates the reported error to the unpaired level (the bug this + replaces read paired σ_m = 0.23 where the per-seed pairs give 0.009).""" + return rng.integers(0, n, size=(B, n)) + + +def harvest_cell(args): + """Six-arm paired estimators for one (resolution, noise) cell. + + ``args`` is a flat tuple so the cell is a single ProcessPoolExecutor task. + Returns (summary dict, raw per-(seed,arm) arrays for the .npz).""" + (cell_index, res_ratio, band, noise, gal_hlr, psf_fwhm, img_size, n_epochs, + gamma, n_seeds, true_noise, wcs_g1, wcs_g2, wcs_theta_deg, + centroid_source) = args + # Built per worker process: galsim WCS objects stay out of the task tuple. + wcs = build_wcs(wcs_g1, wcs_g2, wcs_theta_deg) + # Seeds independent per cell, never recycled across cells. + seeds = [10_000_000 + cell_index * 100_000 + s for s in range(n_seeds)] + + # Per-arm columns of per-seed records (aligned by seed index). + cols = {arm: {k: [] for k in RECORD_KEYS} for arm in ARMS} for seed in seeds: - for arm, (shear, psf_shear) in arms.items(): - ntot += 1 - d, orig = one_stamp(noise, gal_hlr, psf_fwhm, img_size, n_epochs, - shear, psf_shear, seed) - ns = d["noshear"] - if ns["flags"] != 0 or not np.all(np.isfinite(ns["g"])): - nfail += 1 - continue - if arm == "g1": - g1ns.append(ns["g"][0]); R11.append(response(d, 1, 1)) - s2n.append(ns["s2n"]); trat.append(ns["T"] / orig["T_psf"]) - elif arm == "g2": - g2ns.append(ns["g"][1]); R22.append(response(d, 2, 2)) - else: - c1ns.append(ns["g"][0]); cR.append(response(d, 1, 1)) - mR11, mR22 = np.mean(R11), np.mean(R22) - return dict( - gal_hlr=gal_hlr, psf_fwhm=psf_fwhm, noise=noise, - m1=float(np.mean(g1ns) / mR11 / gamma - 1.0), - m2=float(np.mean(g2ns) / mR22 / gamma - 1.0), - c1=float(np.mean(c1ns) / np.mean(cR)), - R11=float(mR11), R22=float(mR22), - mean_s2n=float(np.mean(s2n)), mean_T_ratio=float(np.mean(trat)), - fail_frac=float(nfail / ntot), n_seeds=len(seeds), + for arm, (shear, psf_shear) in ARMS.items(): + rec = one_stamp(noise, gal_hlr, psf_fwhm, img_size, n_epochs, + shear, psf_shear, seed, true_noise, + wcs=wcs, centroid_source=centroid_source) + for k in RECORD_KEYS: + cols[arm][k].append(rec[k]) + arr = {arm: {k: np.asarray(v) for k, v in cols[arm].items()} for arm in ARMS} + + def ok(arm): + """Boolean mask: this arm's fit succeeded (flags==0 and finite g).""" + return (arr[arm]["flags"] == 0) & arr[arm]["finite"].astype(bool) + + def pair_mask(a, b): + """Both arms of the pair usable for this seed (pair-drop policy).""" + return ok(a) & ok(b) + + fail_frac = {arm: float(np.mean(~ok(arm))) for arm in ARMS} + rng = np.random.default_rng(cell_index) # deterministic bootstrap per cell + + summary = dict( + cell_index=cell_index, res_ratio=float(res_ratio), s2n_band=band, + gal_hlr=float(gal_hlr), psf_fwhm=float(psf_fwhm), noise=float(noise), + n_seeds=int(n_seeds), fail_frac=fail_frac, ) + # --- multiplicative bias, paired, per axis --- + for axis, (ap, am, ri, comp) in { + 1: ("g1p", "g1m", "R11", "g1"), 2: ("g2p", "g2m", "R22", "g2"), + }.items(): + mask = pair_mask(ap, am) + npair = int(mask.sum()) + gp, gm = arr[ap][comp][mask], arr[am][comp][mask] + Rp, Rm = arr[ap][ri][mask], arr[am][ri][mask] + Rbar = float(0.5 * (Rp.mean() + Rm.mean())) + # Joint bootstrap: one seed resample per replicate, all arms together — + # the pairing (identical noise across ±γ) must survive the resampling. + idx = _boot_idx(npair, BOOTSTRAP_B, rng) + bR = 0.5 * (Rp[idx].mean(axis=1) + Rm[idx].mean(axis=1)) + seR = float(np.std(bR)) + degenerate = abs(Rbar) < 5 * seR + bm_boot = (gp[idx].mean(axis=1) - gm[idx].mean(axis=1)) / (2 * gamma * bR) - 1.0 + m = None if degenerate else float((np.mean(gp) - np.mean(gm)) / (2 * gamma * Rbar) - 1.0) + m_err = None if degenerate else float(np.std(bm_boot)) + summary[f"m{axis}"] = m + summary[f"m{axis}_err"] = m_err + summary[f"m{axis}_degenerate"] = bool(degenerate) + summary[f"m{axis}_n_pairs"] = npair + summary[f"Rbar{axis}{axis}"] = Rbar + summary[f"Rbar{axis}{axis}_err"] = seR + + # --- additive bias c1/c2 from the PSF-shear pair, leakage alpha --- + mask = pair_mask("psfp", "psfm") + npair = int(mask.sum()) + Rpsf_p, Rpsf_m = arr["psfp"]["R11"][mask], arr["psfm"]["R11"][mask] + Rbar_psf = float(0.5 * (Rpsf_p.mean() + Rpsf_m.mean())) + idx_c = _boot_idx(npair, BOOTSTRAP_B, rng) + bRpsf = 0.5 * (Rpsf_p[idx_c].mean(axis=1) + Rpsf_m[idx_c].mean(axis=1)) + seR_psf = float(np.std(bRpsf)) + degenerate_c = abs(Rbar_psf) < 5 * seR_psf + for axis, comp in ((1, "g1"), (2, "g2")): + gp, gm = arr["psfp"][comp][mask], arr["psfm"][comp][mask] + bc = (gp[idx_c].mean(axis=1) - gm[idx_c].mean(axis=1)) / (2 * bRpsf) + c = None if degenerate_c else float((np.mean(gp) - np.mean(gm)) / (2 * Rbar_psf)) + c_err = None if degenerate_c else float(np.std(bc)) + summary[f"c{axis}"] = c + summary[f"c{axis}_err"] = c_err + summary["c_degenerate"] = bool(degenerate_c) + summary["c_n_pairs"] = npair + summary["Rbar_psf"] = Rbar_psf + summary["Rbar_psf_err"] = seR_psf + summary["alpha"] = None if summary["c1"] is None else float(summary["c1"] / PSF_SHEAR) + summary["alpha_err"] = None if summary["c1_err"] is None else float(summary["c1_err"] / PSF_SHEAR) + + # --- legacy ratio-of-averages (v1 continuity): m = //gamma - 1 --- + for axis, (ap, ri, comp) in {1: ("g1p", "R11", "g1"), 2: ("g2p", "R22", "g2")}.items(): + m = ok(ap) + g, R = arr[ap][comp][m], arr[ap][ri][m] + idx_l = _boot_idx(len(g), BOOTSTRAP_B, rng) # joint: g and R share the seed resample + summary[f"m{axis}_legacy"] = float(np.mean(g) / np.mean(R) / gamma - 1.0) + summary[f"m{axis}_legacy_err"] = float( + np.std(g[idx_l].mean(axis=1) / R[idx_l].mean(axis=1) / gamma - 1.0)) + + # --- diagnostics --- + g1ok = ok("g1p") + summary["mean_s2n"] = float(np.nanmean(arr["g1p"]["s2n"][g1ok])) + summary["mean_T_ratio"] = float(np.nanmean( + arr["g1p"]["T"][g1ok] / arr["g1p"]["T_psf_orig"][g1ok])) + # contamination diagnostic: >1 means the data-inferred sigma exceeds the true + # noise because object flux inflates the whole-stamp MAD. + summary["mean_sigma_ratio"] = float(np.mean( + arr["g1p"]["sigma_mad_est"] / arr["g1p"]["noise_true"])) + + return summary, {f"{arm}__{k}": arr[arm][k] for arm in ARMS for k in RECORD_KEYS} + + +def _git_describe(): + try: + return subprocess.check_output( + ["git", "describe", "--tags", "--always", "--dirty"], + cwd=os.path.dirname(os.path.abspath(__file__)), text=True, + ).strip() + except Exception: + return "unknown" + def main(): p = argparse.ArgumentParser(description=__doc__.split("\n\n")[0]) @@ -130,45 +367,93 @@ def main(): p.add_argument("--psf-fwhm", type=float, default=0.55, help="Moffat PSF FWHM (arcsec).") p.add_argument("--img-size", type=int, default=51) p.add_argument("--n-epochs", type=int, default=1) - p.add_argument("--n-seeds", type=int, default=16, help="realisations averaged per cell.") + p.add_argument("--n-seeds", type=json.loads, default=None, + help='JSON per-band seed counts, e.g. \'{"high":100,"moderate":400,"low":1600}\'.') p.add_argument("--n-res", type=int, default=8, help="resolution-ratio grid points.") p.add_argument("--res-lo", type=float, default=0.15) p.add_argument("--res-hi", type=float, default=1.2) + p.add_argument("--bands", nargs="+", default=list(NOISE), choices=list(NOISE), + help="subset of S/N bands to run (default all three).") + p.add_argument("--res-index", type=int, nargs="+", default=None, + help="subset of resolution-grid indices to run (default all).") + p.add_argument("--true-noise", action="store_true", + help="route the per-pixel true-inverse-variance weight path.") + p.add_argument("--centroid-source", choices=["hsm", "wcs"], default="hsm", + help='ngmix jacobian centre: "hsm" (legacy re-centering) or ' + '"wcs" (production default since 31ae736c; the harness ' + "fabricates the astrometry truth).") + p.add_argument("--wcs-g1", type=float, default=0.0, + help="drawing-WCS jacobian shear g1 (ngmix#72 axis).") + p.add_argument("--wcs-g2", type=float, default=0.0, + help="drawing-WCS jacobian shear g2.") + p.add_argument("--wcs-theta-deg", type=float, default=0.0, + help="drawing-WCS jacobian rotation (degrees).") + p.add_argument("--jobs", type=int, default=4, help="parallel cells (ProcessPoolExecutor).") a = p.parse_args() - seeds = list(range(a.n_seeds)) + n_seeds = a.n_seeds or N_SEEDS_DEFAULT res_grid = np.geomspace(a.res_lo, a.res_hi, a.n_res) + res_idx = a.res_index if a.res_index is not None else list(range(a.n_res)) + # cell_index enumerates the FULL 8x3 grid in (res, band) order so seeds are + # stable regardless of which subset a given run selects. + all_cells = [(ri, bi) for ri in range(a.n_res) for bi, _ in enumerate(NOISE)] + cell_of = {(ri, bi): k for k, (ri, bi) in enumerate(all_cells)} + bands_ordered = list(NOISE) + + tasks = [] + for ri in res_idx: + for band in a.bands: + bi = bands_ordered.index(band) + cell_index = cell_of[(ri, bi)] + gal_hlr = float(res_grid[ri] * a.psf_fwhm) + tasks.append((cell_index, res_grid[ri], band, NOISE[band], gal_hlr, + a.psf_fwhm, a.img_size, a.n_epochs, a.gamma, + n_seeds[band], a.true_noise, a.wcs_g1, a.wcs_g2, + a.wcs_theta_deg, a.centroid_source)) + + os.makedirs(os.path.join(a.output, "cells"), exist_ok=True) cells, t0 = [], time.time() - for r_ratio in res_grid: - gal_hlr = float(r_ratio * a.psf_fwhm) - for band, noise in NOISE.items(): - cell = harvest_cell(noise, gal_hlr, a.psf_fwhm, a.img_size, - a.n_epochs, a.gamma, seeds) - cell.update(res_ratio=float(r_ratio), s2n_band=band) - cells.append(cell) - print(f"r={r_ratio:.2f} {band:9s} s2n={cell['mean_s2n']:7.1f} " - f"m1={cell['m1']:+.4f} m2={cell['m2']:+.4f} " - f"c1={cell['c1']:+.5f} R11={cell['R11']:.3f} " - f"fail={cell['fail_frac']:.2f}", flush=True) + with ProcessPoolExecutor(max_workers=a.jobs) as ex: + for summary, raw in ex.map(harvest_cell, tasks): + cells.append(summary) + np.savez_compressed( + os.path.join(a.output, "cells", + f"{summary['s2n_band']}_r{summary['res_ratio']:.3f}.npz"), + **raw) + m1 = " null" if summary["m1"] is None else f"{summary['m1']:+.4f}" + m1e = " " if summary["m1_err"] is None else f"{summary['m1_err']:.4f}" + c1 = " null" if summary["c1"] is None else f"{summary['c1']:+.5f}" + print(f"r={summary['res_ratio']:.2f} {summary['s2n_band']:9s} " + f"s2n={summary['mean_s2n']:7.1f} m1={m1}+-{m1e} c1={c1} " + f"Rbar11={summary['Rbar11']:.3f} " + f"deg={summary['m1_degenerate']} fail={summary['fail_frac']['g1p']:.2f}", + flush=True) out = dict( meta=dict( - step=STEP, gamma=a.gamma, psf_fwhm=a.psf_fwhm, img_size=a.img_size, - n_epochs=a.n_epochs, n_seeds=a.n_seeds, seeds=seeds, - res_grid=res_grid.tolist(), noise_map=NOISE, - m_tol_resolved=5e-3, runtime_s=time.time() - t0, + version=2, step=STEP, gamma=a.gamma, psf_shear=PSF_SHEAR, + psf_fwhm=a.psf_fwhm, img_size=a.img_size, n_epochs=a.n_epochs, + n_seeds=n_seeds, res_grid=res_grid.tolist(), noise_map=NOISE, + bands=a.bands, res_index=res_idx, true_noise=a.true_noise, + centroid_source=a.centroid_source, + wcs_shear=[a.wcs_g1, a.wcs_g2], wcs_theta_deg=a.wcs_theta_deg, + bootstrap_B=1000, degeneracy_k=5, m_tol_resolved=5e-3, + runtime_s=time.time() - t0, args=vars(a), ), provenance=dict( ngmix_version=ngmix.__version__, galsim_version=galsim.__version__, shapepipe_source=os.path.dirname(__import__("shapepipe").__file__), + git_describe=_git_describe(), ), cells=cells, ) - os.makedirs(a.output, exist_ok=True) - with open(os.path.join(a.output, "breakdown_grid.json"), "w") as fh: - json.dump(out, fh, indent=2) + # v2 JSON is primary; mirror to breakdown_grid.json so the status-page reader + # (build_status.py loads that exact name) stays green with no code change. + for fname in ("breakdown_grid_v2.json", "breakdown_grid.json"): + with open(os.path.join(a.output, fname), "w") as fh: + json.dump(out, fh, indent=2) print(f"DONE {len(cells)} cells in {time.time() - t0:.1f}s " - f"-> {a.output}/breakdown_grid.json") + f"-> {a.output}/breakdown_grid_v2.json (+ breakdown_grid.json mirror)") if __name__ == "__main__": diff --git a/src/shapepipe/testing/simulate.py b/src/shapepipe/testing/simulate.py index 8ab5a06f..4f263508 100644 --- a/src/shapepipe/testing/simulate.py +++ b/src/shapepipe/testing/simulate.py @@ -22,6 +22,8 @@ def make_data( pixel_scale=0.1857, img_size=201, psf_shear=(0.0, 0.0), + wcs=None, + return_centers=False, ): """Simulate an exponential galaxy with Moffat PSF. @@ -55,6 +57,23 @@ def make_data( guardrail needs: the PSF carries a shape that an unbiased deconvolution must remove from a round galaxy. The same sheared PSF is used both to convolve the object and as the PSF model stamp. + wcs : galsim.BaseWCS, optional + Uniform drawing WCS. Default None means ``galsim.PixelScale + (pixel_scale)`` with the legacy world-space sub-pixel shift draw — + bit-identical rng consumption and values to earlier versions of this + function. A non-trivial jacobian (rotation/shear — the ngmix#72 + sensitivity axis) switches the shift draw to PIXEL space (uniform + ±0.5 pixel per axis), guaranteeing the object lands within half a + pixel of the stamp centre — the contract that round-to-nearest-pixel + centroid logic (e.g. ngmix ``centroid_source="wcs"``) depends on. + Note ``psfs_sigmas`` are HSM pixel-unit sigmas of the *WCS-drawn* + PSF stamp; under a sheared WCS they are not circularly symmetric + measures. + return_centers : bool, optional + If True, also return per-epoch true object centres as + ``galsim.PositionD`` in galsim image coordinates (1-based; stamp + centre is ``(img_size + 1) / 2``). Default False keeps the legacy + 6-tuple return. Returns ------- @@ -64,19 +83,44 @@ def make_data( weights : list of numpy.ndarray flags : list of numpy.ndarray jacob_lists : list of galsim.BaseWCS + centers : list of galsim.PositionD, only if ``return_centers`` """ psf_noise = 1.0e-6 scale = pixel_scale - wcs = galsim.PixelScale(scale) + + if wcs is None: + wcs = galsim.PixelScale(scale) + + def draw_shift(): + # Legacy world-space draw — these exact lines (incl. dy-then-dx + # order) keep rng consumption and values bit-identical to the + # pre-wcs version of this function. + dy, dx = rng.uniform(low=-scale / 2, high=scale / 2, size=2) + return dx, dy, dx / scale, dy / scale + + else: + jac = wcs.jacobian() + transform = np.array([[jac.dudx, jac.dudy], [jac.dvdx, jac.dvdy]]) + + def draw_shift(): + # Pixel-space draw (same dy-then-dx order): under a sheared + # jacobian a world-space ±scale/2 box maps to a pixel-space + # parallelogram exceeding ±0.5 pixel, which would silently break + # rounded-centroid logic with off-by-one-pixel centres. + dpy, dpx = rng.uniform(low=-0.5, high=0.5, size=2) + dx, dy = transform @ (dpx, dpy) + return float(dx), float(dy), dpx, dpy if share_shift: - dy, dx = rng.uniform(low=-scale / 2, high=scale / 2, size=2) + dx, dy, dpx, dpy = draw_shift() - gals, psfs, psfs_sigmas, weights, flags, jacob_lists = [], [], [], [], [], [] + gals, psfs, psfs_sigmas, weights, flags, jacob_lists, centers = ( + [], [], [], [], [], [], [] + ) for epoch in range(n_epochs): if not share_shift: - dy, dx = rng.uniform(low=-scale / 2, high=scale / 2, size=2) + dx, dy, dpx, dpy = draw_shift() psf = galsim.Moffat(beta=2.5, fwhm=psf_fwhm).shear( g1=psf_shear[0], g2=psf_shear[1] @@ -101,5 +145,10 @@ def make_data( weights.append(im * 0 + 1.0 / noise ** 2) flags.append(im * 0) jacob_lists.append(wcs.jacobian()) + centers.append(galsim.PositionD( + (img_size + 1) / 2 + dpx, (img_size + 1) / 2 + dpy + )) + if return_centers: + return gals, psfs, psfs_sigmas, weights, flags, jacob_lists, centers return gals, psfs, psfs_sigmas, weights, flags, jacob_lists