DM-41648: Write a Plugin for computing higher-order moments

[arunkannawadi]

Setting up a base branch (+PR) to receive contributions from Tianqing Zhang

[arunkannawadi]

Pinging @timj to weigh in on the COPYRIGHT file.

[timj]

LSST Corp copyrights were handed to AURA so these should be combined with AURA entry.

[jmeyers314]

Not done yet, but wanted you to see comments so far since I was also starring at _calculate_higher_order_moments...

[jmeyers314]

I know this is just a blackening, but personally, I find the following (black-acceptable) format nicer to read:

shearType = pexConfig.ChoiceField[str](
    doc=(
        "The desired method of PSF correction using GalSim. The first three options return an e-type "
        "distortion, whereas the last option returns a g-type shear."
    ),
)

Your call though.

[jmeyers314]

Maybe just Image of source?

[jmeyers314]

I think this implies that the normalized x coordinate lies either on the unnormalized major or minor axis. Can you clarify which (assuming I got that first part right)?

[arunkannawadi]

It does, and I don't think it matters whether it lies along the major or minor axis due to the normalization.

[jmeyers314]

If you swapped x from major -> minor, doesn't that mean that you also swap Mpq to Mqp ? That's what I'm trying to clarify, that the moments are in some sense M_major_minor. Or am I missing something?

[arunkannawadi]

Oh hmm, good point. I think x is along the major axis, but this makes me question if I should think of these numbers are coordinate-independent or not. I had naively assumed it was. If we have to transform them from pixel-coordinates to sky-coordinates, we need to remember to transform these accordingly then.

[jmeyers314]

If you're going to leave commented-out code, then I think you need to attach an explanation as to why it's commented out.

[arunkannawadi]

Oh, these should be removed actually.

[jmeyers314]

Just saw the use_linear_algebra bit appear, so seems you're staring at this too :).

Aside from the 2x2 specialization (which I think should be the only code path BTW), I got a tad bit more performance by caching the power arrays. Something like:

    std_x_pow = {}
    std_y_pow = {}

    pmin = min(p for p, q in pqlist)
    qmin = min(q for p, q in pqlist)
    if pmin == 0:
        std_x_pow[0] = 1.0
    elif pmin == 1:
        std_x_pow[1] = std_x
    else:
        std_x_pow[pmin] = std_x**pmin

    if qmin == 0:
        std_y_pow[0] = 1.0
    elif qmin == 1:
        std_y_pow[1] = std_y
    else:
        std_y_pow[qmin] = std_y**qmin

    for p in range(pmin, max(p for p, q in pqlist)+1):
        std_x_pow[p+1] = std_x_pow[p] * std_x
    for q in range(qmin, max(q for p, q in pqlist)+1):
        std_y_pow[q+1] = std_y_pow[q] * std_y

    results = {}
    for p, q in pqlist:
        tmp = std_x_pow[p] * std_y_pow[q]
        tmp *= image_weight
        results[(p, q)] = np.sum(tmp) / normalization

[arunkannawadi]

Good idea, and I hope this does well even on objects with large footprints.

pmin and qmin are always going to be 0 given how we generate those tuples. So I am going to remove the if-elif-else conditions.

[jmeyers314]

Actually, it occurs to me that np.vander may be your friend here (but I didn't try it).

[arunkannawadi]

I think I'd like to have the option of having the more expensive code path, given that that was what TQ originally used in the previous HSC data, so I consider that validated. I'd like to have an option to resort to the validated code if something doesn't look right. I could add another unit test to check that both code paths result in consistent answers.

[arunkannawadi]

It looks like np.vander is only for 1-d arrays. So I went with your code suggestion, minus the condition checking.

[arunkannawadi]

Summary of changes since you last reviewed:

• We accumulate the powers of standard coordinates instead of computing on the fly
• Instead of making two copies of arrays to make badpix, there is only one copying now.
• Added a test to check the validate method of the configs.
• Added a test that checks that the codepaths with either values for use_linear_algebra produces the same result, even with an asymmetric M matrix.
• Fixed other documentation related comments.

[erykoff]

I think I found the bug

[erykoff]

This needs to be a copy of the sub-image.

[jmeyers314]

Might be clearer if the copy is inside _calculate_higher_order_moments, but yeah, good catch.

[arunkannawadi]

Copying inside _calculate_higher_order_moments would make unnecessary copies for psfImage. I fixed it to make a deep copy on a need-to basis.

[arunkannawadi]

Alternatively, I could make the modification to the weight image and not have to make a deep copy ever.

[jmeyers314]

Right. I just meant we don't want anyone to have to remember that _calculate_higher_order_moments modifies the image. Modifying the weight is a great solution.

[arunkannawadi]

The commit I just pushed looks to me like the cleanest solution - multiply the weight and image with no regards for the badpix and mask this product as configured. It avoids making expensive deep copies, and avoids any division of pixel values (as I currently have) which may become unstable if the pixel values are close to 0.

[erykoff]

Looks much better!

[jmeyers314]

I guess these could be lists now. Probably doesn't matter in practice...

[jmeyers314]

LGTM. Looking forward to seeing some higher order moments analysis in the near future!