[BUG] ngmix module passes the galaxy shape prior to the PSF fit, which is prior-dominated; g_psfo/Tpsf columns record the metacal reconvolution PSF, not the original PSF

[cailmdaley]

This writes up the "PSF fit is prior-dominated" finding from PR #741 (comment, 06-11) as its own thread, with the reproduction code @fabianhervaspeters asked for — and resolves the prior-vs-fitgauss question, including @fabianhervaspeters's full-pipeline result that removing the prior changed nothing. One claim from the original comment is also corrected below.

What the code does

Two distinct PSF "fits" exist in the module, and the discussion conflated them:

1. The PSF-model fit — make_runners builds psf_fitter = ngmix.fitting.Fitter(model='gauss', prior=prior) (ngmix.py:1138) where prior is the galaxy joint prior from get_prior, containing GPriorBA(sigma=0.4) (ngmix.py:79). The same construction exists pre-June (bd60dc8e, line 1061). The PSF Observation is built with no weight map (ngmix.py:1029; ngmix defaults to unit weights) on a unit-flux stamp.
2. The metacal reconvolution kernel — metacal_pars = {'psf': 'fitgauss'} (ngmix.py:1193). ngmix's MetacalFitGaussPSF does a prior-free admom fit and then builds a round galsim.Gaussian, dilated. Isotropic by construction — and by design; that part is a feature.

What we measured

Setup (deterministic, N=30, seeds fixed; container shapepipe_ngmix_v2.0-dev.sif, ngmix 2.4.0): CFIS-like Gaussian PSF, fwhm 0.65″, sheared to g=(0.025, −0.015), drawn at 0.187″/px on a 48×48 unit-flux stamp. Pixel-window dilution makes the stamp's actual ellipticity g=(0.02406, −0.01444) — confirmed independently by HSM FindAdaptiveMom.

Fit of the original PSF stamp	g1	g2	T
stamp truth (HSM)	0.02406	−0.01444	—
module fitter, with galaxy prior (plain Bootstrapper, no metacal anywhere in the path)	9.1×10⁻⁵ ± 1.1×10⁻⁵	−3.7×10⁻⁵ ± 1.6×10⁻⁵	0.1581 ± 0.0013
identical fitter, prior removed	0.024058	−0.014442	0.15854

The shape is suppressed ~265×; the size is fine (the T prior is flat — only g, where GPriorBA is informative, is destroyed).

Why the prior wins: on a unit-flux, unit-weight stamp, forcing the best-fit Gaussian round costs only Δχ²/2 ≈ 1.4×10⁻⁵, while GPriorBA(0.4) rewards it by Δln p ≈ 4.0×10⁻³ — a ~290× pull ratio, which quantitatively predicts the observed suppression (1/(1+290) ≈ 3.4×10⁻³ vs measured 3.8×10⁻³). With the stamp at the star's flux (2800) the likelihood would win by ~10⁴ — the unit-flux/no-weight-map construction is exactly what makes the fit prior-dominated.

Resolving "prior-dominated" vs "fitgauss is isotropic by construction"

Both are true — about different fits:

• The g1 = 9×10⁻⁵ result comes from a plain ngmix.bootstrap.Bootstrapper fit of the original PSF stamp; fitgauss appears nowhere in that code path, so it cannot be the explanation there.
• The discriminating test is the prior removal: a structurally isotropic model could never return g=(0.0241, −0.0144); the same Fitter(model='gauss') without the prior does, matching HSM to 5 decimals. The model has free shape; the galaxy prior is what rounds it.
• fitgauss is isotropic by construction — but it is the reconvolution kernel, which is supposed to be round.

