From e9ee1f0eaa302d855cf4e05367bee75096e8e863 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Wed, 1 Jul 2026 00:16:04 +0200 Subject: [PATCH 1/2] ngmix: give the PSF observation a finite weight map (#749) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit make_ngmix_observation built psf_obs weightless, so ngmix defaulted to unit weight on the unit-flux PSF stamp and the galaxy GPriorBA(0.4) prior handed to the diagnostic PSF fitter swamped the likelihood. The recovered PSF shape collapsed toward the prior (zero ellipticity), corrupting the *_psf_orig diagnostic columns (#749), and the small finite-weight noise budget present in the pre-v2.0 esheldon/aguinot code was lost in the rewrite (the now-retired #774 — the same bug from the noise-budget angle). Build psf_obs with a flat weight = 1/PSF_NOISE**2 (PSF_NOISE = 1e-5, the esheldon/aguinot value Axel's reproduction used), forced float for parity with the galaxy weight — the PSF analogue of the galaxy weight built just below. Weight only; PSFEx/MCCD models already carry a little noise, so no explicit stamp-noise injection is needed (Axel's guidance on #749). The exact value is non-critical: validated flat across 1e-4..1e-6 on the twin. Calibration is untouched: metacal fits this same psf_obs with a prior-free AdmomFitter, which is insensitive to a flat weight's absolute scale, so the calibrated shear is invariant (verified bit-for-bit on the twin). Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_017KQn6mXoL4XMgscZDsaTkB --- src/shapepipe/modules/ngmix_package/ngmix.py | 23 +++- tests/module/test_ngmix_weight_validation.py | 110 +++++++++++++++++++ 2 files changed, 132 insertions(+), 1 deletion(-) diff --git a/src/shapepipe/modules/ngmix_package/ngmix.py b/src/shapepipe/modules/ngmix_package/ngmix.py index cc7ab87b..371264e2 100644 --- a/src/shapepipe/modules/ngmix_package/ngmix.py +++ b/src/shapepipe/modules/ngmix_package/ngmix.py @@ -23,6 +23,15 @@ METACAL_TYPES = ('noshear', '1p', '1m', '2p', '2m') +# Noise budget for the PSF observation's flat weight map (psf_wt = +# 1/PSF_NOISE**2). Mirrors the esheldon/aguinot pattern (sigma ~ 1e-5/1e-6); +# 1e-5 is the value Axel Guinot's #749 reproduction used. The fit is driven +# by the *relative* weighting of the likelihood vs the prior, so the exact +# value is non-critical once it is finite (validated on the digital twin: the +# recovered PSF shape/size are flat across 1e-4..1e-6). See +# make_ngmix_observation. +PSF_NOISE = 1e-5 + class MetacalResult(NamedTuple): """Return of :func:`do_ngmix_metacal`: the metacal fit plus both PSFs. @@ -1091,7 +1100,19 @@ def make_ngmix_observation( col=(psf.shape[1] - 1) / 2, wcs=wcs, ) - psf_obs = Observation(psf, jacobian=psf_jacob) + # A model fit always needs a weight: without one ngmix defaults to unit + # weights on a unit-flux PSF stamp, so the galaxy GPriorBA(0.4) prior + # handed to the diagnostic PSF fitter swamps the (tiny) likelihood and the + # recovered PSF shape collapses toward the prior (#749). PSFEx/MCCD models + # already carry a little noise, so a finite flat weight matching the + # esheldon/aguinot noise budget (psf_noise ~ 1e-5) suffices to restore the + # likelihood — no explicit stamp-noise injection needed (#774). This is the + # PSF analogue of the galaxy weight built just below; force float so an + # integer-typed stamp cannot truncate the weight (matching that build). + # Metacal's own PSF fit is prior-free (AdmomFitter) and so is insensitive + # to this flat weight's scale, leaving the calibration untouched. + psf_wt = np.full_like(psf, 1.0 / PSF_NOISE ** 2, dtype=float) + psf_obs = Observation(psf, weight=psf_wt, jacobian=psf_jacob) gal_masked, weight_map, noise_img = prepare_ngmix_weights( gal, weight, flag, rng, bkg_rms=bkg_rms diff --git a/tests/module/test_ngmix_weight_validation.py b/tests/module/test_ngmix_weight_validation.py index 2dad3a8c..a09a8106 100644 --- a/tests/module/test_ngmix_weight_validation.py +++ b/tests/module/test_ngmix_weight_validation.py @@ -41,6 +41,7 @@ import pytest from shapepipe.modules.ngmix_package.ngmix import ( + PSF_NOISE, Postage_stamp, do_ngmix_metacal, get_prior, @@ -318,3 +319,112 @@ def test_metacal_ivar_beats_binary_under_noise_gradient( f"metacal ivar/binary shape-scatter ratio {ratio:.3f} >= 0.9: the" " rms-weight advantage does not survive the metacal/fixnoise path" ) + + +# --------------------------------------------------------------------------- +# PSF-observation weight (#749 / #774) +# +# make_ngmix_observation built the PSF observation weightless, so ngmix +# defaulted to unit weight on the unit-flux PSF stamp and the galaxy +# GPriorBA(0.4) prior shared with the PSF fitter swamped the likelihood. The +# recovered PSF *shape* then collapsed toward the prior (zero ellipticity). +# With the column grammar (#761) this corrupts the *_psf_orig diagnostic +# columns, fit from the original image PSF by average_original_psf. A finite +# flat weight = 1/PSF_NOISE**2 restores the likelihood. The fix touches only +# the diagnostic original-PSF fit; metacal builds its own high-weight +# reconvolution-kernel PSF observation, so the calibrated shear is unchanged. +# --------------------------------------------------------------------------- + +PSF_SHEAR = (0.05, -0.03) # injected PSF ellipticity for the recovery test + + +def _psf_orig_via_metacal(psf_noise, psf_shear=PSF_SHEAR, seed=7): + """Run the real production path and return (HSM truth, psf_orig, resdict). + + Builds a sheared-PSF stamp with the module's own ``make_data`` simulator, + pushes it through ``do_ngmix_metacal`` (which fits the original image PSF + via ``average_original_psf`` into the ``*_psf_orig`` columns), under a + given ``PSF_NOISE``. ``psf_noise=1.0`` reproduces the pre-fix unit weight. + """ + import shapepipe.modules.ngmix_package.ngmix as ngm + from shapepipe.testing.simulate import make_data + + old = ngm.PSF_NOISE + ngm.PSF_NOISE = psf_noise + try: + rng = np.random.RandomState(seed) + prior = get_prior(PIX_SCALE, rng) + gals, psfs, _, weights, flags, jacobs = make_data( + rng=np.random.RandomState(123), shear=(0.0, 0.0), noise=1e-5, + n_epochs=2, img_size=51, psf_shear=psf_shear, + ) + hsm = galsim.hsm.FindAdaptiveMom(galsim.Image(psfs[0], scale=PIX_SCALE)) + g_truth = np.array([hsm.observed_shape.g1, hsm.observed_shape.g2]) + stamp = Postage_stamp(bkg_sub=False, megacam_flip=False) + (stamp.gals, stamp.psfs, stamp.weights, stamp.flags, stamp.jacobs) = ( + gals, psfs, weights, flags, jacobs, + ) + resdict, _psf_reconv, psf_orig = do_ngmix_metacal( + stamp, prior, 1.0, rng, + ) + return g_truth, psf_orig, resdict + finally: + ngm.PSF_NOISE = old + + +def test_make_ngmix_observation_psf_obs_is_weighted(): + """Structural guard: psf_obs carries the finite 1/PSF_NOISE**2 weight. + + The whole #749 fix is this weight; a regression to the weightless build + (unit weight) is what we must never ship again. + """ + from shapepipe.testing.simulate import make_data + + gals, psfs, _, weights, flags, jacobs = make_data( + rng=np.random.RandomState(123), shear=(0.0, 0.0), noise=1e-5, + n_epochs=1, img_size=51, psf_shear=PSF_SHEAR, + ) + obs = make_ngmix_observation( + gals[0], weights[0], flags[0], psfs[0], jacobs[0], + np.random.RandomState(0), + ) + npt.assert_allclose(obs.psf.weight, 1.0 / PSF_NOISE ** 2) + assert (obs.psf.weight > 1.0).all(), "PSF obs fell back to unit weight" + + +def test_psf_weight_recovers_original_psf_ellipticity(): + """#749 done-condition, end to end: *_psf_orig recovers the PSF shape. + + Through the real ``do_ngmix_metacal`` / ``average_original_psf`` path, the + original-PSF ellipticity matches the injected PSF moment shape with the + weight, and collapses toward the prior (zero) without it — the contrast + that proves the weight is load-bearing for the diagnostic columns. + """ + g_truth, psf_orig_w, _ = _psf_orig_via_metacal(PSF_NOISE) + _, psf_orig_u, _ = _psf_orig_via_metacal(1.0) # pre-fix unit weight + + # The fix: original-PSF column recovers the PSF moment shape. + npt.assert_allclose(psf_orig_w["g_psf"], g_truth, atol=1e-3) + + # The bug it cures: weightless collapses the shape toward the prior. + assert np.linalg.norm(psf_orig_u["g_psf"]) < 0.1 * np.linalg.norm(g_truth), ( + "weightless original-PSF fit did not collapse -- the #749 contrast is" + " gone, so this test no longer proves the weight is load-bearing" + ) + + +def test_psf_weight_leaves_metacal_shear_invariant(): + """The fix must not move the calibration: metacal shear flat in PSF_NOISE. + + Metacal rebuilds the reconvolution-kernel PSF observation at its own high + weight, independent of the diagnostic psf_obs, so every metacal-type shear + is invariant to the diagnostic weight. Asserted bit-for-bit between the + unit weight (pre-fix) and the shipped PSF_NOISE. + """ + _, _, res_w = _psf_orig_via_metacal(PSF_NOISE) + _, _, res_u = _psf_orig_via_metacal(1.0) + for mtype in ("noshear", "1p", "1m", "2p", "2m"): + npt.assert_allclose( + res_w[mtype]["g"], res_u[mtype]["g"], rtol=0, atol=1e-6, + err_msg=f"PSF-obs weight moved the metacal {mtype} shear", + ) From 6d30fc14460cf5e595a9d51d09a1f1c0288f1b7a Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Wed, 1 Jul 2026 00:16:04 +0200 Subject: [PATCH 2/2] tests: guard the PSF-obs weight end-to-end (#749) Three assertions in test_ngmix_weight_validation.py, on the real production path now that this sits on the column grammar (#761): - structural: make_ngmix_observation builds psf_obs with the finite 1/PSF_NOISE**2 weight, never unit weight; - recovery (#749 done-condition): through do_ngmix_metacal / average_original_psf, the *_psf_orig column recovers the injected PSF moment ellipticity to 1e-3 with the weight, and collapses below 10% without it (the contrast that proves the weight is load-bearing for the diagnostic columns); - invariance: every metacal-type shear is bit-for-bit identical between the unit-weight (pre-fix) and shipped PSF_NOISE builds (metacal's prior-free AdmomFitter is insensitive to the flat weight's scale). Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_017KQn6mXoL4XMgscZDsaTkB --- tests/module/test_ngmix_weight_validation.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/tests/module/test_ngmix_weight_validation.py b/tests/module/test_ngmix_weight_validation.py index a09a8106..abc908fd 100644 --- a/tests/module/test_ngmix_weight_validation.py +++ b/tests/module/test_ngmix_weight_validation.py @@ -331,8 +331,9 @@ def test_metacal_ivar_beats_binary_under_noise_gradient( # With the column grammar (#761) this corrupts the *_psf_orig diagnostic # columns, fit from the original image PSF by average_original_psf. A finite # flat weight = 1/PSF_NOISE**2 restores the likelihood. The fix touches only -# the diagnostic original-PSF fit; metacal builds its own high-weight -# reconvolution-kernel PSF observation, so the calibrated shear is unchanged. +# the diagnostic original-PSF fit (run with the galaxy prior); metacal fits the +# same psf_obs with a prior-free AdmomFitter, which is insensitive to a flat +# weight's scale, so the calibrated shear is unchanged. # --------------------------------------------------------------------------- PSF_SHEAR = (0.05, -0.03) # injected PSF ellipticity for the recovery test @@ -416,10 +417,11 @@ def test_psf_weight_recovers_original_psf_ellipticity(): def test_psf_weight_leaves_metacal_shear_invariant(): """The fix must not move the calibration: metacal shear flat in PSF_NOISE. - Metacal rebuilds the reconvolution-kernel PSF observation at its own high - weight, independent of the diagnostic psf_obs, so every metacal-type shear - is invariant to the diagnostic weight. Asserted bit-for-bit between the - unit weight (pre-fix) and the shipped PSF_NOISE. + Metacal fits the same psf_obs with a prior-free AdmomFitter (then dilates + to a round reconvolution kernel); a prior-free moment fit is insensitive to + a flat weight's absolute scale, so every metacal-type shear is invariant to + the diagnostic weight. Asserted bit-for-bit between the unit weight + (pre-fix) and the shipped PSF_NOISE. """ _, _, res_w = _psf_orig_via_metacal(PSF_NOISE) _, _, res_u = _psf_orig_via_metacal(1.0)