diff --git a/src/shapepipe/modules/make_cat_package/make_cat.py b/src/shapepipe/modules/make_cat_package/make_cat.py index 062ac5719..948709520 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,30 @@ 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). 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", -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; # -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 +667,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 +677,15 @@ 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 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 + ) galaxy_psf_cat.close() diff --git a/src/shapepipe/modules/merge_starcat_package/merge_starcat.py b/src/shapepipe/modules/merge_starcat_package/merge_starcat.py index 7bc0cb76b..0d724d210 100644 --- a/src/shapepipe/modules/merge_starcat_package/merge_starcat.py +++ b/src/shapepipe/modules/merge_starcat_package/merge_starcat.py @@ -556,10 +556,17 @@ 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 = [], [], [] + 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 +595,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 +649,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, diff --git a/src/shapepipe/modules/psfex_interp_package/psfex_interp.py b/src/shapepipe/modules/psfex_interp_package/psfex_interp.py index b6da6623f..68019858c 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,178 @@ 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. + + 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 + + +# 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. @@ -355,6 +528,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 (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`. + """ if import_fail: raise ImportError("Galsim is required to get shapes information") @@ -371,16 +552,12 @@ 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( [ - [ - moms.observed_shape.g1, - moms.observed_shape.g2, - moms.moments_sigma, - int(bool(moms.error_message)), - ] - for moms in psf_moms + _hsm_row(psf, moms, wcs) + for psf, moms, wcs in zip(self.interp_PSFs, psf_moms, wcs_list) ] ) @@ -399,10 +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_columns(self.psf_shapes, "PSF"), } else: data = {"VIGNET": self.interp_PSFs} @@ -484,6 +658,14 @@ 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 (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. + """ if import_fail: raise ImportError("Galsim is required to get shapes information") @@ -508,16 +690,14 @@ 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( [ - [ - moms.observed_shape.g1, - moms.observed_shape.g2, - moms.moments_sigma, - int(bool(moms.error_message)), - ] - for moms in star_moms + _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 + ) ] ) @@ -579,6 +759,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 + 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 ---------- star_dict : dict @@ -594,14 +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_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_columns(self.psf_shapes, "PSF"), + **_hsm_columns(self.star_shapes, "STAR"), } data = {**data, **star_dict} @@ -653,6 +832,11 @@ def _interpolate_me(self): Interpolate PSFs for multi-epoch run. + @sc [label:schema] psfex-me-shapes-columns + 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 ------ KeyError @@ -861,22 +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] 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 new file mode 100644 index 000000000..b05c52b69 --- /dev/null +++ b/tests/module/test_hsm_column_seams.py @@ -0,0 +1,71 @@ +"""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 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 ( + _HSM_ROW, + _hsm_columns, +) + + +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 _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 = _written("PSF") | _written("STAR") + 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 = _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 34a1958e8..63c759d2b 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 (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 = { + 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}"], [-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]) diff --git a/tests/module/test_psf_fourth_moments.py b/tests/module/test_psf_fourth_moments.py new file mode 100644 index 000000000..74e732ee7 --- /dev/null +++ b/tests/module/test_psf_fourth_moments.py @@ -0,0 +1,368 @@ +"""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) + + +# --------------------------------------------------------------------------- +# 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. +# --------------------------------------------------------------------------- + + +@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-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]) +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-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(): + """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) + + +# --------------------------------------------------------------------------- +# 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. +# --------------------------------------------------------------------------- + + +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, + ) 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) ]