From a6bf29d5af4cae0130457d4256d8266b7910a54d Mon Sep 17 00:00:00 2001 From: Erin Sheldon Date: Fri, 13 Dec 2024 10:23:15 -0500 Subject: [PATCH] add normalized moments and moment errors by name Add normalized moments by name, such as M20 which is the M20 sums divided by the flux sums Also add uncertainties Higher order moments will be included only if present Since the only failure mode for these is flux < 0, we don't store flags. We instead store the number as nan --- ngmix/moments.py | 38 ++++++++++++++++++++++++++++++++++++++ ngmix/tests/test_gmix.py | 7 +++++++ 2 files changed, 45 insertions(+) diff --git a/ngmix/moments.py b/ngmix/moments.py index b0f0d980..c1e028f5 100644 --- a/ngmix/moments.py +++ b/ngmix/moments.py @@ -534,9 +534,47 @@ def make_mom_result(sums, sums_cov, sums_norm=None): res["flux_flagstr"] = ngmix.flags.get_flags_str(res["flux_flags"]) res["T_flagstr"] = ngmix.flags.get_flags_str(res["T_flags"]) + _add_moments_by_name(res) + return res +def _add_moments_by_name(res): + sums = res['sums'] + sums_cov = res['sums_cov'] + + mf_ind = MOMENTS_NAME_MAP["MF"] + fsum = sums[mf_ind] + fsum_err = np.sqrt(sums_cov[mf_ind, mf_ind]) + + # add in named sums normalized by flux sum (weight * image).sum() + # we don't store flags or errors for these + mkeys = list(MOMENTS_NAME_MAP.keys()) + for name in mkeys: + ind = MOMENTS_NAME_MAP[name] + if ind > sums.size-1: + continue + + err_name = f'{name}_err' + + if name in ['MF', 'M00']: + res[name] = fsum + res[err_name] = fsum_err + else: + if fsum > 0: + res[name] = sums[ind] / fsum + res[err_name] = get_ratio_error( + sums[ind], + sums[mf_ind], + sums_cov[ind, ind], + sums_cov[mf_ind, mf_ind], + sums_cov[ind, mf_ind], + ) + else: + res[name] = np.nan + res[err_name] = np.nan + + def regularize_mom_shapes(res, fwhm_reg): """Apply regularization to the shapes computed from moments sums. diff --git a/ngmix/tests/test_gmix.py b/ngmix/tests/test_gmix.py index 021c25ef..585095f1 100644 --- a/ngmix/tests/test_gmix.py +++ b/ngmix/tests/test_gmix.py @@ -524,10 +524,17 @@ def test_higher_order_smoke(do_higher): for name, ind in MOMENTS_NAME_MAP.items(): # make sure the index is valid res['sums'][ind] + assert name in res and np.isfinite(res[name]) + err_name = f'{name}_err' + assert err_name in res and np.isfinite(res[err_name]) else: assert res['sums'].shape == (6, ) assert res['sums_cov'].shape == (6, 6) + for name, ind in MOMENTS_NAME_MAP.items(): + if ind > 5: + assert name not in res + def test_higher_order(): rng = np.random.RandomState(seed=35)