Add optional high order moments to the HSM output file.

[rmjarvis]

@ztq1996 has made a convincing case that we should look at the high-order moments of the PSF to look for potential biases in the PSF models leading to biased galaxy shapes.

This would be especially easy for many use cases if we add these moment measurements to the standard hsm output catalog we make in Piff. My proposal is to add a fourth_order=True option to the HSMCatalog statistic type, which would turn on these measurements.

To make these more directly useful to galaxy shape error diagnostics, I propose that we output these in terms of the spin-0, spin-2 and spin-4 combinations:

spin-0: M_40 + 2M_22 + M_04.
spin-2: (M_40 − M_04) + i(2M_31 + 2M_13)
spin-4: (M_40 - 6M_22 + M_04) + i(4M_31 - 4M_13)

The spin-4 one probably isn't too useful (e.g. there is no spin-4 quantity in TreeCorr to do proper correlations with it), but the other two definitely are.

TQ calls the first two rho and M^(4)_PSF, respectively, which I don't love as names for these. They are kind of the 4-th order T and e quantities, so maybe T4 and e4 would be more accessible to most users? I'm open to suggestions here.

We also need some name for the spin-4 quantity (TQ didn't name it AFAICT).

There is already a function in Piff to calculate the high-order moments, so I think this will be easy to add. Just call that function if fourth_order is set to true and add those columns to the output file. So the implementation should be pretty straightforward.

Refs:
https://ui.adsabs.harvard.edu/abs/2022MNRAS.510.1978Z/abstract
https://ui.adsabs.harvard.edu/abs/2022arXiv221203257Z/abstract

[ztq1996]

Hi @rmjarvis -

Thanks for pointing this out. I feel like we have touched on this a few times verbally but never written something down.

I just took a look at piff.util.calculate_moments(). It looks pretty good to me for setting up higher moments measurement. I just have a few questions and comments:

1. I have a hard time lining up the index of M_{pq} in piff with its mathematical meaning: for example, for second moments, M11=trace, M20=e1, M02=e2. I don't think we need to change how the indices are done here, as that could break some other code, but maybe we need some explanation for the indices.
2. Here du and dv are coordinated from the centroid. So the higher moments will be correlated to second moments (larger PSF = larger higher radial moments, for example). This is slightly different from the definition that I have been using, which erases the second moment's information. I think both definitions are probably fine for capturing additive biases in shear, but I need to look at their difference with a certain dataset to say that for sure.

As in the name, I agree T4 and e4 are much easier to understand. I would still suggest adding a parenthesis, i.e., T^(4) and e^(4) to avoid confusion with the power of the second moments. I am open to suggestions too.

[rmjarvis]

I agree that the naming scheme for the M's in calculate_moments is a bit non-standard. I think the notation came from @aaronroodman. But it's trying to hew closer to the "spin" combinations. Here is the logic of it.

For p>=q:
M_pq = Re(<(u+iv)^p (u-iv)^q>)
M_qp = Im(<(u+iv)^p (u-iv)^q>)

So M_20 is aka e_1. M_02 is e_2. M11 is T.
At fourth order, M_22 is the spin-0 quantity. M_31 + iM_13 is spin-2. And M_40 + iM04 is spin-4.

[aaronroodman]

Credit Chris Davis for this scheme

…

_________________________________________________________________________ Prof. Aaron Roodman Deputy Director Vera C. Rubin Observatory Dept. of Particle Physics & Astrophysics Kavli Institute for Particle Astrophysics & Cosmology SLAC National Accelerator Laboratory Stanford University SLAC National Accelerator Laboratory E-mail: ***@***.******@***.***> 2575 Sand Hill Rd. Phone: 650-926-2705 MS 29 Menlo Park, CA 94025 URL: http://www.slac.stanford.edu/~roodman

_________________________________________________________________________ From: Mike Jarvis ***@***.***> Date: Monday, February 6, 2023 at 7:46 AM To: rmjarvis/Piff ***@***.***> Cc: Roodman, Aaron ***@***.***>, Mention ***@***.***> Subject: Re: [rmjarvis/Piff] Add optional high order moments to the HSM output file. (Issue #142) I agree that the naming scheme for the M's in calculate_moments is a bit non-standard. I think the notation came from @aaronroodman<https://github.com/aaronroodman>. But it's trying to hew closer to the "spin" combinations. Here is the logic of it. For p>=q: M_pq = Re(<(u+iv)^p (u-iv)^q>) M_qp = Im(<(u+iv)^p (u-iv)^q>) So M_20 is aka e_1. M_02 is e_2. M11 is T. At fourth order, M_22 is the spin-0 quantity. M_31 + iM_13 is spin-2. And M_40 + iM04 is spin-4. — Reply to this email directly, view it on GitHub<#142 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ABARV46NHQNCO6I36677P5TWWEMEXANCNFSM6AAAAAAUBS4ETY>. You are receiving this because you were mentioned.Message ID: ***@***.***>

[rmjarvis]

PR here: #144