diff --git a/CHANGES.md b/CHANGES.md index 2f136baa..cc03a1d1 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -1,10 +1,17 @@ ## v2.1.0 +### 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 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 diff --git a/mdet_tests/test_mdet_regression.py b/mdet_tests/test_mdet_regression.py index e48e565b..d202b523 100644 --- a/mdet_tests/test_mdet_regression.py +++ b/mdet_tests/test_mdet_regression.py @@ -232,7 +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: + 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/gaussmom.py b/ngmix/gaussmom.py index 9ebd663e..75c5d62c 100644 --- a/ngmix/gaussmom.py +++ b/ngmix/gaussmom.py @@ -58,6 +58,11 @@ def _measure_moments(self, obs): res['flux'] *= fac res['flux_err'] *= fac res['pars'][5] *= fac + res["sums"] *= fac + res["sums_cov"] *= fac**2 + res["sums_norm"] *= fac + res["wsum"] *= fac + res["sums_err"] *= fac return res def _set_mompars(self): diff --git a/ngmix/gmix/gmix.py b/ngmix/gmix/gmix.py index 30782710..e0066fc6 100644 --- a/ngmix/gmix/gmix.py +++ b/ngmix/gmix/gmix.py @@ -1239,15 +1239,13 @@ 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] res.update( - moments.make_mom_result(res["sums"], res["sums_cov"]) + moments.make_mom_result(res["sums"].copy(), res["sums_cov"].copy(), res["wsum"]) ) return res diff --git a/ngmix/moments.py b/ngmix/moments.py index c5120d59..81ca936c 100644 --- a/ngmix/moments.py +++ b/ngmix/moments.py @@ -349,15 +349,18 @@ def g2mom(g1, g2, T): return Irr, Irc, Icc -def make_mom_result(mom, mom_cov): - """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. + 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. Returns ------- @@ -365,12 +368,12 @@ def make_mom_result(mom, mom_cov): 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." ) @@ -386,9 +389,10 @@ def make_mom_result(mom, mom_cov): res = {} res["flags"] = 0 res["flagstr"] = "" - res["flux"] = mom[mf_ind] - res["mom"] = mom - res["mom_cov"] = mom_cov + 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 @@ -404,25 +408,25 @@ def make_mom_result(mom, mom_cov): 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 @@ -431,21 +435,21 @@ def make_mom_result(mom, mom_cov): 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"], @@ -454,18 +458,18 @@ def make_mom_result(mom, mom_cov): 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 9ca4361d..135cd7a5 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, sums_norm=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 = kernels["fk00"] 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=knrm, ) @@ -679,6 +681,7 @@ def _gauss_kernels( fkc=fkc, msk=msk, nrm=nrm, + fk00=knrm, ) diff --git a/ngmix/tests/test_prepsfmom.py b/ngmix/tests/test_prepsfmom.py index 99718cfd..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) + 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) + 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) + 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) + 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 @@ -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["sums_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)