PSF fourth-order moments in the sky frame (spin-2 + rho4)

[cailmdaley]

Closes #697. Supersedes #698.

Summary

The PSFEx chain now measures the fourth-order moments of the PSF model and of stars, alongside the second-order shape (g1, g2, T). They are measured directly in the sky frame, as the second moments have been since #812. Three new quantities per object:

Column	Definition	Spin
HSM_M4_1_{PSF,STAR}	M40 − M04	2
HSM_M4_2_{PSF,STAR}	2(M13 + M31)	2
HSM_RHO4_{PSF,STAR}	galsim moments_rho4 (radial kurtosis; 2 for a Gaussian)	0

Here M_pq are the whitened, Gaussian-weighted moments of PSFHOME (Zhang et al. 2022, 2023); see Design.

Problem

• PSF-model errors in the fourth moments also bias shear, and ShapePipe measured only second moments. Diagnosing this needs:
  • star-vs-model fourth moments at star positions;
  • PSF fourth moments at galaxy positions, for leakage regression.
• fourth-order moments #698 ported the PSFHOME computation (Sacha) and plumbed it through the validation chain (Martin). It measured in the pixel frame and converted to the sky afterwards in convert_psf_pix2world.py, treating (M4_1, M4_2) like a shear. HSM: measure PSF/star shapes directly in sky coords (use_sky_coords) #812 retired that script by measuring second moments directly in sky coordinates. This PR does the same for the fourth moments, which also makes the transformation exact (see Design).
• fourth-order moments #698 predates the ngmix v2.0: one column grammar across estimators (PSF source-of-truth, single size, E↔G) #761 column grammar.
• galsim already computes rho4 inside FindAdaptiveMom; ShapePipe discarded it.

Scope

• In: psfex_interp (single-epoch, validation and multi-epoch outputs), merge_starcat, and make_cat's per-epoch PSF slots.
• Out: MCCD, which writes its own HSM columns (PSFEx is the fiducial PSF model).

Design

Measurement

For each PSF or star stamp (_fourth_moments in psfex_interp.py):

1. HSM adaptive moments with use_sky_coords=True give the shape, σ and centroid in world coordinates. The second-moment matrix M is built from them.
2. The local WCS Jacobian maps the pixel grid to world offsets from that centroid.
3. The offsets are whitened with the symmetric M^{-1/2} (a matched elliptical Gaussian becomes a unit circle) and weighted with a unit Gaussian.
4. The centred p+q=4 moments give M4_1 and M4_2. rho4 comes from the same HSM call.

Why measure in the sky frame: under a WCS map A, the whitened coordinates change by Q = M'^{-1/2} A M^{1/2}. That is a pure rotation or flip, and its angle depends on the object's own M. So the pixel-frame values differ from the sky values by an object-dependent rotation, not by a shear. Building M and the grid in world coordinates gives the sky values directly.

Data flow

flowchart LR
  P["psfex_interp<br/>_hsm_columns"] -->|validation catalogue| M["merge_starcat"] --> S[("merged star catalogue<br/>HSM_*_PSF, HSM_*_STAR")]
  P -->|"multi-epoch SHAPES (SAVE_PSF_DATA)"| C["make_cat"] --> F[("final galaxy catalogue<br/>HSM_*_PSF_n per-epoch slots")]
  P -.->|"star_cat_merge (PR 879)"| S2[("full_starcat_{campaign}.hdf5")]
  F -.->|"final_cat_merge (PR 879)"| F2[("final_cat_{campaign}.hdf5")]

Loading

Solid edges carry the new columns in this PR. Dashed edges are follow-ups.

• One column mapper. _HSM_ROW names the shape-row slots once (G1, G2, SIGMA, FLAG, M4_1, M4_2, RHO4). _hsm_columns(shapes, obj) maps them to column names for all three psfex_interp writers, so no writer indexes a shape row by hand.
• merge_starcat reads the six new columns by name, with no fallback. This avoids fourth-order moments #698's bug of filling M_4_PSF_2 from M_4_PSF_1.
• make_cat carries HSM_M4_1/M4_2/RHO4_PSF_n in the same per-epoch slots as HSM_G1_PSF_n, from a single column table.

Missing values

• psfex_interp outputs and the merged star catalogue:
  • On HSM failure: M4_1 = M4_2 = 0, RHO4 = −1 and HSM_FLAG_* = 1. This matches galsim's g = 0 on failure.
  • Select on the flag.
• make_cat per-epoch slots:
  • A slot with no PSFEx fourth-moment measurement holds M4 = −10 and RHO4 = −1, like HSM_G1_PSF_n. That covers an empty slot, an HSM failure, MCCD, or an older producer.
  • For MCCD or an older producer, HSM_FLAG_PSF_n is 0, so select on RHO4 > 0.

Verification

• Against theory (test_psf_fourth_moments.py):
  • An analytic oracle (Isserlis' theorem on a rotated two-component Gaussian) pins the signs, the factor 2 and the weight width.
  • The spin-2 rotation law pins the relative sign of the two components.
  • Invariance under WCS rotation and under world-frame shifts pins the frame and the centroid.
  • Masking and the failure fills are covered too.
• Across modules (test_hsm_column_seams.py): merge_starcat must read exactly the HSM_* set that _hsm_columns writes, and make_cat a subset of it.
• CI is green.

Cost and compatibility

• Cost: the fourth-moment step adds ~1 ms per stamp, about 3× the HSM fit. The multi-epoch run pays it for every galaxy-epoch.
• Compatibility: merge_starcat fails on validation catalogues written before this PR. Regenerate them, as after ngmix v2.0: one column grammar across estimators (PSF source-of-truth, single size, E↔G) #761.

Follow-ups

The #879 campaign merges pick fixed column lists, so the new columns stop there for now:

• Star side: add the six HSM_M4_* / HSM_RHO4_* columns to COLUMNS in workflow/scripts/merge_star_cat.py.
• Galaxy side: the tile config doesn't set SAVE_PSF_DATA, so no HSM_*_PSF_n slots are written. final_cat.param can list them only once make_cat writes a fixed number of slots per tile. Today it writes max(N_EPOCH)+1, which varies by tile.

— Claude (Opus) on behalf of Cail

🤖 Generated with Claude Code

[sachaguer]

This PR looks good to me, great job @cailmdaley !
I did not click approve by mistake but consider it done ;)