From 0e70ec72db159732679185c6102b6671729cc0b0 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Tue, 21 Jul 2026 14:37:33 +0200 Subject: [PATCH 1/8] feat(psfex-interp): add spin-2 fourth moments and rho4 in sky frame Extend the sky-coordinate HSM measurement (PR #812) to also store fourth-order moments for PSF and star: - spin-2 combinations HSM_M4_1_* = M40 - M04 and HSM_M4_2_* = 2*(M13 + M31), ported from the PSFHOME/moment4 math, and - galsim's spin-0 HSM_RHO4_* (moments_rho4). Frame handling is the key change over the moment4 port. Whitening with the symmetric sqrt(inv(M)) leaves the whitened axes aligned with the frame in which M and the pixel grid are built, so the spin-2 combinations are defined relative to those axes. Building them in the pixel frame ties them to the CCD orientation. Instead _fourth_moments builds everything in the world frame: M from the use_sky_coords world ellipticity/size, and the grid via the local WCS Jacobian mapping pixel offsets to world offsets about the world centroid (moments_centroid is relative to the stamp true_center). The result is invariant to a re-orientation of the pixel frame viewing the same sky object. HSM failures fill the spin-2 terms with 0 and pass rho4 through (-1); the FLAG column stays the source of truth. Masked star pixels are zeroed before the moment sum. Columns are added to the single-epoch, validation and multi-epoch output paths. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01K98geWHTuT62whqjKocfCm --- .../psfex_interp_package/psfex_interp.py | 157 +++++++++++++++++- 1 file changed, 154 insertions(+), 3 deletions(-) diff --git a/src/shapepipe/modules/psfex_interp_package/psfex_interp.py b/src/shapepipe/modules/psfex_interp_package/psfex_interp.py index b6da6623f..044ed0092 100644 --- a/src/shapepipe/modules/psfex_interp_package/psfex_interp.py +++ b/src/shapepipe/modules/psfex_interp_package/psfex_interp.py @@ -2,7 +2,7 @@ This module computes the PSFs from a PSFEx model at several galaxy positions. -:Authors: Morgan Schmitz, Axel Guinot, Martin Kilbinger +:Authors: Morgan Schmitz, Axel Guinot, Martin Kilbinger, Sacha Guerrini """ @@ -10,6 +10,7 @@ import re import numpy as np +import scipy.linalg as alg from astropy.io import fits from cs_util import size as cs_size from sqlitedict import SqliteDict @@ -62,6 +63,117 @@ def local_wcs_list(wcs, positions): ] +def _fourth_moments(image, moms, wcs=None): + r"""Fourth-Order Moments. + + Compute the spin-2 fourth-moment combinations and galsim's spin-0 + ``moments_rho4`` for a single object, given its HSM adaptive-moment result. + + The spin-2 combinations follow the PSFHOME convention: whiten the pixel grid + with the object's own second-moment matrix ``M`` (so a matched elliptical + Gaussian becomes a unit circle), Gaussian-weight it, form the centred + ``p + q = 4`` moments ``M_pq`` over the whitened coordinates ``(u, v)``, and + return :: + + fourth_moment_1 = M_40 - M_04 + fourth_moment_2 = 2 * (M_13 + M_31) + + **Frame handling.** Whitening with the *symmetric* matrix square root + ``sqrt(inv(M))`` leaves the whitened axes aligned with the coordinate frame + in which ``M`` and the grid are built (it rescales along the object's + principal axes but does not rotate the frame). The spin-2 combinations are + therefore measured relative to *those* axes. To make them a property of the + sky and not of the CCD orientation, everything must live in one frame — the + world frame: + + * ``moms`` comes from ``FindAdaptiveMom(use_sky_coords=True)``, so + ``observed_shape``, ``moments_sigma`` (arcsec) and ``moments_centroid`` + are already in world ``(u, v)`` coordinates. ``M`` built from them is the + world-frame second-moment matrix. + * The pixel grid is mapped to world offsets about the centroid with the local + WCS Jacobian ``J``: ``world = J @ (pixel - true_center)``. galsim reports + ``moments_centroid`` relative to the stamp ``true_center``, so subtracting + it gives ``J @ (pixel - pixel_centroid)`` — the world offset from the + centroid, in the same arcsec units as ``M``. + + The result is invariant under a re-orientation of the pixel frame (a + different WCS rotation viewing the same sky object) — the science guarantee + that motivates the sky-coordinate measurement. When ``wcs`` is ``None`` the + computation is done consistently in the pixel frame instead (the axes are + then the CCD axes and the result is *not* frame-invariant); this path exists + only for the legacy pixel-frame branch. + + Parameters + ---------- + image : numpy.ndarray + Object postage stamp (PSF or star vignet), masked pixels already zeroed. + moms : galsim.hsm.ShapeData + Adaptive-moment result for ``image``. Measured with + ``use_sky_coords=True`` when ``wcs`` is provided. + wcs : galsim.JacobianWCS, optional + Local WCS at the object position, matching the one attached to the image + for the ``use_sky_coords`` measurement. ``None`` selects the pixel frame. + + Returns + ------- + tuple of float + ``(fourth_moment_1, fourth_moment_2, moments_rho4)``. On an HSM failure + (``moms.error_message`` set) the spin-2 terms are filled with ``0.0`` and + ``moments_rho4`` is passed through (``-1`` on failure); the accompanying + ``FLAG`` column is the source of truth for validity. + + """ + if moms.error_message: + return 0.0, 0.0, moms.moments_rho4 + + ny, nx = image.shape + # 1-indexed pixel-centre coordinates (galsim Image origin is (1, 1)). + y_grid, x_grid = np.mgrid[:ny, :nx] + 1.0 + + # Second-moment matrix M with det(M) = sigma**4, built in the measurement + # frame from the (distortion) ellipticity and size. In sky mode these are + # world-frame quantities, so M is the world-frame second-moment matrix. + e1 = moms.observed_shape.e1 + e2 = moms.observed_shape.e2 + sigma4 = moms.moments_sigma**4 + c = (1 + e1) / (1 - e1) + M = np.zeros((2, 2)) + M[1, 1] = np.sqrt(sigma4 / (c - 0.25 * e2**2 * (1 + c) ** 2)) + M[0, 0] = c * M[1, 1] + M[0, 1] = M[1, 0] = 0.5 * e2 * (M[0, 0] + M[1, 1]) + + if wcs is not None: + # Map the pixel grid to world offsets about the world centroid. + x0 = 0.5 * (nx + 1) + y0 = 0.5 * (ny + 1) + pix = np.array([x_grid - x0, y_grid - y0]) + jac = wcs.jacobian().getMatrix() # [[dudx, dudy], [dvdx, dvdy]] + world = np.einsum("ij,jqp->iqp", jac, pix) + u = world[0] - moms.moments_centroid.x + v = world[1] - moms.moments_centroid.y + else: + # Pixel frame: offsets straight from the pixel-frame centroid. + u = x_grid - moms.moments_centroid.x + v = y_grid - moms.moments_centroid.y + + # Whiten with the symmetric square root of inv(M); M is SPD so the principal + # square root is real (drop the ~1e-16 imaginary residue sqrtm returns). + sqrt_inv_M = np.real(alg.sqrtm(np.linalg.inv(M))) + std_pos = np.einsum("ij,jqp->iqp", sqrt_inv_M, np.array([u, v])) + std_x, std_y = std_pos[0], std_pos[1] + + weight = np.exp(-0.5 * (std_x**2 + std_y**2)) + image_weight = weight * image + normalization = np.sum(image_weight) + + def _m(p, q): + return np.sum(image_weight * std_x**p * std_y**q) / normalization + + fourth_moment_1 = _m(4, 0) - _m(0, 4) + fourth_moment_2 = 2 * (_m(1, 3) + _m(3, 1)) + return fourth_moment_1, fourth_moment_2, moms.moments_rho4 + + class PSFExInterpolator(object): """The PSFEx Interpolator Class. @@ -355,6 +467,14 @@ def _get_psfshapes(self, wcs_list=None): returned are already in sky coordinates and need no post-hoc rotation. When ``None``, moments are measured in the pixel frame. + Notes + ----- + In addition to the second-moment shape, the spin-2 fourth-moment + combinations and galsim's spin-0 ``moments_rho4`` are stored (columns 4, + 5, 6 of ``psf_shapes``). With a ``wcs_list`` these are computed in the + world frame so they are invariant to the CCD orientation; see + :func:`_fourth_moments`. + """ if import_fail: raise ImportError("Galsim is required to get shapes information") @@ -371,6 +491,7 @@ def _get_psfshapes(self, wcs_list=None): hsm.FindAdaptiveMom(Image(psf), strict=False) for psf in self.interp_PSFs ] + wcs_list = [None] * len(self.interp_PSFs) self.psf_shapes = np.array( [ @@ -379,8 +500,9 @@ def _get_psfshapes(self, wcs_list=None): moms.observed_shape.g2, moms.moments_sigma, int(bool(moms.error_message)), + *_fourth_moments(psf, moms, wcs), ] - for moms in psf_moms + for psf, moms, wcs in zip(self.interp_PSFs, psf_moms, wcs_list) ] ) @@ -403,6 +525,9 @@ def _write_output(self): "HSM_G2_PSF": self.psf_shapes[:, 1], "HSM_T_PSF": cs_size.sigma_to_T(self.psf_shapes[:, 2]), "HSM_FLAG_PSF": self.psf_shapes[:, 3].astype(int), + "HSM_M4_1_PSF": self.psf_shapes[:, 4], + "HSM_M4_2_PSF": self.psf_shapes[:, 5], + "HSM_RHO4_PSF": self.psf_shapes[:, 6], } else: data = {"VIGNET": self.interp_PSFs} @@ -484,6 +609,13 @@ def _get_starshapes(self, star_vign, wcs_list=None): provided, moments are measured in world coordinates (``use_sky_coords=True``). + Notes + ----- + As in :meth:`_get_psfshapes`, the spin-2 fourth-moment combinations and + galsim's ``moments_rho4`` are stored (columns 4, 5, 6 of ``star_shapes``), + world-frame when a ``wcs_list`` is given. Masked pixels are zeroed before + the fourth-moment sum so they do not contribute. + """ if import_fail: raise ImportError("Galsim is required to get shapes information") @@ -508,6 +640,7 @@ def _get_starshapes(self, star_vign, wcs_list=None): ) for star, mask in zip(star_vign, masks) ] + wcs_list = [None] * len(star_vign) self.star_shapes = np.array( [ @@ -516,8 +649,11 @@ def _get_starshapes(self, star_vign, wcs_list=None): moms.observed_shape.g2, moms.moments_sigma, int(bool(moms.error_message)), + *_fourth_moments(np.where(mask == 1, 0.0, star), moms, wcs), ] - for moms in star_moms + for star, mask, moms, wcs in zip( + star_vign, masks, star_moms, wcs_list + ) ] ) @@ -598,10 +734,16 @@ def _write_output_validation(self, star_dict, psfex_cat_dict): "HSM_G2_PSF": self.psf_shapes[:, 1], "HSM_T_PSF": cs_size.sigma_to_T(self.psf_shapes[:, 2]), "HSM_FLAG_PSF": self.psf_shapes[:, 3].astype(int), + "HSM_M4_1_PSF": self.psf_shapes[:, 4], + "HSM_M4_2_PSF": self.psf_shapes[:, 5], + "HSM_RHO4_PSF": self.psf_shapes[:, 6], "HSM_G1_STAR": self.star_shapes[:, 0], "HSM_G2_STAR": self.star_shapes[:, 1], "HSM_T_STAR": cs_size.sigma_to_T(self.star_shapes[:, 2]), "HSM_FLAG_STAR": self.star_shapes[:, 3].astype(int), + "HSM_M4_1_STAR": self.star_shapes[:, 4], + "HSM_M4_2_STAR": self.star_shapes[:, 5], + "HSM_RHO4_STAR": self.star_shapes[:, 6], } data = {**data, **star_dict} @@ -874,6 +1016,15 @@ def _interpolate_me(self): shape_dict["HSM_FLAG_PSF"] = final_list[j][2][ where_res[0] ][3] + shape_dict["HSM_M4_1_PSF"] = final_list[j][2][ + where_res[0] + ][4] + shape_dict["HSM_M4_2_PSF"] = final_list[j][2][ + where_res[0] + ][5] + shape_dict["HSM_RHO4_PSF"] = final_list[j][2][ + where_res[0] + ][6] output_dict[id_tmp][final_list[j][3][where_res[0]]][ "SHAPES" ] = shape_dict From df07bd65326ea5e229ace9b7eb0b58e901de3743 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Tue, 21 Jul 2026 14:37:33 +0200 Subject: [PATCH 2/8] test(psfex-interp): analytic and frame-invariance tests for fourth moments Drive the real _get_psfshapes / _get_starshapes on galsim stamps rendered through a nontrivial local WCS (scale + rotation + shear): - round and elliptical Gaussians: whitened spin-2 fourth moments ~ 0 and rho4 ~ 2 (whitening circularises any single sheared circular profile, so an elliptical Moffat would also vanish -- a two-Gaussian composite is used to get a genuine non-zero spin-2 signal), - frame invariance: one sky object rendered under WCS orientations of 0/28/63/-45/90 deg gives the same world-frame fourth moments, on both the PSF and star paths, - PSF/star agreement on a clean stamp, and HSM-failure sentinel handling. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01K98geWHTuT62whqjKocfCm --- tests/module/test_psf_fourth_moments.py | 189 ++++++++++++++++++++++++ 1 file changed, 189 insertions(+) create mode 100644 tests/module/test_psf_fourth_moments.py diff --git a/tests/module/test_psf_fourth_moments.py b/tests/module/test_psf_fourth_moments.py new file mode 100644 index 000000000..0cd2a82bd --- /dev/null +++ b/tests/module/test_psf_fourth_moments.py @@ -0,0 +1,189 @@ +"""UNIT TESTS FOR PSF/STAR FOURTH-ORDER MOMENTS. + +Cover the spin-2 fourth-moment combinations (``HSM_M4_1_*`` / ``HSM_M4_2_*``) +and galsim's spin-0 ``moments_rho4`` (``HSM_RHO4_*``) added to the PSFEx-interp +sky-coordinate shape measurement (:mod:`shapepipe.modules.psfex_interp_package`). + +The measurement whitens the object by its own second-moment matrix, so the +predictions are analytic for elliptically-symmetric profiles: + +* a round or elliptical **Gaussian** matches its own adaptive weight, so the + whitened profile is an isotropic Gaussian: the spin-2 fourth moments vanish + and ``rho4`` equals the Gaussian value 2. +* to get *non-zero* spin-2 fourth moments the profile must not be a single + sheared circular profile (whitening circularises any such profile). We use a + sum of two coaxial Gaussians of different ellipticity, whose isophote shape + changes with radius. + +The science guarantee is **frame invariance**: because the whitening and the +grid are built in the world frame (see ``_fourth_moments``), one sky object +rendered onto two differently oriented pixel grids gives the same fourth +moments. That is the key test below. + +The WCS construction mirrors ``test_hsm_sky_coords.py``: a CD-matrix-like local +Jacobian with a realistic scale, rotation and small shear. +""" + +import galsim +import numpy as np +import numpy.testing as npt +import pytest + +from shapepipe.modules.psfex_interp_package.psfex_interp import ( + PSFExInterpolator, + _fourth_moments, +) + +# Column layout of psf_shapes / star_shapes. +M4_1, M4_2, RHO4 = 4, 5, 6 + +_STAMP = 101 + + +def make_wcs(theta_deg, scale=0.187, g1=0.06, g2=-0.04): + """A CD-matrix-like local WCS: scale, rotation ``theta_deg``, small shear.""" + th = np.deg2rad(theta_deg) + c, s = np.cos(th), np.sin(th) + mat = galsim.Shear(g1=g1, g2=g2).getMatrix() * scale @ np.array( + [[c, -s], [s, c]] + ) + return galsim.JacobianWCS(mat[0, 0], mat[0, 1], mat[1, 0], mat[1, 1]) + + +def render(profile, wcs): + """Render a world-frame profile onto the pixel grid defined by ``wcs``.""" + return profile.drawImage( + nx=_STAMP, ny=_STAMP, wcs=wcs, method="no_pixel" + ).array + + +def psf_row(profile, wcs): + """Run the real ``_get_psfshapes`` on one stamp; return its shape row.""" + interp = object.__new__(PSFExInterpolator) + interp.interp_PSFs = [render(profile, wcs)] + interp._get_psfshapes([wcs]) + return interp.psf_shapes[0] + + +def star_row(profile, wcs, mask=None): + """Run the real ``_get_starshapes`` on one star vignet; return its row.""" + stamp = render(profile, wcs) + if mask is not None: + stamp = np.where(mask, -1e30, stamp) + interp = object.__new__(PSFExInterpolator) + interp._get_starshapes(np.array([stamp]), [wcs]) + return interp.star_shapes[0] + + +# Sum of two coaxial Gaussians of different ellipticity: not elliptically +# symmetric, so its whitened spin-2 fourth moments are non-zero. +def composite(beta_deg=0.0): + prof = galsim.Gaussian(sigma=0.4).shear(e1=0.5) + galsim.Gaussian( + sigma=1.0 + ).shear(e1=0.1) + if beta_deg: + prof = prof.rotate(beta_deg * galsim.degrees) + return prof + + +# --------------------------------------------------------------------------- +# Analytic checks on Gaussians. +# --------------------------------------------------------------------------- + + +def test_round_gaussian_moments_vanish_rho4_is_gaussian(): + """Round Gaussian: spin-2 fourth moments ~ 0, rho4 ~ 2 (Gaussian value).""" + row = psf_row(galsim.Gaussian(sigma=0.6), make_wcs(0.0)) + npt.assert_allclose(row[M4_1], 0.0, atol=1e-4) + npt.assert_allclose(row[M4_2], 0.0, atol=1e-4) + npt.assert_allclose(row[RHO4], 2.0, rtol=1e-3) + + +@pytest.mark.parametrize( + "e1, e2", + [(0.3, 0.0), (0.0, 0.25), (0.3, 0.15), (-0.2, -0.3)], +) +def test_elliptical_gaussian_whitening_vanishes(e1, e2): + """Elliptical Gaussian: whitening by its own 2nd moment circularises it, so + the spin-2 fourth moments vanish regardless of (e1, e2).""" + row = psf_row(galsim.Gaussian(sigma=0.6).shear(e1=e1, e2=e2), make_wcs(0.0)) + npt.assert_allclose(row[M4_1], 0.0, atol=1e-4) + npt.assert_allclose(row[M4_2], 0.0, atol=1e-4) + npt.assert_allclose(row[RHO4], 2.0, rtol=1e-3) + + +def test_composite_has_nonzero_spin2(): + """Coaxial two-Gaussian composite: non-Gaussian radial shape leaves a real + spin-2 fourth moment. Axis-aligned, so only M4_1 is excited (M4_2 ~ 0).""" + row = psf_row(composite(), make_wcs(0.0)) + assert abs(row[M4_1]) > 1e-2 + npt.assert_allclose(row[M4_2], 0.0, atol=1e-4) + + +# --------------------------------------------------------------------------- +# Frame invariance -- the science guarantee. +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize("theta", [28.0, 63.0, -45.0, 90.0]) +def test_psf_fourth_moments_frame_invariant(theta): + """One sky object rendered under two WCS orientations gives the same + world-frame fourth moments (PSF path).""" + prof = composite(beta_deg=30.0) # rotated so both M4_1 and M4_2 are excited + ref = psf_row(prof, make_wcs(0.0)) + rot = psf_row(prof, make_wcs(theta)) + + assert abs(ref[M4_1]) > 1e-2 and abs(ref[M4_2]) > 1e-2 # non-trivial + npt.assert_allclose(rot[M4_1], ref[M4_1], rtol=1e-3, atol=1e-5) + npt.assert_allclose(rot[M4_2], ref[M4_2], rtol=1e-3, atol=1e-5) + npt.assert_allclose(rot[RHO4], ref[RHO4], rtol=1e-4) + + +@pytest.mark.parametrize("theta", [28.0, -45.0]) +def test_star_fourth_moments_frame_invariant(theta): + """Same guarantee on the star path (which also handles a bad-pixel mask).""" + prof = composite(beta_deg=30.0) + ref = star_row(prof, make_wcs(0.0)) + rot = star_row(prof, make_wcs(theta)) + + npt.assert_allclose(rot[M4_1], ref[M4_1], rtol=1e-3, atol=1e-5) + npt.assert_allclose(rot[M4_2], ref[M4_2], rtol=1e-3, atol=1e-5) + npt.assert_allclose(rot[RHO4], ref[RHO4], rtol=1e-4) + + +def test_psf_and_star_paths_agree(): + """PSF and star measurement of the same clean stamp agree.""" + prof = composite(beta_deg=30.0) + wcs = make_wcs(17.0) + p = psf_row(prof, wcs) + s = star_row(prof, wcs) + npt.assert_allclose(s[M4_1], p[M4_1], rtol=1e-4, atol=1e-6) + npt.assert_allclose(s[M4_2], p[M4_2], rtol=1e-4, atol=1e-6) + npt.assert_allclose(s[RHO4], p[RHO4], rtol=1e-4) + + +# --------------------------------------------------------------------------- +# Failure handling. +# --------------------------------------------------------------------------- + + +def test_hsm_failure_fills_sentinels(): + """A flat stamp fails HSM: FLAG is set, spin-2 filled with 0, rho4 = -1.""" + interp = object.__new__(PSFExInterpolator) + interp.interp_PSFs = [np.zeros((21, 21))] + interp._get_psfshapes([make_wcs(0.0)]) + row = interp.psf_shapes[0] + assert int(row[3]) == 1 # FLAG + assert row[M4_1] == 0.0 + assert row[M4_2] == 0.0 + assert row[RHO4] == -1.0 + + +def test_fourth_moments_failure_guard(): + """``_fourth_moments`` short-circuits on a ShapeData carrying an error.""" + failed = galsim.hsm.ShapeData(error_message="boom", moments_rho4=-1.0) + assert _fourth_moments(np.zeros((5, 5)), failed, make_wcs(0.0)) == ( + 0.0, + 0.0, + -1.0, + ) From 11e2f4ce34d4f24c3eb99968a52179b9dd0b9212 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Tue, 21 Jul 2026 14:53:38 +0200 Subject: [PATCH 3/8] fix(merge-starcat): carry HSM fourth moments through PSFEx merge MergeStarCatPSFEX read and wrote only HSM_G1/G2/T/FLAG for PSF and STAR, so the spin-2 fourth moments and rho4 written per-tile by psfex_interp were dropped at the merge step. The merged full_starcat is the catalogue PSF validation consumes, so the feature was inert past the per-tile stage. Read and write HSM_M4_1/HSM_M4_2/HSM_RHO4 for both PSF and STAR, matching the per-tile column grammar. Avoids the moment4 reference's copy-paste bug (m4_2_psf now reads the _2 column) and adds RHO4, which the reference omitted. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01K98geWHTuT62whqjKocfCm --- .../modules/merge_starcat_package/merge_starcat.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/src/shapepipe/modules/merge_starcat_package/merge_starcat.py b/src/shapepipe/modules/merge_starcat_package/merge_starcat.py index 7bc0cb76b..58fdd3767 100644 --- a/src/shapepipe/modules/merge_starcat_package/merge_starcat.py +++ b/src/shapepipe/modules/merge_starcat_package/merge_starcat.py @@ -559,7 +559,9 @@ def process(self): """ x, y, ra, dec = [], [], [], [] g1_psf, g2_psf, size_psf = [], [], [] + m4_1_psf, m4_2_psf, rho4_psf = [], [], [] g1, g2, size = [], [], [] + m4_1_star, m4_2_star, rho4_star = [], [], [] flag_psf, flag_star = [], [] mag, snr, psfex_acc = [], [], [] ccd_nb = [] @@ -588,9 +590,15 @@ def process(self): g1_psf += list(data_j["HSM_G1_PSF"]) g2_psf += list(data_j["HSM_G2_PSF"]) size_psf += list(data_j["HSM_T_PSF"]) + m4_1_psf += list(data_j["HSM_M4_1_PSF"]) + m4_2_psf += list(data_j["HSM_M4_2_PSF"]) + rho4_psf += list(data_j["HSM_RHO4_PSF"]) g1 += list(data_j["HSM_G1_STAR"]) g2 += list(data_j["HSM_G2_STAR"]) size += list(data_j["HSM_T_STAR"]) + m4_1_star += list(data_j["HSM_M4_1_STAR"]) + m4_2_star += list(data_j["HSM_M4_2_STAR"]) + rho4_star += list(data_j["HSM_RHO4_STAR"]) # flags flag_psf += list(data_j["HSM_FLAG_PSF"]) @@ -636,9 +644,15 @@ def process(self): "HSM_G1_PSF": g1_psf, "HSM_G2_PSF": g2_psf, "HSM_T_PSF": size_psf, + "HSM_M4_1_PSF": m4_1_psf, + "HSM_M4_2_PSF": m4_2_psf, + "HSM_RHO4_PSF": rho4_psf, "HSM_G1_STAR": g1, "HSM_G2_STAR": g2, "HSM_T_STAR": size, + "HSM_M4_1_STAR": m4_1_star, + "HSM_M4_2_STAR": m4_2_star, + "HSM_RHO4_STAR": rho4_star, "HSM_FLAG_PSF": flag_psf, "HSM_FLAG_STAR": flag_star, "MAG": mag, From ce5ed528ff85e0cc83fe93ae697ed778a1884a19 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Tue, 21 Jul 2026 14:53:52 +0200 Subject: [PATCH 4/8] test(psfex-interp): cover bad-pixel masking and centroid transform Two paths the summary claimed but the suite left unexercised: - Bad-pixel masking: every star fixture had an all-zero mask, so zeroing masked pixels before the fourth-moment sum was a no-op in the suite. Add a test that plants the -1e30 sentinel in outskirt pixels, asserts the masked vignet reproduces the clean result, and asserts the same sentinels fed through un-zeroed blow the fourth moments up by ~5 orders (mutation: removing the zeroing makes M4_1 go -0.03 -> -4278). - Centroid transform: every fixture was drawn centered, so moments_centroid ~ 0 and the sky-frame recentering was dead. Add a world-frame-shift case; the fourth moments stay invariant only because the transformed centroid is subtracted (mutation: dropping the subtraction fails both angles). Also tighten the frame-invariance tolerances from rtol=1e-3/atol=1e-5 to rtol=1e-5/atol=1e-6 (measured cross-orientation agreement is ~1e-8), so a subtle partial-frame error at the 1e-4 level now bites. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01K98geWHTuT62whqjKocfCm --- tests/module/test_psf_fourth_moments.py | 78 +++++++++++++++++++++++-- 1 file changed, 72 insertions(+), 6 deletions(-) diff --git a/tests/module/test_psf_fourth_moments.py b/tests/module/test_psf_fourth_moments.py index 0cd2a82bd..4782ca5f4 100644 --- a/tests/module/test_psf_fourth_moments.py +++ b/tests/module/test_psf_fourth_moments.py @@ -134,9 +134,9 @@ def test_psf_fourth_moments_frame_invariant(theta): rot = psf_row(prof, make_wcs(theta)) assert abs(ref[M4_1]) > 1e-2 and abs(ref[M4_2]) > 1e-2 # non-trivial - npt.assert_allclose(rot[M4_1], ref[M4_1], rtol=1e-3, atol=1e-5) - npt.assert_allclose(rot[M4_2], ref[M4_2], rtol=1e-3, atol=1e-5) - npt.assert_allclose(rot[RHO4], ref[RHO4], rtol=1e-4) + npt.assert_allclose(rot[M4_1], ref[M4_1], rtol=1e-5, atol=1e-6) + npt.assert_allclose(rot[M4_2], ref[M4_2], rtol=1e-5, atol=1e-6) + npt.assert_allclose(rot[RHO4], ref[RHO4], rtol=1e-5) @pytest.mark.parametrize("theta", [28.0, -45.0]) @@ -146,9 +146,31 @@ def test_star_fourth_moments_frame_invariant(theta): ref = star_row(prof, make_wcs(0.0)) rot = star_row(prof, make_wcs(theta)) - npt.assert_allclose(rot[M4_1], ref[M4_1], rtol=1e-3, atol=1e-5) - npt.assert_allclose(rot[M4_2], ref[M4_2], rtol=1e-3, atol=1e-5) - npt.assert_allclose(rot[RHO4], ref[RHO4], rtol=1e-4) + npt.assert_allclose(rot[M4_1], ref[M4_1], rtol=1e-5, atol=1e-6) + npt.assert_allclose(rot[M4_2], ref[M4_2], rtol=1e-5, atol=1e-6) + npt.assert_allclose(rot[RHO4], ref[RHO4], rtol=1e-5) + + +@pytest.mark.parametrize("theta", [0.0, 40.0]) +def test_fourth_moments_invariant_under_world_shift(theta): + """Off-center object: the sky-frame centroid subtraction makes the fourth + moments invariant under a world-frame translation of the source. + + Every other fixture is drawn centered, so ``moments_centroid`` is ~0 and the + recentering in ``_fourth_moments`` is a no-op. Real star vignets are only + approximately centered, so this case shifts the source in the *world* frame + (arcsec) and asserts the fourth moments are unchanged -- exercising the + centroid transform that is otherwise dead in the suite. + """ + prof = composite(beta_deg=30.0) + wcs = make_wcs(theta) + ref = psf_row(prof, wcs) + shifted = psf_row(prof.shift(0.9, -1.3), wcs) # world-frame shift, arcsec + + assert abs(ref[M4_1]) > 1e-2 and abs(ref[M4_2]) > 1e-2 # non-trivial + npt.assert_allclose(shifted[M4_1], ref[M4_1], rtol=1e-4, atol=1e-6) + npt.assert_allclose(shifted[M4_2], ref[M4_2], rtol=1e-4, atol=1e-6) + npt.assert_allclose(shifted[RHO4], ref[RHO4], rtol=1e-4) def test_psf_and_star_paths_agree(): @@ -162,6 +184,50 @@ def test_psf_and_star_paths_agree(): npt.assert_allclose(s[RHO4], p[RHO4], rtol=1e-4) +# --------------------------------------------------------------------------- +# Bad-pixel masking on the star path. +# --------------------------------------------------------------------------- + + +def test_masked_pixels_are_zeroed_before_fourth_moment_sum(): + """Star path: pixels carrying the ``-1e30`` bad-pixel sentinel are zeroed + before the fourth-moment sum, so a masked vignet reproduces the clean + result -- and the zeroing is load-bearing. + + Without the zeroing the sum runs over the raw ``-1e30`` sentinels and the + fourth moments become garbage (the sum over an un-zeroed stamp is order + ``-1e30``). We mask low-flux outskirt pixels: dropping them leaves the clean + measurement essentially unchanged, while feeding the same sentinels through + un-zeroed blows the result up by many orders of magnitude. + """ + prof = composite(beta_deg=30.0) + wcs = make_wcs(17.0) + clean = star_row(prof, wcs) + + # A block of outskirt pixels, offset off-axis so it does not cancel in M4_1. + mask = np.zeros((_STAMP, _STAMP), dtype=bool) + c = _STAMP // 2 + mask[c + 18 : c + 22, c + 12 : c + 15] = True + + masked = star_row(prof, wcs, mask=mask) + + # Zeroing => the masked measurement reproduces the clean one. + npt.assert_allclose(masked[M4_1], clean[M4_1], rtol=1e-5, atol=1e-6) + npt.assert_allclose(masked[M4_2], clean[M4_2], rtol=1e-5, atol=1e-6) + npt.assert_allclose(masked[RHO4], clean[RHO4], rtol=1e-5) + + # Contrast: the same sentinel stamp with the zeroing removed is garbage. + sentinel = np.where(mask, -1e30, render(prof, wcs)) + moms = galsim.hsm.FindAdaptiveMom( + galsim.Image(sentinel, wcs=wcs), + badpix=galsim.Image(mask.astype(float)), + strict=False, + use_sky_coords=True, + ) + raw = _fourth_moments(sentinel, moms, wcs) # no zeroing applied + assert abs(raw[0] - clean[M4_1]) > 1.0 + + # --------------------------------------------------------------------------- # Failure handling. # --------------------------------------------------------------------------- From 6864b20cf99834022627eb180048039bab6818f7 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 26 Sep 2026 02:20:34 +0200 Subject: [PATCH 5/8] feat(make_cat): carry PSF fourth moments into the per-epoch HSM_*_PSF_n slots psfex_interp's multi-epoch SHAPES dict already holds HSM_M4_1/M4_2/RHO4_PSF; _save_psf_data now copies them alongside G1/G2/T/FLAG, driven by one column table (name, fill, dtype) instead of six hand-rolled slot initialisers. A SHAPES dict from an older producer leaves the slot at its fill. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_0151inhpkat9cPPdMTG65Rwm --- .../modules/make_cat_package/make_cat.py | 93 +++++++------------ tests/module/test_make_cat.py | 51 ++++++++-- 2 files changed, 74 insertions(+), 70 deletions(-) diff --git a/src/shapepipe/modules/make_cat_package/make_cat.py b/src/shapepipe/modules/make_cat_package/make_cat.py index 062ac5719..303384049 100644 --- a/src/shapepipe/modules/make_cat_package/make_cat.py +++ b/src/shapepipe/modules/make_cat_package/make_cat.py @@ -619,6 +619,12 @@ def _save_psf_data(self, galaxy_psf_path): Save the PSF catalogue into the final one. + @sc [label:schema] psf-epoch-slot-columns + The per-epoch families are ``psf_shape_cols`` plus ``EXP_ID``/``CCD``, + all slot-aligned by the producer's epoch enumeration; the names must + be a subset of what ``PSFExInterpolator._interpolate_me`` writes into + ``SHAPES`` (``test_hsm_column_seams``). + Parameters ---------- galaxy_psf_path : str @@ -628,49 +634,28 @@ def _save_psf_data(self, galaxy_psf_path): galaxy_psf_cat = SqliteDict(galaxy_psf_path) max_epoch = np.max(self._final_cat_file.get_data()["N_EPOCH"]) + 1 + n_obj = len(self._obj_id) - self._output_dict = { - f"HSM_G1_PSF_{idx + 1}": np.ones(len(self._obj_id)) * -10.0 - for idx in range(max_epoch) - } - self._output_dict = { - **self._output_dict, - **{ - f"HSM_G2_PSF_{idx + 1}": np.ones(len(self._obj_id)) * -10.0 - for idx in range(max_epoch) - }, - } - self._output_dict = { - **self._output_dict, - **{ - f"HSM_T_PSF_{idx + 1}": np.zeros(len(self._obj_id)) - for idx in range(max_epoch) - }, - } - self._output_dict = { - **self._output_dict, - **{ - f"HSM_FLAG_PSF_{idx + 1}": np.ones( - len(self._obj_id), dtype="int16" - ) - for idx in range(max_epoch) - }, - } + # Per-epoch PSF shape columns copied from the producer's SHAPES dict: + # (column, empty-slot fill, dtype). HSM_T_PSF already holds T + # (sigma_to_T applied at the producer's _interpolate_me). + psf_shape_cols = [ + ("HSM_G1_PSF", -10.0, float), + ("HSM_G2_PSF", -10.0, float), + ("HSM_T_PSF", 0.0, float), + ("HSM_FLAG_PSF", 1, "int16"), + ("HSM_M4_1_PSF", 0.0, float), + ("HSM_M4_2_PSF", 0.0, float), + ("HSM_RHO4_PSF", -1.0, float), + ] # Per-epoch exposure ID and CCD number, slot-aligned with HSM_*_PSF_n; # -1 marks an empty slot (the CCD_N sentinel convention). + epoch_id_cols = [("EXP_ID", -1, "int32"), ("CCD", -1, "int32")] + self._output_dict = { - **self._output_dict, - **{ - f"EXP_ID_{idx + 1}": np.ones(len(self._obj_id), dtype="int32") * -1 - for idx in range(max_epoch) - }, - } - self._output_dict = { - **self._output_dict, - **{ - f"CCD_{idx + 1}": np.ones(len(self._obj_id), dtype="int32") * -1 - for idx in range(max_epoch) - }, + f"{name}_{idx + 1}": np.full(n_obj, fill, dtype=dtype) + for name, fill, dtype in psf_shape_cols + epoch_id_cols + for idx in range(max_epoch) } for idx, id_tmp in enumerate(self._obj_id): @@ -680,7 +665,7 @@ def _save_psf_data(self, galaxy_psf_path): for epoch, key in enumerate(galaxy_psf_cat[str(id_tmp)].keys()): - gpc_data = galaxy_psf_cat[str(id_tmp)][key] + shapes = galaxy_psf_cat[str(id_tmp)][key]["SHAPES"] # `key` is "-"; reading it in the enumeration that # assigns `epoch` aligns EXP_ID_n/CCD_n with HSM_*_PSF_n by @@ -690,28 +675,14 @@ def _save_psf_data(self, galaxy_psf_path): self._add2dict(f"EXP_ID_{epoch + 1}", int(exp_name), idx) self._add2dict(f"CCD_{epoch + 1}", int(ccd_n), idx) - if gpc_data["SHAPES"]["HSM_FLAG_PSF"] != 0: + if shapes["HSM_FLAG_PSF"] != 0: continue - self._add2dict( - f"HSM_G1_PSF_{epoch + 1}", - gpc_data["SHAPES"]["HSM_G1_PSF"], idx - ) - self._add2dict( - f"HSM_G2_PSF_{epoch + 1}", - gpc_data["SHAPES"]["HSM_G2_PSF"], idx - ) - - # HSM_T_PSF already holds T (sigma_to_T applied at the - # producer's _interpolate_me); read straight through. - self._add2dict( - f"HSM_T_PSF_{epoch + 1}", - gpc_data["SHAPES"]["HSM_T_PSF"], idx - ) - - self._add2dict( - f"HSM_FLAG_PSF_{epoch + 1}", - gpc_data["SHAPES"]["HSM_FLAG_PSF"], idx - ) + for name, fill, _ in psf_shape_cols: + # A SHAPES dict from an older producer may lack the + # fourth-moment keys; the slot then keeps its fill. + self._add2dict( + f"{name}_{epoch + 1}", shapes.get(name, fill), idx + ) galaxy_psf_cat.close() diff --git a/tests/module/test_make_cat.py b/tests/module/test_make_cat.py index 34a1958e8..85ae16de7 100644 --- a/tests/module/test_make_cat.py +++ b/tests/module/test_make_cat.py @@ -313,16 +313,21 @@ def test_save_ngmix_data_matches_module_serialised_catalogue(tmp_path): # --- _save_psf_data: EXP_ID_n / CCD_n alignment with HSM_*_PSF_n (#890) --- -def _psf_epoch(g1, g2, t, flag=0): - """One epoch's interpolated-PSF HSM shape entry (psfex_interp's SHAPES dict).""" - return { - "SHAPES": { - "HSM_G1_PSF": g1, - "HSM_G2_PSF": g2, - "HSM_T_PSF": t, - "HSM_FLAG_PSF": flag, - }, +def _psf_epoch(g1, g2, t, flag=0, m4=None): + """One epoch's interpolated-PSF HSM shape entry (psfex_interp's SHAPES dict). + + ``m4`` is an optional ``(M4_1, M4_2, RHO4)`` triple; omitted, the entry + mimics a producer that predates the fourth-moment columns. + """ + shapes = { + "HSM_G1_PSF": g1, + "HSM_G2_PSF": g2, + "HSM_T_PSF": t, + "HSM_FLAG_PSF": flag, } + if m4 is not None: + shapes.update(zip(("HSM_M4_1_PSF", "HSM_M4_2_PSF", "HSM_RHO4_PSF"), m4)) + return {"SHAPES": shapes} def _write_galaxy_psf_cat(path, per_obj): @@ -448,3 +453,31 @@ def test_save_psf_data_fills_sentinel_for_absent_epochs(tmp_path): assert out[col][0] == -1, col for col in ("EXP_ID_1", "CCD_1", "EXP_ID_2", "CCD_2", "EXP_ID_3", "CCD_3"): assert out[col][1] == -1, col + + +def test_save_psf_data_carries_fourth_moments_per_epoch(tmp_path): + """HSM_M4_1/M4_2/RHO4_PSF_n ride the same slots as HSM_G1_PSF_n. + + The fourth-moment columns psfex_interp writes into SHAPES (shapepipe#697) + land per epoch; a SHAPES dict without them (an older producer) leaves + the slot at its fill, as does an unused slot. + """ + galaxy_psf_path = tmp_path / "galaxy_psf.sqlite" + per_obj = { + 101: { + "2113864-7": _psf_epoch(0.01, 0.02, 0.5, m4=(0.11, -0.22, 2.05)), + "2113865-3": _psf_epoch(0.03, 0.04, 0.6), + }, + } + _write_galaxy_psf_cat(galaxy_psf_path, per_obj) + + out = _run_save_psf(galaxy_psf_path, [101], n_epoch=[2]) + + npt.assert_allclose(out["HSM_M4_1_PSF_1"], [0.11]) + npt.assert_allclose(out["HSM_M4_2_PSF_1"], [-0.22]) + npt.assert_allclose(out["HSM_RHO4_PSF_1"], [2.05]) + for n in (2, 3): + npt.assert_allclose(out[f"HSM_M4_1_PSF_{n}"], [0.0]) + npt.assert_allclose(out[f"HSM_M4_2_PSF_{n}"], [0.0]) + npt.assert_allclose(out[f"HSM_RHO4_PSF_{n}"], [-1.0]) + npt.assert_allclose(out["HSM_G1_PSF_2"], [0.03]) From b40dfc2fe8c2252cd3988d4719449804010e2fed Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 26 Sep 2026 02:20:34 +0200 Subject: [PATCH 6/8] test: HSM column seam checks + @sc contracts for the fourth-moment chain MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Static AST check that MergeStarCatPSFEX reads exactly the HSM_* set _write_output_validation writes (both read by name with no fallback — the seam behind the #698 copy-paste bug), and that make_cat's per-epoch columns are a subset of what _interpolate_me writes into SHAPES. @sc contracts at the four declarations name the invariants; the frozen-grammar regex learns the M4_1/M4_2/RHO4 tokens. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_0151inhpkat9cPPdMTG65Rwm --- log_-i | 1 + .../merge_starcat_package/merge_starcat.py | 5 ++ .../psfex_interp_package/psfex_interp.py | 17 +++++ tests/module/test_hsm_column_seams.py | 63 +++++++++++++++++++ tests/module/test_psf_grammar_properties.py | 6 +- 5 files changed, 91 insertions(+), 1 deletion(-) create mode 100644 log_-i create mode 100644 tests/module/test_hsm_column_seams.py diff --git a/log_-i b/log_-i new file mode 100644 index 000000000..989bf77a3 --- /dev/null +++ b/log_-i @@ -0,0 +1 @@ +-i /tmp/pytest-of-cdaley/pytest-80/test_build_map_pole_guard_skip0/corners.txt -o /tmp/pytest-of-cdaley/pytest-80/test_build_map_pole_guard_skip0/cov.hsp -c 32 -n 1024 diff --git a/src/shapepipe/modules/merge_starcat_package/merge_starcat.py b/src/shapepipe/modules/merge_starcat_package/merge_starcat.py index 58fdd3767..0d724d210 100644 --- a/src/shapepipe/modules/merge_starcat_package/merge_starcat.py +++ b/src/shapepipe/modules/merge_starcat_package/merge_starcat.py @@ -556,6 +556,11 @@ def process(self): Process merging. + @sc [label:schema] psfex-starcat-columns-strict + Every ``HSM_*`` column is read by name with no fallback, so the set + read here equals the set ``PSFExInterpolator._write_output_validation`` + writes (``test_hsm_column_seams``). + """ x, y, ra, dec = [], [], [], [] g1_psf, g2_psf, size_psf = [], [], [] diff --git a/src/shapepipe/modules/psfex_interp_package/psfex_interp.py b/src/shapepipe/modules/psfex_interp_package/psfex_interp.py index 044ed0092..63ab5a1bb 100644 --- a/src/shapepipe/modules/psfex_interp_package/psfex_interp.py +++ b/src/shapepipe/modules/psfex_interp_package/psfex_interp.py @@ -66,6 +66,13 @@ def local_wcs_list(wcs, positions): def _fourth_moments(image, moms, wcs=None): r"""Fourth-Order Moments. + @sc [label:frame] fourth-moments-single-frame + ``moms`` and the pixel-grid mapping must live in one frame. With ``wcs``, + ``moms`` comes from ``FindAdaptiveMom(use_sky_coords=True)`` on an image + carrying that same local WCS, and the whitening matrix is built from those + world-frame moments; mixing frames silently breaks the CCD-orientation + invariance that ``test_psf_fourth_moments_frame_invariant`` enforces. + Compute the spin-2 fourth-moment combinations and galsim's spin-0 ``moments_rho4`` for a single object, given its HSM adaptive-moment result. @@ -715,6 +722,11 @@ def _write_output_validation(self, star_dict, psfex_cat_dict): Save computed PSFs and stars to fits file. + @sc [label:schema] psfex-validation-hsm-columns + The ``HSM_*`` columns written here are the exact set + ``MergeStarCatPSFEX.process`` reads without fallback; add or rename on + both sides together (``test_hsm_column_seams`` checks the two agree). + Parameters ---------- star_dict : dict @@ -795,6 +807,11 @@ def _interpolate_me(self): Interpolate PSFs for multi-epoch run. + @sc [label:schema] psfex-me-shapes-columns + The per-epoch ``SHAPES`` dict must carry every column + ``make_cat._save_psf_data`` copies into ``HSM_*_PSF_n`` + (``test_hsm_column_seams`` checks the superset). + Raises ------ KeyError diff --git a/tests/module/test_hsm_column_seams.py b/tests/module/test_hsm_column_seams.py new file mode 100644 index 000000000..c888d0a58 --- /dev/null +++ b/tests/module/test_hsm_column_seams.py @@ -0,0 +1,63 @@ +"""Static seam checks between HSM column producers and consumers. + +Contracts ``psfex-validation-hsm-columns`` / ``psfex-starcat-columns-strict`` +(psfex_interp ↔ merge_starcat) and ``psfex-me-shapes-columns`` / +``psf-epoch-slot-columns`` (psfex_interp ↔ make_cat): the column names are +string literals on both sides, read with no fallback, so the two sides are +compared as sets of ``HSM_*`` literals collected from each function's AST. +""" + +import ast +import inspect +import textwrap + +from shapepipe.modules.make_cat_package.make_cat import SaveCatalogue +from shapepipe.modules.merge_starcat_package.merge_starcat import ( + MergeStarCatPSFEX, +) +from shapepipe.modules.psfex_interp_package.psfex_interp import ( + PSFExInterpolator, +) + + +def _hsm_literals(func, subscript_of=None): + """Set of ``HSM_*`` string literals in ``func``'s source. + + With ``subscript_of``, only literals used as ``["HSM_..."]`` reads + count, so a column that is still *written* under the same name cannot + mask a dropped read. + """ + tree = ast.parse(textwrap.dedent(inspect.getsource(func))) + if subscript_of is None: + nodes = (n for n in ast.walk(tree) if isinstance(n, ast.Constant)) + else: + nodes = ( + n.slice + for n in ast.walk(tree) + if isinstance(n, ast.Subscript) + and isinstance(n.value, ast.Name) + and n.value.id == subscript_of + and isinstance(n.slice, ast.Constant) + ) + return { + n.value + for n in nodes + if isinstance(n.value, str) and n.value.startswith("HSM_") + } + + +def test_merge_starcat_reads_exactly_what_psfex_validation_writes(): + """psfex-validation-hsm-columns == psfex-starcat-columns-strict.""" + written = _hsm_literals(PSFExInterpolator._write_output_validation) + read = _hsm_literals(MergeStarCatPSFEX.process, subscript_of="data_j") + assert written == read, { + "written_not_read": sorted(written - read), + "read_not_written": sorted(read - written), + } + + +def test_make_cat_epoch_columns_come_from_psfex_me_shapes(): + """psf-epoch-slot-columns ⊆ psfex-me-shapes-columns.""" + produced = _hsm_literals(PSFExInterpolator._interpolate_me) + consumed = _hsm_literals(SaveCatalogue._save_psf_data) + assert consumed <= produced, sorted(consumed - produced) diff --git a/tests/module/test_psf_grammar_properties.py b/tests/module/test_psf_grammar_properties.py index 6bf9e4713..6c7553f86 100644 --- a/tests/module/test_psf_grammar_properties.py +++ b/tests/module/test_psf_grammar_properties.py @@ -417,7 +417,7 @@ def test_param_file_ngmix_tokens_are_producible(param_path, obj_ids): r"|NGMIXm?_(?:MCAL_FLAGS|MCAL_TYPES_FAIL|N_EPOCH|NEIGHBOUR_FLAG)" # HSM: g-type, explicit PSF/STAR object, singular FLAG; the multi-epoch # sink in make_cat._save_psf_data appends a bare epoch index. - r"|HSM_(?:G1|G2|T)_(?:PSF|STAR)(?:_\d+)?" + r"|HSM_(?:G1|G2|T|M4_1|M4_2|RHO4)_(?:PSF|STAR)(?:_\d+)?" r"|HSM_FLAG_(?:PSF|STAR)(?:_\d+)?" r")$" ) @@ -438,6 +438,9 @@ def test_param_file_ngmix_tokens_are_producible(param_path, obj_ids): "HSM_G1_PSF_3", # make_cat multi-epoch sink "HSM_T_PSF_2", "HSM_FLAG_STAR_1", + "HSM_M4_1_PSF", # spin-2 fourth moments + rho4 (shapepipe#697) + "HSM_RHO4_STAR", + "HSM_M4_2_PSF_2", ] # Pre-#761 / off-grammar columns the rename REMOVES — each violates exactly @@ -452,6 +455,7 @@ def test_param_file_ngmix_tokens_are_producible(param_path, obj_ids): "HSM_FLAGS_PSF", # plural — HSM is singular FLAG "NGMIX_ELL_PSF_ORIG_NOSHEAR", # packed ellipticity, not G1/G2 "SPREAD_MODEL", # removed entirely, not renamed + "M_4_PSF_1", # pre-#859 fourth-moment naming (shapepipe#698) ] From 9f5b8de2195455ea8c94c541609db13565f83621 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 26 Sep 2026 02:32:40 +0200 Subject: [PATCH 7/8] refactor(psfex-interp): one _hsm_columns writer; out-of-range fourth-moment fills in make_cat _HSM_ROW names the shape-row slots once and _hsm_columns maps them to the HSM_*_{PSF,STAR} grammar (sigma -> T, FLAG -> int), so the single-epoch, validation and multi-epoch writers no longer each index [4],[5],[6] by hand; the seam test takes the produced set from the helper. make_cat fills unmeasured HSM_M4_*_PSF_n slots with -10 (like G1/G2) instead of 0.0, which is a physically valid value and would let an MCCD or pre-fourth-moment SHAPES dict pass for a measurement. Drops the stray log_-i. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_0151inhpkat9cPPdMTG65Rwm --- log_-i | 1 - .../modules/make_cat_package/make_cat.py | 15 +- .../psfex_interp_package/psfex_interp.py | 137 +++++++++--------- tests/module/test_hsm_column_seams.py | 20 ++- tests/module/test_make_cat.py | 8 +- 5 files changed, 97 insertions(+), 84 deletions(-) delete mode 100644 log_-i diff --git a/log_-i b/log_-i deleted file mode 100644 index 989bf77a3..000000000 --- a/log_-i +++ /dev/null @@ -1 +0,0 @@ --i /tmp/pytest-of-cdaley/pytest-80/test_build_map_pole_guard_skip0/corners.txt -o /tmp/pytest-of-cdaley/pytest-80/test_build_map_pole_guard_skip0/cov.hsp -c 32 -n 1024 diff --git a/src/shapepipe/modules/make_cat_package/make_cat.py b/src/shapepipe/modules/make_cat_package/make_cat.py index 303384049..948709520 100644 --- a/src/shapepipe/modules/make_cat_package/make_cat.py +++ b/src/shapepipe/modules/make_cat_package/make_cat.py @@ -637,15 +637,17 @@ def _save_psf_data(self, galaxy_psf_path): n_obj = len(self._obj_id) # Per-epoch PSF shape columns copied from the producer's SHAPES dict: - # (column, empty-slot fill, dtype). HSM_T_PSF already holds T - # (sigma_to_T applied at the producer's _interpolate_me). + # (column, empty-slot fill, dtype). Fills are out of physical range + # so an unmeasured slot cannot pass for a measurement. HSM_T_PSF + # already holds T (sigma_to_T applied at the producer's + # _interpolate_me). psf_shape_cols = [ ("HSM_G1_PSF", -10.0, float), ("HSM_G2_PSF", -10.0, float), ("HSM_T_PSF", 0.0, float), ("HSM_FLAG_PSF", 1, "int16"), - ("HSM_M4_1_PSF", 0.0, float), - ("HSM_M4_2_PSF", 0.0, float), + ("HSM_M4_1_PSF", -10.0, float), + ("HSM_M4_2_PSF", -10.0, float), ("HSM_RHO4_PSF", -1.0, float), ] # Per-epoch exposure ID and CCD number, slot-aligned with HSM_*_PSF_n; @@ -679,8 +681,9 @@ def _save_psf_data(self, galaxy_psf_path): continue for name, fill, _ in psf_shape_cols: - # A SHAPES dict from an older producer may lack the - # fourth-moment keys; the slot then keeps its fill. + # A SHAPES dict without the fourth-moment keys (MCCD, or + # a producer predating them) leaves those slots at their + # out-of-range fill while FLAG still reports the fit. self._add2dict( f"{name}_{epoch + 1}", shapes.get(name, fill), idx ) diff --git a/src/shapepipe/modules/psfex_interp_package/psfex_interp.py b/src/shapepipe/modules/psfex_interp_package/psfex_interp.py index 63ab5a1bb..68019858c 100644 --- a/src/shapepipe/modules/psfex_interp_package/psfex_interp.py +++ b/src/shapepipe/modules/psfex_interp_package/psfex_interp.py @@ -181,6 +181,60 @@ def _m(p, q): return fourth_moment_1, fourth_moment_2, moms.moments_rho4 +# One HSM shape row per object, in this order; ``_hsm_columns`` names them. +_HSM_ROW = ("G1", "G2", "SIGMA", "FLAG", "M4_1", "M4_2", "RHO4") + + +def _hsm_row(image, moms, wcs=None): + """HSM shape row for one object, ordered as ``_HSM_ROW``.""" + return [ + moms.observed_shape.g1, + moms.observed_shape.g2, + moms.moments_sigma, + int(bool(moms.error_message)), + *_fourth_moments(image, moms, wcs), + ] + + +def _hsm_columns(shapes, obj): + """Named ``HSM_*_`` columns from HSM shape rows. + + @sc [label:schema] hsm-column-grammar + The only place the ``_HSM_ROW`` slot order meets the column grammar: + sigma is stored as ``T`` (``cs_util.size.sigma_to_T``), ``FLAG`` as int, + the fourth moments under ``M4_1``/``M4_2``/``RHO4``. Every writer in this + module (single-epoch, validation, multi-epoch) goes through here. + + Parameters + ---------- + shapes : numpy.ndarray + Shape rows, ``(n_obj, len(_HSM_ROW))`` or a single ``(len(_HSM_ROW),)`` + row; a single row yields Python scalars. + obj : str + Object token, ``"PSF"`` or ``"STAR"``. + + Returns + ------- + dict + Column name -> values. + + """ + shapes = np.asarray(shapes) + col = {name: shapes[..., i] for i, name in enumerate(_HSM_ROW)} + out = { + f"HSM_G1_{obj}": col["G1"], + f"HSM_G2_{obj}": col["G2"], + f"HSM_T_{obj}": cs_size.sigma_to_T(col["SIGMA"]), + f"HSM_FLAG_{obj}": col["FLAG"].astype(int), + f"HSM_M4_1_{obj}": col["M4_1"], + f"HSM_M4_2_{obj}": col["M4_2"], + f"HSM_RHO4_{obj}": col["RHO4"], + } + if shapes.ndim == 1: + out = {key: np.asarray(val).item() for key, val in out.items()} + return out + + class PSFExInterpolator(object): """The PSFEx Interpolator Class. @@ -477,8 +531,8 @@ def _get_psfshapes(self, wcs_list=None): Notes ----- In addition to the second-moment shape, the spin-2 fourth-moment - combinations and galsim's spin-0 ``moments_rho4`` are stored (columns 4, - 5, 6 of ``psf_shapes``). With a ``wcs_list`` these are computed in the + combinations and galsim's spin-0 ``moments_rho4`` are stored (the + ``M4_1``, ``M4_2``, ``RHO4`` slots of ``_HSM_ROW``). With a ``wcs_list`` these are computed in the world frame so they are invariant to the CCD orientation; see :func:`_fourth_moments`. @@ -502,13 +556,7 @@ def _get_psfshapes(self, wcs_list=None): self.psf_shapes = np.array( [ - [ - moms.observed_shape.g1, - moms.observed_shape.g2, - moms.moments_sigma, - int(bool(moms.error_message)), - *_fourth_moments(psf, moms, wcs), - ] + _hsm_row(psf, moms, wcs) for psf, moms, wcs in zip(self.interp_PSFs, psf_moms, wcs_list) ] ) @@ -528,13 +576,7 @@ def _write_output(self): if self._compute_shape: data = { "VIGNET": self.interp_PSFs, - "HSM_G1_PSF": self.psf_shapes[:, 0], - "HSM_G2_PSF": self.psf_shapes[:, 1], - "HSM_T_PSF": cs_size.sigma_to_T(self.psf_shapes[:, 2]), - "HSM_FLAG_PSF": self.psf_shapes[:, 3].astype(int), - "HSM_M4_1_PSF": self.psf_shapes[:, 4], - "HSM_M4_2_PSF": self.psf_shapes[:, 5], - "HSM_RHO4_PSF": self.psf_shapes[:, 6], + **_hsm_columns(self.psf_shapes, "PSF"), } else: data = {"VIGNET": self.interp_PSFs} @@ -619,7 +661,8 @@ def _get_starshapes(self, star_vign, wcs_list=None): Notes ----- As in :meth:`_get_psfshapes`, the spin-2 fourth-moment combinations and - galsim's ``moments_rho4`` are stored (columns 4, 5, 6 of ``star_shapes``), + galsim's ``moments_rho4`` are stored (the ``M4_1``, ``M4_2``, ``RHO4`` + slots of ``_HSM_ROW``), world-frame when a ``wcs_list`` is given. Masked pixels are zeroed before the fourth-moment sum so they do not contribute. @@ -651,13 +694,7 @@ def _get_starshapes(self, star_vign, wcs_list=None): self.star_shapes = np.array( [ - [ - moms.observed_shape.g1, - moms.observed_shape.g2, - moms.moments_sigma, - int(bool(moms.error_message)), - *_fourth_moments(np.where(mask == 1, 0.0, star), moms, wcs), - ] + _hsm_row(np.where(mask == 1, 0.0, star), moms, wcs) for star, mask, moms, wcs in zip( star_vign, masks, star_moms, wcs_list ) @@ -723,9 +760,9 @@ def _write_output_validation(self, star_dict, psfex_cat_dict): Save computed PSFs and stars to fits file. @sc [label:schema] psfex-validation-hsm-columns - The ``HSM_*`` columns written here are the exact set - ``MergeStarCatPSFEX.process`` reads without fallback; add or rename on - both sides together (``test_hsm_column_seams`` checks the two agree). + Writes ``_hsm_columns`` for both ``PSF`` and ``STAR`` — the exact + ``HSM_*`` set ``MergeStarCatPSFEX.process`` reads without fallback; + add or rename on both sides together (``test_hsm_column_seams``). Parameters ---------- @@ -742,20 +779,8 @@ def _write_output_validation(self, star_dict, psfex_cat_dict): ) data = { - "HSM_G1_PSF": self.psf_shapes[:, 0], - "HSM_G2_PSF": self.psf_shapes[:, 1], - "HSM_T_PSF": cs_size.sigma_to_T(self.psf_shapes[:, 2]), - "HSM_FLAG_PSF": self.psf_shapes[:, 3].astype(int), - "HSM_M4_1_PSF": self.psf_shapes[:, 4], - "HSM_M4_2_PSF": self.psf_shapes[:, 5], - "HSM_RHO4_PSF": self.psf_shapes[:, 6], - "HSM_G1_STAR": self.star_shapes[:, 0], - "HSM_G2_STAR": self.star_shapes[:, 1], - "HSM_T_STAR": cs_size.sigma_to_T(self.star_shapes[:, 2]), - "HSM_FLAG_STAR": self.star_shapes[:, 3].astype(int), - "HSM_M4_1_STAR": self.star_shapes[:, 4], - "HSM_M4_2_STAR": self.star_shapes[:, 5], - "HSM_RHO4_STAR": self.star_shapes[:, 6], + **_hsm_columns(self.psf_shapes, "PSF"), + **_hsm_columns(self.star_shapes, "STAR"), } data = {**data, **star_dict} @@ -808,9 +833,9 @@ def _interpolate_me(self): Interpolate PSFs for multi-epoch run. @sc [label:schema] psfex-me-shapes-columns - The per-epoch ``SHAPES`` dict must carry every column - ``make_cat._save_psf_data`` copies into ``HSM_*_PSF_n`` - (``test_hsm_column_seams`` checks the superset). + The per-epoch ``SHAPES`` dict is ``_hsm_columns(row, "PSF")`` and + must carry every column ``make_cat._save_psf_data`` copies into + ``HSM_*_PSF_n`` (``test_hsm_column_seams`` checks the superset). Raises ------ @@ -1020,31 +1045,9 @@ def _interpolate_me(self): "VIGNET" ] = final_list[j][1][where_res[0]] if self._compute_shape: - shape_dict = {} - shape_dict["HSM_G1_PSF"] = final_list[j][2][ - where_res[0] - ][0] - shape_dict["HSM_G2_PSF"] = final_list[j][2][ - where_res[0] - ][1] - shape_dict["HSM_T_PSF"] = cs_size.sigma_to_T( - final_list[j][2][where_res[0]][2] - ) - shape_dict["HSM_FLAG_PSF"] = final_list[j][2][ - where_res[0] - ][3] - shape_dict["HSM_M4_1_PSF"] = final_list[j][2][ - where_res[0] - ][4] - shape_dict["HSM_M4_2_PSF"] = final_list[j][2][ - where_res[0] - ][5] - shape_dict["HSM_RHO4_PSF"] = final_list[j][2][ - where_res[0] - ][6] output_dict[id_tmp][final_list[j][3][where_res[0]]][ "SHAPES" - ] = shape_dict + ] = _hsm_columns(final_list[j][2][where_res[0]], "PSF") counter += 1 if counter == 0: output_dict[id_tmp] = "empty" diff --git a/tests/module/test_hsm_column_seams.py b/tests/module/test_hsm_column_seams.py index c888d0a58..b05c52b69 100644 --- a/tests/module/test_hsm_column_seams.py +++ b/tests/module/test_hsm_column_seams.py @@ -2,21 +2,25 @@ Contracts ``psfex-validation-hsm-columns`` / ``psfex-starcat-columns-strict`` (psfex_interp ↔ merge_starcat) and ``psfex-me-shapes-columns`` / -``psf-epoch-slot-columns`` (psfex_interp ↔ make_cat): the column names are -string literals on both sides, read with no fallback, so the two sides are -compared as sets of ``HSM_*`` literals collected from each function's AST. +``psf-epoch-slot-columns`` (psfex_interp ↔ make_cat). The producer side is +``_hsm_columns`` (``hsm-column-grammar``), which every psfex_interp writer +goes through; the consumers read column names as string literals with no +fallback, collected from each function's AST. """ import ast import inspect import textwrap +import numpy as np + from shapepipe.modules.make_cat_package.make_cat import SaveCatalogue from shapepipe.modules.merge_starcat_package.merge_starcat import ( MergeStarCatPSFEX, ) from shapepipe.modules.psfex_interp_package.psfex_interp import ( - PSFExInterpolator, + _HSM_ROW, + _hsm_columns, ) @@ -46,9 +50,13 @@ def _hsm_literals(func, subscript_of=None): } +def _written(obj): + return set(_hsm_columns(np.zeros((2, len(_HSM_ROW))), obj)) + + def test_merge_starcat_reads_exactly_what_psfex_validation_writes(): """psfex-validation-hsm-columns == psfex-starcat-columns-strict.""" - written = _hsm_literals(PSFExInterpolator._write_output_validation) + written = _written("PSF") | _written("STAR") read = _hsm_literals(MergeStarCatPSFEX.process, subscript_of="data_j") assert written == read, { "written_not_read": sorted(written - read), @@ -58,6 +66,6 @@ def test_merge_starcat_reads_exactly_what_psfex_validation_writes(): def test_make_cat_epoch_columns_come_from_psfex_me_shapes(): """psf-epoch-slot-columns ⊆ psfex-me-shapes-columns.""" - produced = _hsm_literals(PSFExInterpolator._interpolate_me) + produced = _written("PSF") consumed = _hsm_literals(SaveCatalogue._save_psf_data) assert consumed <= produced, sorted(consumed - produced) diff --git a/tests/module/test_make_cat.py b/tests/module/test_make_cat.py index 85ae16de7..63c759d2b 100644 --- a/tests/module/test_make_cat.py +++ b/tests/module/test_make_cat.py @@ -459,8 +459,8 @@ def test_save_psf_data_carries_fourth_moments_per_epoch(tmp_path): """HSM_M4_1/M4_2/RHO4_PSF_n ride the same slots as HSM_G1_PSF_n. The fourth-moment columns psfex_interp writes into SHAPES (shapepipe#697) - land per epoch; a SHAPES dict without them (an older producer) leaves - the slot at its fill, as does an unused slot. + land per epoch; a SHAPES dict without them (MCCD, or an older producer) + leaves the slot at its out-of-range fill, as does an unused slot. """ galaxy_psf_path = tmp_path / "galaxy_psf.sqlite" per_obj = { @@ -477,7 +477,7 @@ def test_save_psf_data_carries_fourth_moments_per_epoch(tmp_path): npt.assert_allclose(out["HSM_M4_2_PSF_1"], [-0.22]) npt.assert_allclose(out["HSM_RHO4_PSF_1"], [2.05]) for n in (2, 3): - npt.assert_allclose(out[f"HSM_M4_1_PSF_{n}"], [0.0]) - npt.assert_allclose(out[f"HSM_M4_2_PSF_{n}"], [0.0]) + npt.assert_allclose(out[f"HSM_M4_1_PSF_{n}"], [-10.0]) + npt.assert_allclose(out[f"HSM_M4_2_PSF_{n}"], [-10.0]) npt.assert_allclose(out[f"HSM_RHO4_PSF_{n}"], [-1.0]) npt.assert_allclose(out["HSM_G1_PSF_2"], [0.03]) From 490bda7a0a97c60c3ce08521012f712b3f4f98eb Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 26 Sep 2026 02:35:06 +0200 Subject: [PATCH 8/8] test(psfex-interp): analytic oracle + spin-2 rotation law for the fourth moments The suite was null and invariance tests only, so a factor-1 or sign flip on M4_2, swapped whitened axes, or a wrong weight width all passed. The oracle predicts M4_1/M4_2 for a rotated two-component coaxial Gaussian from Isserlis fourth moments of each whitened, weighted component (agreement ~5e-8, tolerance 1e-6); the rotation test checks (M4_1 + i M4_2) picks up e^{2i dbeta}. Each of the four mutants now fails at least one of them. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_0151inhpkat9cPPdMTG65Rwm --- tests/module/test_psf_fourth_moments.py | 113 ++++++++++++++++++++++++ 1 file changed, 113 insertions(+) diff --git a/tests/module/test_psf_fourth_moments.py b/tests/module/test_psf_fourth_moments.py index 4782ca5f4..74e732ee7 100644 --- a/tests/module/test_psf_fourth_moments.py +++ b/tests/module/test_psf_fourth_moments.py @@ -120,6 +120,119 @@ def test_composite_has_nonzero_spin2(): npt.assert_allclose(row[M4_2], 0.0, atol=1e-4) +# --------------------------------------------------------------------------- +# Analytic oracle and spin-2 transformation -- pin the PSFHOME convention. +# --------------------------------------------------------------------------- + +# Pixel-frame composite (wcs=None): sigmas >= 3 px, so the pixel-centre +# sampling of each Gaussian reproduces its continuous moments essentially +# exactly and the oracle below needs no pixelisation correction. +_PIX_COMPONENTS = [(1.0, 3.0, 0.5), (1.0, 6.0, 0.1)] # (flux, sigma, e1) + + +def pixel_composite(beta_deg): + """Coaxial two-Gaussian composite in pixel units, rotated by ``beta_deg``.""" + prof = galsim.Add( + [ + galsim.Gaussian(flux=f, sigma=s).shear(e1=e1) + for f, s, e1 in _PIX_COMPONENTS + ] + ) + return prof.rotate(beta_deg * galsim.degrees) + + +def pixel_measure(prof): + """Draw ``prof`` at unit pixel scale; return ``(moms, _fourth_moments)``.""" + image = prof.drawImage(nx=_STAMP, ny=_STAMP, scale=1.0, method="no_pixel") + moms = galsim.hsm.FindAdaptiveMom(galsim.Image(image.array)) + return moms, _fourth_moments(image.array, moms) + + +def oracle_fourth_moments(moms, beta_deg): + """Analytic weighted fourth moments of :func:`pixel_composite`. + + Component k is a Gaussian of flux F_k and covariance + ``C_k = s_k**2 R A_k A_k R^T`` (A_k the unit-determinant shear matrix, + R the rotation). Whitening by ``S = sqrtm(inv(M))``, with ``M`` rebuilt + independently from HSM's ``(sigma, e1, e2)`` as ``sigma**2 A A^T``, gives + ``C_k' = S C_k S``. The weight ``exp(-r^T r / 2)`` turns it into a Gaussian + of covariance ``Sig_k = inv(I + inv(C_k'))`` and integrated amplitude + ``a_k = F_k sqrt(det Sig_k / det C_k')``. Isserlis then gives + ``E[u^4] = 3 Sig_uu^2``, ``E[v^4] = 3 Sig_vv^2``, + ``E[u^3 v] = 3 Sig_uu Sig_uv``, ``E[u v^3] = 3 Sig_vv Sig_uv``, and each + ``M_pq`` is the ``a_k``-weighted mean over components. + """ + shape = moms.observed_shape + a_obj = galsim.Shear(e1=shape.e1, e2=shape.e2).getMatrix() + M = moms.moments_sigma**2 * a_obj @ a_obj.T + w, V = np.linalg.eigh(M) + S = V @ np.diag(w**-0.5) @ V.T + + th = np.deg2rad(beta_deg) + R = np.array([[np.cos(th), -np.sin(th)], [np.sin(th), np.cos(th)]]) + + amp, m40, m04, m31, m13 = [], [], [], [], [] + for flux, sigma, e1 in _PIX_COMPONENTS: + a_k = galsim.Shear(e1=e1).getMatrix() + C = S @ (sigma**2 * R @ a_k @ a_k.T @ R.T) @ S + Sig = np.linalg.inv(np.eye(2) + np.linalg.inv(C)) + amp.append(flux * np.sqrt(np.linalg.det(Sig) / np.linalg.det(C))) + suu, svv, suv = Sig[0, 0], Sig[1, 1], Sig[0, 1] + m40.append(3 * suu**2) + m04.append(3 * svv**2) + m31.append(3 * suu * suv) + m13.append(3 * svv * suv) + + amp = np.array(amp) + + def mean(x): + return np.sum(amp * np.array(x)) / np.sum(amp) + + return mean(m40) - mean(m04), 2 * (mean(m13) + mean(m31)) + + +@pytest.mark.parametrize("beta", [30.0, -55.0]) +def test_fourth_moments_match_analytic_oracle(beta): + """Pixel-frame composite: ``_fourth_moments`` equals the Isserlis oracle. + + Pins the PSFHOME convention absolutely -- whitening by the object's own + adaptive second moments, weight ``exp(-r^2/2)`` in whitened coordinates, + flux-normalised ``M4_1 = M40 - M04`` and ``M4_2 = 2 (M13 + M31)``. The + rotation puts signal in both spin-2 components. + + Tolerance ``rtol=1e-6``: with sigmas >= 3 px the pixel-centre sampling + error is exponentially small and the stamp holds > 8 sigma, so the only + residual is HSM's centroid/moment convergence. Measured agreement is + ~5e-8 relative on both components (beta = 30, -55 deg). + """ + moms, (m4_1, m4_2, _) = pixel_measure(pixel_composite(beta)) + o4_1, o4_2 = oracle_fourth_moments(moms, beta) + + assert abs(o4_1) > 1e-2 and abs(o4_2) > 1e-2 # both components excited + npt.assert_allclose(m4_1, o4_1, rtol=1e-6) + npt.assert_allclose(m4_2, o4_2, rtol=1e-6) + + +@pytest.mark.parametrize("delta", [25.0, 70.0, -40.0]) +def test_fourth_moments_rotate_as_spin2(delta): + """Rotating the object by ``delta`` multiplies ``M4_1 + i M4_2`` by + ``exp(2 i delta)``. + + Independent of the oracle, this fixes the relative sign and factor + between the two components and their spin: a flipped sign, a lost factor + of 2, or swapped whitened axes each break the transformation law. + """ + beta = 30.0 + _, (a1, a2, _) = pixel_measure(pixel_composite(beta)) + _, (b1, b2, _) = pixel_measure(pixel_composite(beta + delta)) + + z_ref = complex(a1, a2) + assert abs(a1) > 1e-2 and abs(a2) > 1e-2 # both components excited + expected = np.exp(2j * np.deg2rad(delta)) * z_ref + npt.assert_allclose(b1, expected.real, rtol=1e-3, atol=1e-3 * abs(z_ref)) + npt.assert_allclose(b2, expected.imag, rtol=1e-3, atol=1e-3 * abs(z_ref)) + + # --------------------------------------------------------------------------- # Frame invariance -- the science guarantee. # ---------------------------------------------------------------------------