Repository navigation
[BUG] ngmix module passes the galaxy shape prior to the PSF fit, which is prior-dominated #749
Description
Activity
Hello,
To follow up on the discussion in #741 I made a very simple example that shows that the prior has no impact when measuring the PSF. I don't see any code for the results presented above so I don't how you find those results.
Here are my findings:True psf shape: 0.02 -0.01 Galsim FindAdaptiveMome: 0.019099739248397904 -0.009599278364509202 NGMix psf_runner with galaxy prior: 0.019098474073373506 -0.00959589621223074
Here the python script to reproduce them:
import numpy as np import ngmix import galsim def get_prior(pixel_scale, rng, T_range=None, F_range=None): """Build ngmix joint prior for a 6-parameter galaxy model. Parameters ---------- pixel_scale : float Pixel scale in arcsec (sets centroid prior width). rng : numpy.random.RandomState Random state for all priors. T_range : tuple of float, optional (min, max) for flat size prior; default (-1, 1e3). F_range : tuple of float, optional (min, max) for flat flux prior; default (-100, 1e9). Returns ------- ngmix.joint_prior.PriorSimpleSep """ if T_range is None: T_range = [-1.0, 1.0e3] if F_range is None: F_range = [-100.0, 1.0e9] cen_prior = ngmix.priors.CenPrior( cen1=0.0, cen2=0.0, sigma1=pixel_scale, sigma2=pixel_scale, rng=rng, ) g_prior = ngmix.priors.GPriorBA(sigma=0.4, rng=rng) T_prior = ngmix.priors.FlatPrior( minval=T_range[0], maxval=T_range[1], rng=rng ) F_prior = ngmix.priors.FlatPrior( minval=F_range[0], maxval=F_range[1], rng=rng ) return ngmix.joint_prior.PriorSimpleSep( cen_prior=cen_prior, g_prior=g_prior, T_prior=T_prior, F_prior=F_prior, ) def make_runners(prior, flux_guess, rng): """Build the module's galaxy and PSF runners. Single source of truth for the fitter configuration (Gaussian galaxy and PSF models, guessers, retry counts), shared by the metacal bootstrap below and by validation tests that fit module-built observations directly. Parameters ---------- prior : ngmix.joint_prior.PriorSimpleSep Priors for the fitting parameters. flux_guess : float Initial flux guess. rng : numpy.random.RandomState Random state for the guessers. Returns ------- tuple (runner, psf_runner) : ngmix.runners.Runner, ngmix.runners.PSFRunner """ fitter = ngmix.fitting.Fitter(model="gauss", prior=prior) guesser = ngmix.guessers.TPSFFluxAndPriorGuesser( rng=rng, T=0.25, prior=prior ) psf_fitter = ngmix.fitting.Fitter(model="gauss", prior=prior) psf_guesser = ngmix.guessers.TFluxGuesser( rng=rng, T=0.25, prior=prior, flux=flux_guess ) return ( ngmix.runners.Runner(fitter=fitter, guesser=guesser, ntry=5), ngmix.runners.PSFRunner( fitter=psf_fitter, guesser=psf_guesser, ntry=2 ), ) def main(): rng = np.random.RandomState(seed=42) img_size = 51 g1_true = 0.02 g2_true = -0.01 print("True psf shape:", g1_true, g2_true) psf = ( galsim.Moffat(fwhm=0.65, beta=2.5) .shear(g1=g1_true, g2=g2_true) .withFlux(1.0) ) gal_img = psf.drawImage(nx=img_size, ny=img_size, scale=0.2) img = gal_img.array img += rng.normal(size=img.shape) * 1e-5 wgt = np.ones_like(img) / (1e-5**2) jac = ngmix.DiagonalJacobian( row=(img_size - 1) / 2, col=(img_size - 1) / 2, scale=0.2, ) obs = ngmix.Observation( image=img, weight=wgt, jacobian=jac, ) s = galsim.hsm.FindAdaptiveMom(gal_img) print("Galsim FindAdaptiveMome:", s.observed_shape.g1, s.observed_shape.g2) rng_shape = np.random.RandomState(seed=42) prior = get_prior(pixel_scale=0.2, rng=rng_shape) gal_runner, psf_runner = make_runners( prior=prior, flux_guess=1.0, rng=rng_shape ) res_psf_runner = psf_runner.go(obs) print( "NGMix psf_runner with galaxy prior:", res_psf_runner["g"][0], res_psf_runner["g"][1], ) if __name__ == "__main__": main()
- added a parent issue
on Jun 20, 2026 Thank you for looking into this @aguinot. I was able to reproduce your results, the difference is how the weight map is defined and used. In your example, you do:
wgt = np.ones_like(img) / (1e-5**2) obs = ngmix.Observation( image=img, weight=wgt, jacobian=jac, )
whereas in the library code (both on the previous develop state and the ngmix branch), no weight map is used for the PSF fit, so I think ngmix defaults to unit weights:
shapepipe/src/shapepipe/modules/ngmix_package/ngmix.py
Lines 1012 to 1022 in fc656da
If I take your example and drop the
weight=wgtargument, I get:True PSF shape: g1=+0.02000 g2=-0.01000 (A) Axel's weight=wgt (~1e10): g1=+0.019098 g2=-0.009596 (B) weight= line removed (unit): g1=+0.000045 g2=-0.000027in addition, if I instead unset the prior e.g.
psf_fitter = ngmix.fitting.Fitter(model="gauss", prior=None), I get(C) weight=weight= line removed (unit), prior=None: g1=+0.019098 g2=-0.009596So my understanding is that there are two ways to recover the input shear: pass a proper weight map, or don't apply the galaxy prior to the PSF fit. The latter seems more straightforward to me. Tell me what you think!
- changed the title
[-][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[/-][+][BUG] ngmix module passes the galaxy shape prior to the PSF fit, which is prior-dominated[/+]on Jun 21, 2026 Interesting! Is it the case that either option of adding the weight or applying a prior regularises the problem such that the input shear is recovered? It sounds like using the weight and removing the prior makes more sense.
Not quite, either adding weights (B above) or dropping the galaxy prior (C above) addresses the problem. With a unit weight map, which is the default if a weight map is not supplied, the PSF fit is dominated by the galaxy prior. I think the simplest solution is to not use a prior, I don't understand why one would use the same prior for the galaxy and PSF in the first place.
A weight should always be provided when doing model fitting, this was a mistake on my part not to include it in the original code. Actually it is also common to add noise on the PSF to stabilize the fit (like gaussian noise with sigma = 1e-5 or 1e-6). Given that the PSF models from PSFEx already contain a little bit of noise just adding a weight should be fine.
The prior is not ideal to fit the PSF but it is not completely wrong either. I would put that as a lower priority given that it still works and simplify things. Otherwise you can define a better PSF fitter like they do in metadetect here.OK it sounds that adding the weight will solve the bug, let´s implement it!
- added 4 commits that reference this issue
on Jun 30, 2026 - linked a pull request that will close this issuengmix: give the PSF observation a finite weight map (#749) #779
on Jul 1, 2026 - added a commit that references this issue
on Jul 1, 2026 Closed by #779, thanks all!
- added 6 commits that reference this issue
on Jul 21, 2026
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-
fitgaussquestion, 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:
make_runnersbuildspsf_fitter = ngmix.fitting.Fitter(model='gauss', prior=prior)(ngmix.py:1138) whereprioris the galaxy joint prior fromget_prior, containingGPriorBA(sigma=0.4)(ngmix.py:79). The same construction exists pre-June (bd60dc8e, line 1061). The PSFObservationis built with no weight map (ngmix.py:1029; ngmix defaults to unit weights) on a unit-flux stamp.metacal_pars = {'psf': 'fitgauss'}(ngmix.py:1193). ngmix'sMetacalFitGaussPSFdoes a prior-free admom fit and then builds a roundgalsim.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 HSMFindAdaptiveMom.Bootstrapper, no metacal anywhere in the path)The shape is suppressed ~265×; the size is fine (the T prior is flat — only g, where
GPriorBAis 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:
ngmix.bootstrap.Bootstrapperfit of the original PSF stamp;fitgaussappears nowhere in that code path, so it cannot be the explanation there.Fitter(model='gauss')without the prior does, matching HSM to 5 decimals. The model has free shape; the galaxy prior is what rounds it.fitgaussis isotropic by construction — but it is the reconvolution kernel, which is supposed to be round.The two effects stack at the catalog level:
metacal_bootstrapruns the module'spsf_runneron the metacal observations, whose PSF is the already-round, dilatedfitgausskernel.average_multiepoch_psfthen propagates that fit intog_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⁻⁶ andTpsf= 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_fitterandpsf_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_bootstrapruns the module'spsf_runneron the metacal observations — whose PSF is thefitgaussreconvolution 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 directBootstrapperpath), 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 thepsf_runnersee 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.Deconvolveof 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
c5a9cc4dthrough HEAD5a51dd65) and onbd60dc8e(06-05): the PSF-fit numbers agree to the digit. Not introduced by the June fixes. Full reruns are byte-identical.Suggested next steps
make_runners.g_psfo/Tpsfshould 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.gaussvsfitgaussreconvolution comparison — @fabianhervaspeters reported a 20–30% improvement in star response withpsf='gauss'; worth a controlled A/B on the same stamps, which the experiment script makes easy.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.pystar_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 viaPYTHONPATH.tip_A.json/old_A.json— raw numbers quoted above (results.direct.psf_g1etc.);rerun_*.jsonbyte-identical replays./n17data/cdaley/unions/scratch-wf/shapepipe-fabian-sha(bd60dc8e).— Claude, on behalf of Cail.