review: address PR #761 code review

- do_ngmix_metacal: seed the original-PSF pre-fit from a snapshot of the
  metacal RNG so it no longer advances the rng boot.go consumes — galaxy/shear
  results are now bit-identical to the pre-PR base (only the *_PSF_ORIG columns
  are added). Conservative: galaxy science unchanged vs base.
- remove the redundant galaxy.sigma_to_fwhm wrapper; inline
  cs_size.sigma_to_fwhm(sigma) * pixel_scale at call sites (spread_model);
  cs_util unchanged.
- drop leftover psf_fit_prior knob references (docstrings/comments/tests).
- canonicalize the ESTIMATOR_COMPONENT_OBJECT grammar spec string across all sites.
- trim over-explained copy-PSF / leakage / weighting prose to one canonical spot;
  document that bit-identity rests on BOTH pristine gal_obs.psf AND the snapshot RNG.
- add tests: PSF columns survive a failed-fit NaN-fill; all-epochs-failed raises;
  grammar token check runs over both shipped param files. Fix a stale cross-ref.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

Author: cailmdaley

src/shapepipe/modules/make_cat_package/make_cat.py

@@ -329,18 +329,16 @@ def _save_ngmix_data(self, ngmix_cat_path, moments=False):
 
         Save the NGMIX catalogue into the final one.
 
