Merge pull request #219 from esheldon/wmom-bugz

ENH add moment weight sums to output

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

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:

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):

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

ngmix/moments.py

@@ -349,28 +349,31 @@ 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
     -------
     res : dict
         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

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,
     )