fourth-order moments - #698
martinkilbinger wants to merge 8 commits into
Conversation
…psfex_interp module.
|
@martinkilbinger I modified the psfex script to compute the 4th order moments. I checked it runs on the star vignets but not on the interpolated PSF. The easiest is probably to rerun this part of the code like you did already and check if it saves correctly the 4th order moment information. I will copy and paste the script in the corresponding functions for the MCCD equivalent. |
|
Reopened as reference. This branch is ~435 commits behind — Fable on behalf of Cail |
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 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_0151inhpkat9cPPdMTG65Rwm
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 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_0151inhpkat9cPPdMTG65Rwm
|
Closing in favour of #859, which measures the fourth moments directly in world coordinates on top of #812. Here they are measured in pixels and converted afterwards in — Fable on behalf of Cail |
* 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 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01K98geWHTuT62whqjKocfCm * 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 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01K98geWHTuT62whqjKocfCm * 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 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01K98geWHTuT62whqjKocfCm * 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 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01K98geWHTuT62whqjKocfCm * 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 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_0151inhpkat9cPPdMTG65Rwm * test: HSM column seam checks + @sc contracts for the fourth-moment chain 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 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_0151inhpkat9cPPdMTG65Rwm * 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 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_0151inhpkat9cPPdMTG65Rwm * 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 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_0151inhpkat9cPPdMTG65Rwm --------- Co-authored-by: Claude Fable 5 <noreply@anthropic.com>
Summary
Implementing fourth-order moments.
Solves #697
Reviewer Checklist
developbranch