The two effects stack at the catalog level: metacal_bootstrap runs the module's psf_runner on the metacal observations, whose PSF is the already-round, dilated fitgauss kernel. average_multiepoch_psf then propagates that fit into g_psfo/g_psf/Tpsf — documented as "the original image psf from psfex or mccd" (ngmix.py:354). In the catalog layer of our run: g1_psfo_ngmix = −1.5×10⁻⁶ and Tpsf = 0.17407 vs 0.1581 for the actual (pixel-convolved) PSF — the size of the dilated reconvolution kernel, ~10% high. So the "PSF original" columns currently record neither the original PSF's shape nor its size.

Directly addressing Fabian's prior-removal result

@fabianhervaspeters — your full-pipeline test (removing the prior on psf_fitter and psf_guess, and g1/g2 still ≈0) is exactly what the two-fits picture predicts, and it's a good piece of evidence, not a puzzle. In the full pipeline, metacal_bootstrap runs the module's psf_runner on the metacal observations — whose PSF is the fitgauss reconvolution kernel, round by construction. So the thing your fitter sees is already isotropic: there's no ellipticity in it for a prior to suppress, and removing the prior can't recover one that isn't there. Your intuition was right for your path — half a round well earned.

The prior-domination only bites the other fit: a direct fit of the original elliptical PSF stamp (verify_psf_prior.py / the direct Bootstrapper path), which the full-pipeline catalog columns never perform. So the two observations are consistent: prior-removal is a no-op where you tested it (round kernel), and decisive where we tested it (original PSF). To confirm we're looking at the same code path: was metacal on in your run, and did the psf_runner see the original PSF observation or the metacal-generated one? If it was the metacal one, we're fully aligned.

One correction to the PR #741 comment

I wrote there that "the deconvolution kernel is rounder than the true PSF". That is wrong: ngmix deconvolves by the actual PSF image (galsim.Deconvolve of the interpolated stamp), not by the fitted model. Metacal shear itself is therefore plausibly unaffected by the prior-dominated fit. What is affected: the PSF catalog columns (hence any ρ-statistics / PSF-leakage diagnostics built from them), and any non-metacal use of the module's PSF fit, where the galaxy model is convolved with the prior-rounded PSF gmix.

Both SHAs

Run identically on current tip (module unchanged from c5a9cc4d through HEAD 5a51dd65) and on bd60dc8e (06-05): the PSF-fit numbers agree to the digit. Not introduced by the June fixes. Full reruns are byte-identical.

Suggested next steps

1. Stop passing the galaxy prior to the PSF fitter — either no prior (our prior-free run shows the fit is then accurate to ~10⁻⁵) or a PSF-appropriate non-informative one. One-line change in make_runners.
2. Decide what g_psfo/Tpsf should mean. If "original PSF", fit the pre-metacal PSF observation (and do it after fix 1, or the prior will silently re-round it); if "reconvolution kernel", rename and document. Right now the columns match neither their docstring nor the original PSF.
3. gauss vs fitgauss reconvolution comparison — @fabianhervaspeters reported a 20–30% improvement in star response with psf='gauss'; worth a controlled A/B on the same stamps, which the experiment script makes easy.
4. Optionally fold the point-source case into the validation harness (test_ngmix_weight_validation.py) with an assertion on recovered PSF ellipticity, so this can't regress silently.

Reproduction (candide)

Everything in /n17data/cdaley/unions/scratch-wf/star-experiment/:

• verify_psf_prior.py — the 19-line discriminating test (HSM + prior-free fit). apptainer exec --bind /n17data /n17data/cdaley/containers/shapepipe_ngmix_v2.0-dev.sif python verify_psf_prior.py
• star_experiment.py — the full two-path experiment (module metacal path + direct Bootstrapper path + catalog layer); feature-detects the module API so the same file runs on both SHAs via PYTHONPATH.
• tip_A.json / old_A.json — raw numbers quoted above (results.direct.psf_g1 etc.); rerun_*.json byte-identical replays.
• Old-SHA worktree: /n17data/cdaley/unions/scratch-wf/shapepipe-fabian-sha (bd60dc8e).

— Claude, on behalf of Cail.