Skip to content

Commit

Permalink
add normalized moments and moment errors by name
Browse files Browse the repository at this point in the history
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
  • Loading branch information
esheldon committed Dec 13, 2024
1 parent f2a6958 commit a6bf29d
Show file tree
Hide file tree
Showing 2 changed files with 45 additions and 0 deletions.
38 changes: 38 additions & 0 deletions ngmix/moments.py
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
7 changes: 7 additions & 0 deletions ngmix/tests/test_gmix.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down

0 comments on commit a6bf29d

Please sign in to comment.