From 03624e9588dc83a3b5ce91fd849fe903b7195d69 Mon Sep 17 00:00:00 2001 From: beckermr Date: Tue, 19 Jul 2022 11:41:45 -0500 Subject: [PATCH 01/14] BUG fix bug in wmom moms --- CHANGES.md | 4 ++++ ngmix/gmix/gmix.py | 14 ++++++++++---- 2 files changed, 14 insertions(+), 4 deletions(-) diff --git a/CHANGES.md b/CHANGES.md index 2f136baa..3a0c2f98 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -1,5 +1,9 @@ ## v2.1.0 +### bug fixes + + - Fixed bug in computing `mom` and `mom_cov` fields for weighted Gaussian moments. + ### new features - Added `fwhm_smooth` keyword to pre-PSF moments routines to allow for extra diff --git a/ngmix/gmix/gmix.py b/ngmix/gmix/gmix.py index 30782710..bb34e465 100644 --- a/ngmix/gmix/gmix.py +++ b/ngmix/gmix/gmix.py @@ -1239,15 +1239,21 @@ def get_weighted_moments_stats(ares): res = {} for n in ares.dtype.names: - if n == "sums": - res[n] = ares[n].copy() - elif n == "sums_cov": + if n in ["sums", "sums_cov"]: res[n] = ares[n].copy() else: res[n] = ares[n] + scales = np.array([res["sums"][5]] * 5 + [1.0]) + res["mom"] = res["sums"].copy() + res["mom"] /= scales + res["mom_cov"] = res["sums_cov"].copy() + for i in range(6): + for j in range(6): + res["mom_cov"][i, j] = res["mom_cov"][i, j]/scales[i]/scales[j] + res.update( - moments.make_mom_result(res["sums"], res["sums_cov"]) + moments.make_mom_result(res["mom"], res["mom_cov"]) ) return res From 86bde0f8d08dc0309aa4dfbbddfa5e5ac5e3c7a2 Mon Sep 17 00:00:00 2001 From: beckermr Date: Tue, 19 Jul 2022 11:50:13 -0500 Subject: [PATCH 02/14] TST lower precision due to change in moments --- mdet_tests/test_mdet_regression.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/mdet_tests/test_mdet_regression.py b/mdet_tests/test_mdet_regression.py index e48e565b..f3f83c7e 100644 --- a/mdet_tests/test_mdet_regression.py +++ b/mdet_tests/test_mdet_regression.py @@ -232,7 +232,8 @@ def test_mdet_regression(fname, write=False): else: old_data = fitsio.read(fname) for col in old_data.dtype.names: - if ('psf_' in col or 'ratio' in col) and '1.3.9' in fname: + if ('psf_' in col or 'ratio' in col) and ( + '1.3.9' in fname or '2.0' in fname): rtol = 1.0e-4 atol = 1.0e-4 else: From b69335c9bd3cbbcc8c8f2f7c8e12f46b35f09f6f Mon Sep 17 00:00:00 2001 From: beckermr Date: Tue, 19 Jul 2022 12:37:22 -0500 Subject: [PATCH 03/14] ENH record moment normalization --- mdet_tests/test_mdet_regression.py | 6 ++++-- ngmix/gmix/gmix.py | 10 +--------- ngmix/moments.py | 5 ++++- ngmix/prepsfmom.py | 9 ++++++--- ngmix/tests/test_prepsfmom.py | 8 ++++---- 5 files changed, 19 insertions(+), 19 deletions(-) diff --git a/mdet_tests/test_mdet_regression.py b/mdet_tests/test_mdet_regression.py index f3f83c7e..d202b523 100644 --- a/mdet_tests/test_mdet_regression.py +++ b/mdet_tests/test_mdet_regression.py @@ -232,8 +232,10 @@ def test_mdet_regression(fname, write=False): else: old_data = fitsio.read(fname) for col in old_data.dtype.names: - if ('psf_' in col or 'ratio' in col) and ( - '1.3.9' in fname or '2.0' in fname): + if ( + ('psf_' in col or 'ratio' in col) + and ('1.3.9' in fname or '2.0' in fname) + ): rtol = 1.0e-4 atol = 1.0e-4 else: diff --git a/ngmix/gmix/gmix.py b/ngmix/gmix/gmix.py index bb34e465..dc5f70ac 100644 --- a/ngmix/gmix/gmix.py +++ b/ngmix/gmix/gmix.py @@ -1244,16 +1244,8 @@ def get_weighted_moments_stats(ares): else: res[n] = ares[n] - scales = np.array([res["sums"][5]] * 5 + [1.0]) - res["mom"] = res["sums"].copy() - res["mom"] /= scales - res["mom_cov"] = res["sums_cov"].copy() - for i in range(6): - for j in range(6): - res["mom_cov"][i, j] = res["mom_cov"][i, j]/scales[i]/scales[j] - res.update( - moments.make_mom_result(res["mom"], res["mom_cov"]) + moments.make_mom_result(res["sums"], res["sums_cov"], res["wsum"]) ) return res diff --git a/ngmix/moments.py b/ngmix/moments.py index c5120d59..b8e98ba4 100644 --- a/ngmix/moments.py +++ b/ngmix/moments.py @@ -349,7 +349,7 @@ def g2mom(g1, g2, T): return Irr, Irc, Icc -def make_mom_result(mom, mom_cov): +def make_mom_result(mom, mom_cov, mom_norm): """Make a fitting results dict from a set of moments. Parameters @@ -358,6 +358,8 @@ def make_mom_result(mom, mom_cov): The array of moments in the order [Mv, Mu, M1, M2, MT, MF]. mom_cov : np.ndarray The array of moment covariances. + mom_norm : float + The sum of the moment weight function itself. Returns ------- @@ -389,6 +391,7 @@ def make_mom_result(mom, mom_cov): res["flux"] = mom[mf_ind] res["mom"] = mom res["mom_cov"] = mom_cov + res["mom_norm"] = mom_norm res["flux_flags"] = 0 res["flux_flagstr"] = "" res["T_flags"] = 0 diff --git a/ngmix/prepsfmom.py b/ngmix/prepsfmom.py index 9ca4361d..51f0699f 100644 --- a/ngmix/prepsfmom.py +++ b/ngmix/prepsfmom.py @@ -160,11 +160,11 @@ def _meas(self, obs, psf_obs, return_kernels): tot_var = np.sum(1.0 / obs.weight[msk]) # run the actual measurements and return - mom, mom_cov = _measure_moments_fft( + mom, mom_cov, mom_norm = _measure_moments_fft( kim, kpsf_im, tot_var, eff_pad_factor, kernels, im_row - psf_im_row, im_col - psf_im_col, ) - res = make_mom_result(mom, mom_cov) + res = make_mom_result(mom, mom_cov, mom_norm) if res['flags'] != 0: logger.debug("pre-psf moments failed: %s" % res['flagstr']) @@ -291,6 +291,7 @@ def _measure_moments_fft( fkp = kernels["fkp"] fkc = kernels["fkc"] + mom_norm = np.sum(kernels["fk00"].real) * df2 mf = np.sum((kim * fkf).real) * df2 mr = np.sum((kim * fkr).real) * df2 mp = np.sum((kim * fkp).real) * df2 @@ -324,7 +325,7 @@ def _measure_moments_fft( mom = np.array([np.nan, np.nan, mp, mc, mr, mf]) - return mom, m_cov + return mom, m_cov, mom_norm @njit @@ -580,6 +581,7 @@ def _ksigma_kernels( fkc=fkc, msk=msk, nrm=nrm, + fk00=fkf[0, 0], ) @@ -679,6 +681,7 @@ def _gauss_kernels( fkc=fkc, msk=msk, nrm=nrm, + fk00=fkf[0, 0], ) diff --git a/ngmix/tests/test_prepsfmom.py b/ngmix/tests/test_prepsfmom.py index 99718cfd..a7b6bea6 100644 --- a/ngmix/tests/test_prepsfmom.py +++ b/ngmix/tests/test_prepsfmom.py @@ -762,7 +762,7 @@ def test_moments_make_mom_result_flags(): for i in range(2, 6): _mom_cov = mom_cov.copy() _mom_cov[i, i] = -1 - res = make_mom_result(mom, _mom_cov) + res = make_mom_result(mom, _mom_cov, 1) assert (res["flags"] & ngmix.flags.NONPOS_VAR) != 0 assert ngmix.flags.NAME_MAP[ngmix.flags.NONPOS_VAR] in res["flagstr"] if i == 5: @@ -782,7 +782,7 @@ def test_moments_make_mom_result_flags(): # neg flux _mom = mom.copy() _mom[5] = -1 - res = make_mom_result(_mom, mom_cov) + res = make_mom_result(_mom, mom_cov, 1) assert (res["flags"] & ngmix.flags.NONPOS_FLUX) != 0 assert ngmix.flags.NAME_MAP[ngmix.flags.NONPOS_FLUX] in res["flagstr"] assert res["flux_flags"] == 0 @@ -793,7 +793,7 @@ def test_moments_make_mom_result_flags(): # neg T _mom = mom.copy() _mom[4] = -1 - res = make_mom_result(_mom, mom_cov) + res = make_mom_result(_mom, mom_cov, 1) assert (res["flags"] & ngmix.flags.NONPOS_SIZE) != 0 assert ngmix.flags.NAME_MAP[ngmix.flags.NONPOS_SIZE] in res["flagstr"] assert res["flux_flags"] == 0 @@ -806,7 +806,7 @@ def test_moments_make_mom_result_flags(): _mom_cov = mom_cov.copy() _mom_cov[4, i] = np.nan _mom_cov[i, 4] = np.nan - res = make_mom_result(mom, _mom_cov) + res = make_mom_result(mom, _mom_cov, 1) assert (res["flags"] & ngmix.flags.NONPOS_SHAPE_VAR) != 0 assert ngmix.flags.NAME_MAP[ngmix.flags.NONPOS_SHAPE_VAR] in res["flagstr"] assert res["flux_flags"] == 0 From ae5e3528136c39b507858ed6abac03615dc7de41 Mon Sep 17 00:00:00 2001 From: beckermr Date: Tue, 19 Jul 2022 12:38:13 -0500 Subject: [PATCH 04/14] DOC more accurate changelog --- CHANGES.md | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/CHANGES.md b/CHANGES.md index 3a0c2f98..b6e0de4e 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -1,14 +1,11 @@ ## v2.1.0 -### bug fixes - - - Fixed bug in computing `mom` and `mom_cov` fields for weighted Gaussian moments. - ### new features - Added `fwhm_smooth` keyword to pre-PSF moments routines to allow for extra smoothing of the profile before the moments are measured. - Added caching of FFTs in metacal and pre-PSF moment rountines. + - Added moments weight function normalization. ## v2.0.6 From 680aa53fb4fca09f9bf858292f9100cac1388fe1 Mon Sep 17 00:00:00 2001 From: beckermr Date: Tue, 19 Jul 2022 12:43:32 -0500 Subject: [PATCH 05/14] BUG remove breaking API change --- ngmix/moments.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/ngmix/moments.py b/ngmix/moments.py index b8e98ba4..93d5d1f1 100644 --- a/ngmix/moments.py +++ b/ngmix/moments.py @@ -349,7 +349,7 @@ def g2mom(g1, g2, T): return Irr, Irc, Icc -def make_mom_result(mom, mom_cov, mom_norm): +def make_mom_result(mom, mom_cov, mom_norm=None): """Make a fitting results dict from a set of moments. Parameters @@ -358,8 +358,9 @@ def make_mom_result(mom, mom_cov, mom_norm): The array of moments in the order [Mv, Mu, M1, M2, MT, MF]. mom_cov : np.ndarray The array of moment covariances. - mom_norm : float - The sum of the moment weight function itself. + mom_norm : float, optional + The sum of the moment weight function itself. This is added to the output data. + The default of None puts in NaN. Returns ------- @@ -391,7 +392,7 @@ def make_mom_result(mom, mom_cov, mom_norm): res["flux"] = mom[mf_ind] res["mom"] = mom res["mom_cov"] = mom_cov - res["mom_norm"] = mom_norm + res["mom_norm"] = mom_norm if mom_norm is not None else np.nan res["flux_flags"] = 0 res["flux_flagstr"] = "" res["T_flags"] = 0 From 23a07f6fa979c22570b0ac237d8fafe70a78f31b Mon Sep 17 00:00:00 2001 From: beckermr Date: Tue, 19 Jul 2022 12:53:56 -0500 Subject: [PATCH 06/14] BUG only one dim here --- ngmix/prepsfmom.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ngmix/prepsfmom.py b/ngmix/prepsfmom.py index 51f0699f..e8db79d9 100644 --- a/ngmix/prepsfmom.py +++ b/ngmix/prepsfmom.py @@ -581,7 +581,7 @@ def _ksigma_kernels( fkc=fkc, msk=msk, nrm=nrm, - fk00=fkf[0, 0], + fk00=fkf[0], ) @@ -681,7 +681,7 @@ def _gauss_kernels( fkc=fkc, msk=msk, nrm=nrm, - fk00=fkf[0, 0], + fk00=fkf[0], ) From 7b0635b442595a8c3f7fc231dec03794a43f3701 Mon Sep 17 00:00:00 2001 From: beckermr Date: Tue, 19 Jul 2022 12:56:52 -0500 Subject: [PATCH 07/14] BUG use actual value --- ngmix/prepsfmom.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ngmix/prepsfmom.py b/ngmix/prepsfmom.py index e8db79d9..b66673a0 100644 --- a/ngmix/prepsfmom.py +++ b/ngmix/prepsfmom.py @@ -581,7 +581,7 @@ def _ksigma_kernels( fkc=fkc, msk=msk, nrm=nrm, - fk00=fkf[0], + fk00=knrm, ) @@ -681,7 +681,7 @@ def _gauss_kernels( fkc=fkc, msk=msk, nrm=nrm, - fk00=fkf[0], + fk00=knrm, ) From 570ec869e44f8d04c5d53ab37904b79671769564 Mon Sep 17 00:00:00 2001 From: "Matthew R. Becker" Date: Tue, 19 Jul 2022 13:03:09 -0500 Subject: [PATCH 08/14] Update ngmix/tests/test_prepsfmom.py --- ngmix/tests/test_prepsfmom.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ngmix/tests/test_prepsfmom.py b/ngmix/tests/test_prepsfmom.py index a7b6bea6..ccc3bcd0 100644 --- a/ngmix/tests/test_prepsfmom.py +++ b/ngmix/tests/test_prepsfmom.py @@ -806,7 +806,7 @@ def test_moments_make_mom_result_flags(): _mom_cov = mom_cov.copy() _mom_cov[4, i] = np.nan _mom_cov[i, 4] = np.nan - res = make_mom_result(mom, _mom_cov, 1) + res = make_mom_result(mom, _mom_cov, mom_norm=1) assert (res["flags"] & ngmix.flags.NONPOS_SHAPE_VAR) != 0 assert ngmix.flags.NAME_MAP[ngmix.flags.NONPOS_SHAPE_VAR] in res["flagstr"] assert res["flux_flags"] == 0 From 7b57811e922ef4dc4e6e90a4cedcd3648e778965 Mon Sep 17 00:00:00 2001 From: "Matthew R. Becker" Date: Tue, 19 Jul 2022 13:03:35 -0500 Subject: [PATCH 09/14] Apply suggestions from code review --- ngmix/gmix/gmix.py | 2 +- ngmix/prepsfmom.py | 2 +- ngmix/tests/test_prepsfmom.py | 6 +++--- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/ngmix/gmix/gmix.py b/ngmix/gmix/gmix.py index dc5f70ac..eb2a92c3 100644 --- a/ngmix/gmix/gmix.py +++ b/ngmix/gmix/gmix.py @@ -1245,7 +1245,7 @@ def get_weighted_moments_stats(ares): res[n] = ares[n] res.update( - moments.make_mom_result(res["sums"], res["sums_cov"], res["wsum"]) + moments.make_mom_result(res["sums"], res["sums_cov"], mom_norm=res["wsum"]) ) return res diff --git a/ngmix/prepsfmom.py b/ngmix/prepsfmom.py index b66673a0..8d8bbb69 100644 --- a/ngmix/prepsfmom.py +++ b/ngmix/prepsfmom.py @@ -164,7 +164,7 @@ def _meas(self, obs, psf_obs, return_kernels): kim, kpsf_im, tot_var, eff_pad_factor, kernels, im_row - psf_im_row, im_col - psf_im_col, ) - res = make_mom_result(mom, mom_cov, mom_norm) + res = make_mom_result(mom, mom_cov, mom_norm=mom_norm) if res['flags'] != 0: logger.debug("pre-psf moments failed: %s" % res['flagstr']) diff --git a/ngmix/tests/test_prepsfmom.py b/ngmix/tests/test_prepsfmom.py index ccc3bcd0..c701f6cf 100644 --- a/ngmix/tests/test_prepsfmom.py +++ b/ngmix/tests/test_prepsfmom.py @@ -762,7 +762,7 @@ def test_moments_make_mom_result_flags(): for i in range(2, 6): _mom_cov = mom_cov.copy() _mom_cov[i, i] = -1 - res = make_mom_result(mom, _mom_cov, 1) + res = make_mom_result(mom, _mom_cov, mom_norm=1) assert (res["flags"] & ngmix.flags.NONPOS_VAR) != 0 assert ngmix.flags.NAME_MAP[ngmix.flags.NONPOS_VAR] in res["flagstr"] if i == 5: @@ -782,7 +782,7 @@ def test_moments_make_mom_result_flags(): # neg flux _mom = mom.copy() _mom[5] = -1 - res = make_mom_result(_mom, mom_cov, 1) + res = make_mom_result(_mom, mom_cov, mom_norm=1) assert (res["flags"] & ngmix.flags.NONPOS_FLUX) != 0 assert ngmix.flags.NAME_MAP[ngmix.flags.NONPOS_FLUX] in res["flagstr"] assert res["flux_flags"] == 0 @@ -793,7 +793,7 @@ def test_moments_make_mom_result_flags(): # neg T _mom = mom.copy() _mom[4] = -1 - res = make_mom_result(_mom, mom_cov, 1) + res = make_mom_result(_mom, mom_cov, mom_norm=1) assert (res["flags"] & ngmix.flags.NONPOS_SIZE) != 0 assert ngmix.flags.NAME_MAP[ngmix.flags.NONPOS_SIZE] in res["flagstr"] assert res["flux_flags"] == 0 From d94ba5b39c03da2c280813738f014fd9d8367be5 Mon Sep 17 00:00:00 2001 From: beckermr Date: Thu, 21 Jul 2022 06:03:24 -0500 Subject: [PATCH 10/14] BUG fix prepsf moment norms --- ngmix/gmix/gmix.py | 2 +- ngmix/prepsfmom.py | 4 ++-- ngmix/tests/test_prepsfmom.py | 36 ++++++++++++++++++++++++++++++++++- 3 files changed, 38 insertions(+), 4 deletions(-) diff --git a/ngmix/gmix/gmix.py b/ngmix/gmix/gmix.py index dc5f70ac..e0066fc6 100644 --- a/ngmix/gmix/gmix.py +++ b/ngmix/gmix/gmix.py @@ -1245,7 +1245,7 @@ def get_weighted_moments_stats(ares): res[n] = ares[n] res.update( - moments.make_mom_result(res["sums"], res["sums_cov"], res["wsum"]) + moments.make_mom_result(res["sums"].copy(), res["sums_cov"].copy(), res["wsum"]) ) return res diff --git a/ngmix/prepsfmom.py b/ngmix/prepsfmom.py index b66673a0..31b4c10a 100644 --- a/ngmix/prepsfmom.py +++ b/ngmix/prepsfmom.py @@ -164,7 +164,7 @@ def _meas(self, obs, psf_obs, return_kernels): kim, kpsf_im, tot_var, eff_pad_factor, kernels, im_row - psf_im_row, im_col - psf_im_col, ) - res = make_mom_result(mom, mom_cov, mom_norm) + res = make_mom_result(mom, mom_cov, mom_norm=mom_norm) if res['flags'] != 0: logger.debug("pre-psf moments failed: %s" % res['flagstr']) @@ -291,7 +291,7 @@ def _measure_moments_fft( fkp = kernels["fkp"] fkc = kernels["fkc"] - mom_norm = np.sum(kernels["fk00"].real) * df2 + mom_norm = kernels["fk00"] mf = np.sum((kim * fkf).real) * df2 mr = np.sum((kim * fkr).real) * df2 mp = np.sum((kim * fkp).real) * df2 diff --git a/ngmix/tests/test_prepsfmom.py b/ngmix/tests/test_prepsfmom.py index a7b6bea6..cf62de44 100644 --- a/ngmix/tests/test_prepsfmom.py +++ b/ngmix/tests/test_prepsfmom.py @@ -909,12 +909,46 @@ def test_prepsfmom_gauss_true_flux( assert_allclose(flux_true, 400, atol=0, rtol=5e-3) +@pytest.mark.parametrize('pixel_scale', [0.25, 0.125]) +@pytest.mark.parametrize('image_size', [107]) +@pytest.mark.parametrize('pad_factor', [3.5, 4]) +@pytest.mark.parametrize('mom_fwhm', [2, 2.5]) +@pytest.mark.parametrize('cls', [PGaussMom, KSigmaMom]) +def test_prepsfmom_mom_norm( + pad_factor, image_size, pixel_scale, mom_fwhm, cls, +): + rng = np.random.RandomState(seed=100) + + cen = (image_size - 1)/2 + gs_wcs = galsim.ShearWCS( + pixel_scale, galsim.Shear(g1=-0.1, g2=0.06)).jacobian() + scale = np.sqrt(gs_wcs.pixelArea()) + shift = rng.uniform(low=-scale/2, high=scale/2, size=2) + xy = gs_wcs.toImage(galsim.PositionD(shift)) + + jac = Jacobian( + y=cen + xy.y, x=cen + xy.x, + dudx=gs_wcs.dudx, dudy=gs_wcs.dudy, + dvdx=gs_wcs.dvdx, dvdy=gs_wcs.dvdy) + + obs = Observation( + image=np.ones((image_size, image_size)), + jacobian=jac, + ) + res = cls( + fwhm=mom_fwhm, pad_factor=pad_factor, + ).go( + obs=obs, no_psf=True, + ) + assert_allclose(res["mom_norm"], res["flux"], atol=0, rtol=2e-4) + + @pytest.mark.parametrize('pixel_scale', [0.25, 0.125]) @pytest.mark.parametrize('fwhm', [2, 0.5]) @pytest.mark.parametrize('image_size', [107]) @pytest.mark.parametrize('pad_factor', [3.5, 4]) @pytest.mark.parametrize('mom_fwhm', [2, 2.5]) -def test_prepsfmom_comp_to_gaussmom( +def test_prepsfmom_comp_to_gaussmom_simple( pad_factor, image_size, fwhm, pixel_scale, mom_fwhm ): rng = np.random.RandomState(seed=100) From 03e6cbdaa54d23425b3d715f3b9c0e48b2870d77 Mon Sep 17 00:00:00 2001 From: beckermr Date: Thu, 21 Jul 2022 06:24:37 -0500 Subject: [PATCH 11/14] BUG use area norm for all gauss mom things --- ngmix/gaussmom.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/ngmix/gaussmom.py b/ngmix/gaussmom.py index 9ebd663e..70faccce 100644 --- a/ngmix/gaussmom.py +++ b/ngmix/gaussmom.py @@ -58,6 +58,13 @@ def _measure_moments(self, obs): res['flux'] *= fac res['flux_err'] *= fac res['pars'][5] *= fac + res["mom"] *= fac + res["mom_cov"] *= fac**2 + res["mom_norm"] *= fac + res["sums"] *= fac + res["sums_cov"] *= fac**2 + res["wsum"] *= fac + res["mom_err"] *= fac return res def _set_mompars(self): From dfcc4fe5e17db16c6234f059d6598781242a61fd Mon Sep 17 00:00:00 2001 From: beckermr Date: Thu, 21 Jul 2022 06:26:41 -0500 Subject: [PATCH 12/14] DOC add changelog entry --- CHANGES.md | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/CHANGES.md b/CHANGES.md index 2f136baa..2b178ac3 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -1,5 +1,9 @@ ## v2.1.0 +### bug fixes + + - Fixed 1/area normalization for all moments fields in `GaussMom` results. + ### new features - Added `fwhm_smooth` keyword to pre-PSF moments routines to allow for extra From d087bba34604d3c853f9b15f31847b26e21071f5 Mon Sep 17 00:00:00 2001 From: beckermr Date: Thu, 21 Jul 2022 12:47:16 -0500 Subject: [PATCH 13/14] REF rename mom to sums everywhere --- ngmix/gaussmom.py | 6 +-- ngmix/moments.py | 84 +++++++++++++++++------------------ ngmix/prepsfmom.py | 2 +- ngmix/tests/test_prepsfmom.py | 38 ++++++++-------- 4 files changed, 64 insertions(+), 66 deletions(-) diff --git a/ngmix/gaussmom.py b/ngmix/gaussmom.py index 70faccce..75c5d62c 100644 --- a/ngmix/gaussmom.py +++ b/ngmix/gaussmom.py @@ -58,13 +58,11 @@ def _measure_moments(self, obs): res['flux'] *= fac res['flux_err'] *= fac res['pars'][5] *= fac - res["mom"] *= fac - res["mom_cov"] *= fac**2 - res["mom_norm"] *= fac res["sums"] *= fac res["sums_cov"] *= fac**2 + res["sums_norm"] *= fac res["wsum"] *= fac - res["mom_err"] *= fac + res["sums_err"] *= fac return res def _set_mompars(self): diff --git a/ngmix/moments.py b/ngmix/moments.py index 93d5d1f1..81ca936c 100644 --- a/ngmix/moments.py +++ b/ngmix/moments.py @@ -349,16 +349,16 @@ def g2mom(g1, g2, T): return Irr, Irc, Icc -def make_mom_result(mom, mom_cov, mom_norm=None): - """Make a fitting results dict from a set of moments. +def make_mom_result(sums, sums_cov, sums_norm=None): + """Make a fitting results dict from a set of unnormalized moments. Parameters ---------- - mom : np.ndarray - The array of moments in the order [Mv, Mu, M1, M2, MT, MF]. - mom_cov : np.ndarray - The array of moment covariances. - mom_norm : float, optional + sums : np.ndarray + The array of unnormalized moments in the order [Mv, Mu, M1, M2, MT, MF]. + sums_cov : np.ndarray + The array of unnormalized moment covariances. + sums_norm : float, optional The sum of the moment weight function itself. This is added to the output data. The default of None puts in NaN. @@ -368,12 +368,12 @@ def make_mom_result(mom, mom_cov, mom_norm=None): A dictionary of results. """ # for safety... - if len(mom) != 6: + if len(sums) != 6: raise ValueError( - "You must pass exactly 6 moments in the order [Mv, Mu, M1, M2, MT, MF] " - "for ngmix.moments.make_mom_result." + "You must pass exactly 6 unnormalized moments in the order " + "[Mv, Mu, M1, M2, MT, MF] for ngmix.moments.make_mom_result." ) - if mom_cov.shape != (6, 6): + if sums_cov.shape != (6, 6): raise ValueError( "You must pass a 6x6 matrix for ngmix.moments.make_mom_result." ) @@ -389,10 +389,10 @@ def make_mom_result(mom, mom_cov, mom_norm=None): res = {} res["flags"] = 0 res["flagstr"] = "" - res["flux"] = mom[mf_ind] - res["mom"] = mom - res["mom_cov"] = mom_cov - res["mom_norm"] = mom_norm if mom_norm is not None else np.nan + res["flux"] = sums[mf_ind] + res["sums"] = sums + res["sums_cov"] = sums_cov + res["sums_norm"] = sums_norm if sums_norm is not None else np.nan res["flux_flags"] = 0 res["flux_flagstr"] = "" res["T_flags"] = 0 @@ -408,25 +408,25 @@ def make_mom_result(mom, mom_cov, mom_norm=None): res["e"] = np.array([np.nan, np.nan]) res["e_err"] = np.array([np.nan, np.nan]) res["e_cov"] = np.diag([np.nan, np.nan]) - res["mom_err"] = np.array([np.nan, np.nan, np.nan, np.nan, np.nan, np.nan]) + res["sums_err"] = np.array([np.nan, np.nan, np.nan, np.nan, np.nan, np.nan]) # handle flux-only - if mom_cov[mf_ind, mf_ind] > 0: - res["flux_err"] = np.sqrt(mom_cov[mf_ind, mf_ind]) + if sums_cov[mf_ind, mf_ind] > 0: + res["flux_err"] = np.sqrt(sums_cov[mf_ind, mf_ind]) res["s2n"] = res["flux"] / res["flux_err"] else: res["flux_flags"] |= ngmix.flags.NONPOS_VAR # handle flux+T only - if mom_cov[mf_ind, mf_ind] > 0 and mom_cov[mt_ind, mt_ind] > 0: - if mom[mf_ind] > 0: - res["T"] = mom[mt_ind] / mom[mf_ind] + if sums_cov[mf_ind, mf_ind] > 0 and sums_cov[mt_ind, mt_ind] > 0: + if sums[mf_ind] > 0: + res["T"] = sums[mt_ind] / sums[mf_ind] res["T_err"] = get_ratio_error( - mom[mt_ind], - mom[mf_ind], - mom_cov[mt_ind, mt_ind], - mom_cov[mf_ind, mf_ind], - mom_cov[mt_ind, mf_ind], + sums[mt_ind], + sums[mf_ind], + sums_cov[mt_ind, mt_ind], + sums_cov[mf_ind, mf_ind], + sums_cov[mt_ind, mf_ind], ) else: # flux <= 0.0 @@ -435,21 +435,21 @@ def make_mom_result(mom, mom_cov, mom_norm=None): res["T_flags"] |= ngmix.flags.NONPOS_VAR # now handle full flags - if np.all(np.diagonal(mom_cov) > 0): - res["mom_err"] = np.sqrt(np.diagonal(mom_cov)) + if np.all(np.diagonal(sums_cov) > 0): + res["sums_err"] = np.sqrt(np.diagonal(sums_cov)) else: res["flags"] |= ngmix.flags.NONPOS_VAR if res["flags"] == 0: if res["flux"] > 0: if res["T"] > 0: - res["e1"] = mom[m1_ind] / mom[mt_ind] - res["e2"] = mom[m2_ind] / mom[mt_ind] + res["e1"] = sums[m1_ind] / sums[mt_ind] + res["e2"] = sums[m2_ind] / sums[mt_ind] res["e"] = np.array([res["e1"], res["e2"]]) res["pars"] = np.array([ - mom[mv_ind], - mom[mu_ind], + sums[mv_ind], + sums[mu_ind], res["e1"], res["e2"], res["T"], @@ -458,18 +458,18 @@ def make_mom_result(mom, mom_cov, mom_norm=None): e_err = np.zeros(2) e_err[0] = get_ratio_error( - mom[m1_ind], - mom[mt_ind], - mom_cov[m1_ind, m1_ind], - mom_cov[mt_ind, mt_ind], - mom_cov[m1_ind, mt_ind], + sums[m1_ind], + sums[mt_ind], + sums_cov[m1_ind, m1_ind], + sums_cov[mt_ind, mt_ind], + sums_cov[m1_ind, mt_ind], ) e_err[1] = get_ratio_error( - mom[m2_ind], - mom[mt_ind], - mom_cov[m2_ind, m2_ind], - mom_cov[mt_ind, mt_ind], - mom_cov[m2_ind, mt_ind] + sums[m2_ind], + sums[mt_ind], + sums_cov[m2_ind, m2_ind], + sums_cov[mt_ind, mt_ind], + sums_cov[m2_ind, mt_ind] ) if np.all(np.isfinite(e_err)): res["e_err"] = e_err diff --git a/ngmix/prepsfmom.py b/ngmix/prepsfmom.py index 31b4c10a..135cd7a5 100644 --- a/ngmix/prepsfmom.py +++ b/ngmix/prepsfmom.py @@ -164,7 +164,7 @@ def _meas(self, obs, psf_obs, return_kernels): kim, kpsf_im, tot_var, eff_pad_factor, kernels, im_row - psf_im_row, im_col - psf_im_col, ) - res = make_mom_result(mom, mom_cov, mom_norm=mom_norm) + res = make_mom_result(mom, mom_cov, sums_norm=mom_norm) if res['flags'] != 0: logger.debug("pre-psf moments failed: %s" % res['flagstr']) diff --git a/ngmix/tests/test_prepsfmom.py b/ngmix/tests/test_prepsfmom.py index 31e12767..bf2f494c 100644 --- a/ngmix/tests/test_prepsfmom.py +++ b/ngmix/tests/test_prepsfmom.py @@ -371,8 +371,8 @@ def test_prepsfmom_gauss( _report_info("T", res["T"], T_true, np.mean(res["T_err"])) _report_info("g1", res["e"][:, 0], g1_true, np.mean(res["e_err"][0])) _report_info("g2", res["e"][:, 1], g2_true, np.mean(res["e_err"][1])) - mom_cov = np.cov(res["mom"].T) - print("mom cov ratio:\n", np.mean(res["mom_cov"], axis=0)/mom_cov, flush=True) + mom_cov = np.cov(res["sums"].T) + print("sums cov ratio:\n", np.mean(res["sums_cov"], axis=0)/mom_cov, flush=True) assert_allclose( np.abs(np.mean(res["flux"]) - flux_true)/np.mean(res["flux_err"]), 0, @@ -496,10 +496,10 @@ def test_prepsfmom_mn_cov_psf( _report_info("T", res["T"], T_true, np.mean(res["T_err"])) _report_info("g1", res["e"][:, 0], g1_true, np.mean(res["e_err"][0])) _report_info("g2", res["e"][:, 1], g2_true, np.mean(res["e_err"][1])) - mom_cov = np.cov(res["mom"].T) - print("mom cov ratio:\n", np.mean(res["mom_cov"], axis=0)/mom_cov, flush=True) - print("mom cov meas:\n", mom_cov, flush=True) - print("mom cov pred:\n", np.mean(res["mom_cov"], axis=0), flush=True) + mom_cov = np.cov(res["sums"].T) + print("sums cov ratio:\n", np.mean(res["sums_cov"], axis=0)/mom_cov, flush=True) + print("sums cov meas:\n", mom_cov, flush=True) + print("sums cov pred:\n", np.mean(res["sums_cov"], axis=0), flush=True) assert_allclose(np.mean(res["flux"]), flux_true, atol=0, rtol=1e-2) assert_allclose(np.mean(res["T"]), T_true, atol=0, rtol=1e-2) @@ -515,14 +515,14 @@ def test_prepsfmom_mn_cov_psf( assert_allclose( mom_cov[2:, 2:], - np.mean(res["mom_cov"][:, 2:, 2:], axis=0), + np.mean(res["sums_cov"][:, 2:, 2:], axis=0), atol=2.5e-1, rtol=0, ) assert_allclose( np.diagonal(mom_cov[2:, 2:]), - np.diagonal(np.mean(res["mom_cov"][:, 2:, 2:], axis=0)), + np.diagonal(np.mean(res["sums_cov"][:, 2:, 2:], axis=0)), atol=0, rtol=2e-2, ) @@ -722,10 +722,10 @@ def test_prepsfmom_mn_cov_nopsf( _report_info("T", res["T"], T_true, np.mean(res["T_err"])) _report_info("g1", res["e"][:, 0], g1_true, np.mean(res["e_err"][0])) _report_info("g2", res["e"][:, 1], g2_true, np.mean(res["e_err"][1])) - mom_cov = np.cov(res["mom"].T) - print("mom cov ratio:\n", np.mean(res["mom_cov"], axis=0)/mom_cov, flush=True) - print("mom cov meas:\n", mom_cov, flush=True) - print("mom cov pred:\n", np.mean(res["mom_cov"], axis=0), flush=True) + mom_cov = np.cov(res["sums"].T) + print("sums cov ratio:\n", np.mean(res["sums_cov"], axis=0)/mom_cov, flush=True) + print("sums cov meas:\n", mom_cov, flush=True) + print("sums cov pred:\n", np.mean(res["sums_cov"], axis=0), flush=True) assert_allclose(np.mean(res["flux"]), flux_true, atol=0, rtol=1e-2) assert_allclose(np.mean(res["T"]), T_true, atol=0, rtol=1e-2) @@ -741,14 +741,14 @@ def test_prepsfmom_mn_cov_nopsf( assert_allclose( mom_cov[2:, 2:], - np.mean(res["mom_cov"][:, 2:, 2:], axis=0), + np.mean(res["sums_cov"][:, 2:, 2:], axis=0), atol=2.5e-1, rtol=0, ) assert_allclose( np.diagonal(mom_cov[2:, 2:]), - np.diagonal(np.mean(res["mom_cov"][:, 2:, 2:], axis=0)), + np.diagonal(np.mean(res["sums_cov"][:, 2:, 2:], axis=0)), atol=0, rtol=2e-2, ) @@ -762,7 +762,7 @@ def test_moments_make_mom_result_flags(): for i in range(2, 6): _mom_cov = mom_cov.copy() _mom_cov[i, i] = -1 - res = make_mom_result(mom, _mom_cov, mom_norm=1) + res = make_mom_result(mom, _mom_cov, sums_norm=1) assert (res["flags"] & ngmix.flags.NONPOS_VAR) != 0 assert ngmix.flags.NAME_MAP[ngmix.flags.NONPOS_VAR] in res["flagstr"] if i == 5: @@ -782,7 +782,7 @@ def test_moments_make_mom_result_flags(): # neg flux _mom = mom.copy() _mom[5] = -1 - res = make_mom_result(_mom, mom_cov, mom_norm=1) + res = make_mom_result(_mom, mom_cov, sums_norm=1) assert (res["flags"] & ngmix.flags.NONPOS_FLUX) != 0 assert ngmix.flags.NAME_MAP[ngmix.flags.NONPOS_FLUX] in res["flagstr"] assert res["flux_flags"] == 0 @@ -793,7 +793,7 @@ def test_moments_make_mom_result_flags(): # neg T _mom = mom.copy() _mom[4] = -1 - res = make_mom_result(_mom, mom_cov, mom_norm=1) + res = make_mom_result(_mom, mom_cov, sums_norm=1) assert (res["flags"] & ngmix.flags.NONPOS_SIZE) != 0 assert ngmix.flags.NAME_MAP[ngmix.flags.NONPOS_SIZE] in res["flagstr"] assert res["flux_flags"] == 0 @@ -806,7 +806,7 @@ def test_moments_make_mom_result_flags(): _mom_cov = mom_cov.copy() _mom_cov[4, i] = np.nan _mom_cov[i, 4] = np.nan - res = make_mom_result(mom, _mom_cov, mom_norm=1) + res = make_mom_result(mom, _mom_cov, sums_norm=1) assert (res["flags"] & ngmix.flags.NONPOS_SHAPE_VAR) != 0 assert ngmix.flags.NAME_MAP[ngmix.flags.NONPOS_SHAPE_VAR] in res["flagstr"] assert res["flux_flags"] == 0 @@ -940,7 +940,7 @@ def test_prepsfmom_mom_norm( ).go( obs=obs, no_psf=True, ) - assert_allclose(res["mom_norm"], res["flux"], atol=0, rtol=2e-4) + assert_allclose(res["sums_norm"], res["flux"], atol=0, rtol=2e-4) @pytest.mark.parametrize('pixel_scale', [0.25, 0.125]) From e731ddf6b3924e76a096722155a175acc56661a2 Mon Sep 17 00:00:00 2001 From: "Matthew R. Becker" Date: Thu, 21 Jul 2022 15:34:05 -0500 Subject: [PATCH 14/14] Update CHANGES.md --- CHANGES.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CHANGES.md b/CHANGES.md index d3f4754f..cc03a1d1 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -3,7 +3,9 @@ ### bug fixes - Fixed 1/area normalization for all moments fields in `GaussMom` results. +### changes + - Changed "mom*" fields in the moments outputs to "sums*". ### new features - Added `fwhm_smooth` keyword to pre-PSF moments routines to allow for extra