-        Column grammar: ``NGMIX[m]_<COMPONENT>[_ERR][_<OBJECT>][_<SHEAR>]``.
-        ``OBJECT`` and ``SHEAR`` are both optional: the per-object metadata
-        columns (``NGMIX_MCAL_FLAGS``, ``NGMIX_N_EPOCH``,
-        ``NGMIX_MCAL_TYPES_FAIL``) carry neither. When present, ``OBJECT`` is
-        one of ``GAL`` (galaxy), ``PSF_ORIG`` (the original psfex/mccd image
-        PSF, fit by ``ngmix.average_original_psf`` — its true ellipticity and
-        size, feeding object-wise PSF leakage), or ``PSF_RECONV`` (the metacal
-        reconvolution kernel, round and enlarged by construction, fit by
-        ``ngmix.average_multiepoch_psf`` — the Tgal / Tpsf size cut and a g~0
-        sanity check). ``PSF_ORIG`` and ``PSF_RECONV`` are independent fits of
-        different PSFs, no longer the single aliased value of the pre-fix code
-        (shapepipe#749).
+        Column grammar: ``NGMIX[m]_<COMPONENT>[_ERR]_<OBJECT>_<SHEAR>``, plus
+        three OBJECT/SHEAR-less per-object metadata columns
+        (``NGMIX[m]_MCAL_FLAGS``, ``NGMIX_N_EPOCH``,
+        ``NGMIX_MCAL_TYPES_FAIL``). ``OBJECT`` is one of ``GAL`` (galaxy),
+        ``PSF_ORIG`` (the original image PSF, fit by
+        :func:`ngmix.average_original_psf`), or ``PSF_RECONV`` (the metacal
+        reconvolution kernel, fit by :func:`ngmix.average_multiepoch_psf`);
+        see those functions for what each PSF family IS. ``PSF_ORIG`` and
+        ``PSF_RECONV`` are independent fits of different PSFs, no longer the
+        single aliased value of the pre-fix code (shapepipe#749).
 
         Parameters
         ----------
@@ -390,9 +388,11 @@ def _save_ngmix_data(self, ngmix_cat_path, moments=False):
         prefix = f"NGMIX{m}"
 
         # Galaxy (GAL), original image PSF (PSF_ORIG) and metacal
-        # reconvolution kernel (PSF_RECONV). G1/G2 are scalar reduced-shear
-        # components, not a 2-vector. Sentinels: sizes/fluxes/mags/flags 0,
-        # *_ERR fluxes/mags -1, ellipticities -10, *_ERR sizes 1e30.
+        # reconvolution kernel (PSF_RECONV); see ngmix.average_original_psf /
+        # average_multiepoch_psf for what each PSF family is. G1/G2 are scalar
+        # reduced-shear components, not a 2-vector. Sentinels:
+        # sizes/fluxes/mags/flags 0, *_ERR fluxes/mags -1, ellipticities -10,
+        # *_ERR sizes 1e30.
         for key_str in (
             f"{prefix}_T_GAL_",
             f"{prefix}_SNR_GAL_",
@@ -484,8 +484,7 @@ def _save_ngmix_data(self, ngmix_cat_path, moments=False):
                         ncf_data["flags"][ind[0]], idx
                     )
 
-                    # Original image PSF (psfex/mccd) — its genuine,
-                    # un-reconvolved ellipticity and size.
+                    # Original image PSF (see ngmix.average_original_psf).
                     self._add2dict(
                         f"{prefix}_G1_PSF_ORIG_{key}",
                         ncf_data["g1_psf_orig"][ind[0]], idx
@@ -511,8 +510,8 @@ def _save_ngmix_data(self, ngmix_cat_path, moments=False):
                         ncf_data["T_err_psf_orig"][ind[0]], idx
                     )
 
-                    # Metacal reconvolution kernel — round, enlarged by
-                    # construction; the Tgal/Tpsf cut and g~0 sanity check.
+                    # Metacal reconvolution kernel (see
+                    # ngmix.average_multiepoch_psf).
                     self._add2dict(
                         f"{prefix}_G1_PSF_RECONV_{key}",
                         ncf_data["g1_psf_reconv"][ind[0]], idx

src/shapepipe/modules/ngmix_package/ngmix.py

@@ -380,16 +380,13 @@ def compile_results(self, results):
             Compiled results ready to be written to a file.
 
             Two PSF column families — each carrying ellipticity *and* size,
-            for *different* PSFs (shapepipe#749):
+            for *different* PSFs (shapepipe#749). See the function docstrings
+            for what each family IS:
 
-            * ``*_psf_orig`` (``g1``/``g2`` + ``*_err``, ``T``) — the
-              ORIGINAL image PSF (the psfex/mccd model stamp), fit before
-              metacal reconvolution by :func:`average_original_psf`. This is
-              the PSF whose true ellipticity and size enter object-wise
-              PSF-leakage diagnostics.
-            * ``*_psf_reconv`` — the metacal RECONVOLUTION kernel (round and
-              enlarged by construction, used for the Tgal/Tpsf size cut and a
-              g~0 sanity check), fit by :func:`average_multiepoch_psf`.
+            * ``*_psf_orig`` (``g1``/``g2`` + ``*_err``, ``T``) — the original
+              image PSF, fit by :func:`average_original_psf`.
+            * ``*_psf_reconv`` — the metacal reconvolution kernel, fit by
+              :func:`average_multiepoch_psf`.
 
         Raises
         ------
@@ -1232,23 +1229,18 @@ def average_original_psf(gal_obs_list, psf_runner):
     """Fit and average the *original* image PSF over epochs.
 
     The original PSF is the psfex/mccd model stamp handed to ngmix
-    (``gal_obs.psf``), fit here with the same ``psf_runner`` machinery — and
-    so the same ``psf_fit_prior`` — used inside metacal, but on the PSF
-    *before* metacal's reconvolution. Exported to the original-PSF columns
+    (``gal_obs.psf``), fit here with the same ``psf_runner`` (hence the same
+    fit prior and guesser) used inside metacal, but on the PSF *before*
+    metacal's reconvolution. Exported to the original-PSF columns
     (``NGMIX_G1/G2_PSF_ORIG``, ``NGMIX_T_PSF_ORIG``). Distinct from the
     reconvolution-kernel fit (:func:`average_multiepoch_psf`): the original
     PSF retains its true ellipticity and size, whereas the reconvolution
     kernel is round and enlarged by construction. This is the PSF whose true
     shape and size enter object-wise PSF-leakage diagnostics.
 
-    Epochs are weighted by the *galaxy* inverse-variance weight
-    (``gal_obs.weight.sum()``). This is the same *form* as
-    :func:`average_multiepoch_psf` (``weight.sum()`` per epoch, skipping
-    ``flags != 0`` fits) but not the identical weight: the reconvolution path
-    weights by the fixnoise-combined inverse variance of the noshear metacal
-    image, whereas this path uses the raw galaxy inverse variance. So the two
-    PSF families share an averaging *scheme* and differ in which PSF is fit
-    and in the precise per-epoch weighting factor.
+    Weighted by the raw galaxy inverse variance (``gal_obs.weight.sum()`` per
+    epoch); :func:`average_multiepoch_psf` uses the same scheme but the
+    fixnoise-combined metacal-image weight instead.
 
     The fit runs on a *copy* of each PSF observation so ``gal_obs.psf`` —
     the object metacal later deep-copies and consumes via
@@ -1258,8 +1250,11 @@ def average_original_psf(gal_obs_list, psf_runner):
     would survive metacal's deep copy and be reused as the
     ``MetacalFitGaussPSF`` fallback when its own admom+ML PSF fits both fail,
     silently rescuing objects the base branch dropped (``BootPSFFailure``) and
-    changing the galaxy/shear result set. Fitting a copy keeps this add-column
-    refactor bit-identical on the galaxy results.
+    changing the galaxy/shear result set. Fitting a copy closes this
+    PSF-aliasing channel; combined with :func:`do_ngmix_metacal` seeding this
+    pre-fit from a *snapshot* of the metacal RNG (so the pre-fit does not
+    advance it), the add-column refactor stays bit-identical on the galaxy
+    results.
 
     Parameters
     ----------
@@ -1268,19 +1263,17 @@ def average_original_psf(gal_obs_list, psf_runner):
         (pre-metacal) PSF observation to fit, with no further ``.psf`` of
         its own so the runner fits the stamp itself. Left pristine.
     psf_runner : ngmix.runners.PSFRunner
-        The module's PSF runner, carrying the resolved ``psf_fit_prior``.
+        The module's PSF runner (built by :func:`make_runners` from the
+        shared ``prior``).
 
     Returns
     -------
     dict
         Same keys as :func:`average_multiepoch_psf`.
     """
     def fit(gal_obs):
-        # Fit a COPY of the PSF observation: PSFRunner.go fits the PSF stamp
-        # and sets its .meta['result'] (and .gmix on success). Reading from
-        # the copy leaves gal_obs.psf pristine — no gmix to leak through
-        # metacal's deep copy into the MetacalFitGaussPSF fallback. Failed
-        # fits keep flags != 0 and are dropped by _average_psf_fits.
+        # Fit a COPY so gal_obs.psf stays pristine for metacal — see docstring.
+        # Failed fits keep flags != 0 and are dropped by _average_psf_fits.
         psf_obs = gal_obs.psf.copy()
         psf_runner.go(psf_obs)
         return psf_obs.meta['result'], gal_obs.weight.sum()
@@ -1379,14 +1372,15 @@ def do_ngmix_metacal(stamp, prior, flux_guess, rng, centroid_source="hsm"):
 
     runner, psf_runner = make_runners(prior, flux_guess, rng)
 
-    # Fit the ORIGINAL (psfex/mccd) PSF before metacal reconvolves it, using
-    # the same psf_runner (so the same psf_fit_prior and centroid). This is
-    # the PSF_ORIG family; metacal below produces the PSF_RECONV family. Run
-    # first so the original PSF is fit on its own stamp, distinct from the
-    # round, enlarged kernel metacal convolves back in. average_original_psf
-    # fits a COPY of each gal_obs.psf, so gal_obs_list reaches boot.go below
-    # pristine — no stray gmix to leak into MetacalFitGaussPSF's fallback.
-    psf_orig_res = average_original_psf(gal_obs_list, psf_runner)
+    # Fit the ORIGINAL (psfex/mccd) PSF before metacal reconvolves it. Use a
+    # psf_runner seeded from a snapshot of rng's state so this prefit does NOT
+    # advance the rng that boot.go consumes below — keeping the galaxy/shear
+    # results bit-identical to the no-prefit branch. average_original_psf fits a
+    # COPY of each gal_obs.psf, so gal_obs_list also reaches boot.go pristine.
+    prefit_rng = np.random.RandomState()
+    prefit_rng.set_state(rng.get_state())
+    _, prefit_psf_runner = make_runners(prior, flux_guess, prefit_rng)
+    psf_orig_res = average_original_psf(gal_obs_list, prefit_psf_runner)
 
     metacal_pars = {
         'types': ['noshear', '1p', '1m', '2p', '2m'],

src/shapepipe/modules/spread_model_package/spread_model.py

@@ -8,10 +8,10 @@
 
 import galsim
 import numpy as np
+from cs_util import size as cs_size
 from sqlitedict import SqliteDict
 
 from shapepipe.pipeline import file_io
-from shapepipe.utilities import galaxy
 
 
 def get_sm(obj_vign, psf_vign, model_vign, weight_vign):
@@ -102,7 +102,7 @@ def get_model(sigma, flux, img_shape, pixel_scale=0.186):
 
     """
     # Get scale radius
-    scale_radius = 1 / 16 * galaxy.sigma_to_fwhm(sigma, pixel_scale=pixel_scale)
+    scale_radius = 1 / 16 * cs_size.sigma_to_fwhm(sigma) * pixel_scale
 
     # Get galaxy model
     gal_obj = galsim.Exponential(scale_radius=scale_radius, flux=flux)

src/shapepipe/utilities/galaxy.py

@@ -5,94 +5,3 @@
 :Author: Martin Kilbinger <martin.kilbinger@cea.fr>
 
 """
-
-import numpy as np
-from cs_util import size as cs_size
-
-
-def sigma_to_fwhm(sigma, pixel_scale=1.0):
-    r"""Convert Sigma to FWHM.
-
-    Transform standard deviation of a 1D Gaussian, sigma, to FWHM
-    (Full Width Half Maximum).
-
-    Thin wrapper over :func:`cs_util.size.sigma_to_fwhm`, the single
-    source of truth for the bare ``FWHM = 2 sqrt(2 ln 2) sigma``
-    conversion; this layer adds the optional ``pixel_scale`` rescaling
-    and ShapePipe's input validation.
-
-    Parameters
-    ----------
-    sigma : numpy.ndarray
-        input standard deviation(s)
-    pixel_scale : float, optional, default=1
-        pixel size in arcsec, set to 1 if no scaling
-        required
-
-    Returns
-    -------
-    fwhm : (array of) float
-        output fwhm(s)
-
-    Raises
-    ------
-    TypeError
-        If ``sigma`` is not of type numpy array or float
-    TypeError
-        If ``sigma`` array values are not of type float
-    TypeError
-        If ``pixel_scale`` is not of type float
-    ValueError
-        If ``sigma`` array values are not greater than 0.0
-    ValueError
-        If ``sigma`` is not greater than 0.0
-    ValueError
-        If ``pixel_scale`` is not greater than 0.0
-
-    Notes
-    -----
-    To compute the FWHMh for a 1D Gaussian N(x), solve the equation
-
-    .. math::
-
-        N(x) = (\sigma \sqrt{2\pi})^{-1} \exp[x^2/2\sigma^2] = \frac 1 2 N(x)
-
-
-    for :math:`x`. The FWHM is :math:`x + (-x) = 2x`. The solution is
-
-    .. math::
-
-        \textrm{FWHM} = 2 \sqrt(2 \ln 2) \sigma \approx 2.355 \sigma
-
-    """
-    if not isinstance(sigma, (np.ndarray, float)):
-        raise TypeError(
-            f"Sigma must be of type numpy array or float, not {type(sigma)}."
-        )
-    elif isinstance(sigma, np.ndarray) and sigma.dtype != np.float64:
-        raise TypeError(
-            f"Sigma array values must be of type float, not {sigma.dtype}."
-        )
-
-    if not isinstance(pixel_scale, float):
-        raise TypeError(
-            f"The pixel scale must of type float, not {type(pixel_scale)}."
-        )
-
-    if isinstance(sigma, np.ndarray) and np.any(sigma <= 0.0):
-        raise ValueError(
-            f"Found {sigma[sigma <=0].size} invalid standard deviation array "
-            + "values, all elements must to be greater than 0.0."
-        )
-    elif isinstance(sigma, float) and sigma <= 0.0:
-        raise ValueError(
-            f"Invalid standard deviation {sigma}, needs to be greater than "
-            + "0.0."
-        )
-
-    if pixel_scale <= 0.0:
-        raise ValueError(
-            f"Invalid pixel scale {pixel_scale}, needs to be greater than 0.0."
-        )
-
-    return cs_size.sigma_to_fwhm(sigma) * pixel_scale

tests/module/test_make_cat.py

@@ -2,10 +2,11 @@
 
 Drives ``SaveCatalogue._save_ngmix_data`` against a synthetic ngmix catalogue
 to lock in the renamed PSF-column grammar
-``NGMIX_<COMPONENT>[_ERR]_<OBJECT>[_<SHEAR>]`` (shapepipe#749): the original
-image PSF (``PSF_ORIG``) and the metacal reconvolution kernel (``PSF_RECONV``)
-are independent fits of *different* PSFs, no longer the single aliased value of
-the pre-fix code.
+``NGMIX[m]_<COMPONENT>[_ERR]_<OBJECT>_<SHEAR>`` (shapepipe#749), plus the three
+OBJECT/SHEAR-less metadata columns (``NGMIX[m]_MCAL_FLAGS``, ``NGMIX_N_EPOCH``,
+``NGMIX_MCAL_TYPES_FAIL``): the original image PSF (``PSF_ORIG``) and the metacal
+reconvolution kernel (``PSF_RECONV``) are independent fits of *different* PSFs,
+no longer the single aliased value of the pre-fix code.
 """
 
 import numpy as np

tests/module/test_ngmix.py

@@ -238,7 +238,7 @@ def test_compile_results_psf_families_are_unaliased():
     came from one ``average_multiepoch_psf`` result and so were byte-
     identical; here they must differ, each tracing its own source. The
     companion size un-aliasing is checked in
-    ``test_compile_results_size_columns_are_half_light_radii``.
+    ``test_compile_results_size_columns_are_unaliased``.
     """
     from shapepipe.modules.ngmix_package.ngmix import Ngmix
 
@@ -296,6 +296,14 @@ def test_compile_results_nan_fills_failed_fit_types():
     assert failed["flags"] == [0x8]
     assert failed["mcal_flags"] == [0x8] and failed["mcal_flags"][0] != 0
 
+    # The object-level PSF columns survive on a failed fit type: they are copied
+    # outside the flags==0 guard, so the failed type keeps a full-length PSF row
+    # (gating that copy on the galaxy flags would shorten it -> save crash).
+    npt.assert_allclose(failed["T_psf_orig"], [ORIG_PSF_T])
+    npt.assert_allclose(failed["T_psf_reconv"], [0.09])
+    npt.assert_allclose(failed["g1_psf_orig"], [ORIG_PSF_G[0]])
+    npt.assert_allclose(failed["g1_psf_reconv"], [RECONV_PSF_G[0]])
+
     # successful types are untouched, but share the object's mcal_flags
     npt.assert_allclose(out["noshear"]["flux"], [100.0])
     assert out["noshear"]["flags"] == [0]

tests/module/test_psf_averaging_properties.py

@@ -20,6 +20,7 @@
 
 import numpy as np
 import numpy.testing as npt
+import pytest
 from astropy.io import fits
 from hypothesis import given, settings
 from hypothesis import strategies as st
@@ -199,6 +200,31 @@ def test_single_survivor_returns_its_own_value(survivor, flagged_specs):
     assert out["n_epoch"] == 1
 
 
+@given(
+    st.lists(
+        st.tuples(_weight, st.integers(min_value=1, max_value=255)),
+        min_size=1, max_size=5,
+    )
+)
+@settings(deadline=None, max_examples=25)
+def test_all_epochs_failed_raises_zero_division(flagged_specs):
+    """All-epochs-failed is the contract that wsum == 0 raises, not returns 0/NaN.
+
+    An object whose PSF fit failed on every epoch reaches both PSF-family
+    wrappers (reconv and orig) with zero survivors, so ``wsum`` stays 0. The
+    core must raise ``ZeroDivisionError`` rather than silently divide — a
+    regression returning NaN/0 would be caught here. Every input epoch is
+    flagged (``flags != 0``) with an arbitrary positive weight, so no survivor
+    can rescue the sum.
+    """
+    flagged = [
+        (_result(np.nan, np.nan, 1e6, np.nan, np.nan, np.nan, flags=f), w)
+        for w, f in flagged_specs
+    ]
+    with pytest.raises(ZeroDivisionError):
+        _average_psf_fits(flagged)
+
+
 # --------------------------------------------------------------------------- #
 # (d) make_cat: an obj_id absent from the ngmix cat keeps its sentinel fill.
 # --------------------------------------------------------------------------- #