Skip to content

fourth-order moments - #698

Closed
martinkilbinger wants to merge 8 commits into
developfrom
moment4
Closed

martinkilbinger wants to merge 8 commits into
developfrom
moment4

Conversation

@martinkilbinger

@martinkilbinger martinkilbinger commented May 12, 2025 •

Copy link
Copy Markdown
Contributor

Summary

Implementing fourth-order moments.

Solves #697

Reviewer Checklist

Reviewers should tick the following boxes before approving and merging the PR.

  • The PR targets the develop branch
  • The PR is assigned to the developer
  • The PR has appropriate labels
  • The PR is included in appropriate projects and/or milestones
  • The PR includes a clear description of the proposed changes
  • If the PR addresses an open issue the description includes "closes #"
  • The code and documentation style match the current standards
  • Documentation has been added/updated consistently with the code
  • All CI tests are passing
  • API docs have been built and checked at least once (if relevant)
  • All changed files have been checked and comments provided to the developer
  • All of the reviewer's comments have been satisfactorily addressed by the developer

@sachaguer

Copy link
Copy Markdown

@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.

@martinkilbinger
martinkilbinger deleted the moment4 branch January 16, 2026 07:21
@cailmdaley cailmdaley reopened this Jul 21, 2026
@cailmdaley

Copy link
Copy Markdown
Contributor

Reopened as reference. This branch is ~435 commits behind develop and edits the same two functions PR #812 rewrites for sky-coordinate HSM measurement. The plan (see #697) is to re-implement on top of #812 rather than rebase this branch: spin-2 fourth-moment combinations in the sky frame, plus galsim's moments_rho4.

— Fable on behalf of Cail

cailmdaley added a commit that referenced this pull request Sep 26, 2026
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
cailmdaley added a commit that referenced this pull request Sep 26, 2026
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
@cailmdaley

cailmdaley commented Sep 26, 2026 •

Copy link
Copy Markdown
Contributor

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 convert_psf_pix2world.py as a shear. Under the WCS the whitened frame rotates by an object-dependent angle rather than shearing, so measuring in the sky frame is the exact route. Nothing from this branch is missing there.

— Fable on behalf of Cail

@cailmdaley cailmdaley closed this Sep 26, 2026
cailmdaley added a commit that referenced this pull request Sep 28, 2026
* 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>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants