From 3d687dbfdd506dfcb6b8c30a0e9d3c22b1d9d837 Mon Sep 17 00:00:00 2001
From: conda-forge-admin <pelson.pub+conda-forge@gmail.com>
Date: Sun, 19 Jun 2022 07:17:22 -0500
Subject: [PATCH 01/15] WIP extra smoothing in prepsf moments

---
 ngmix/prepsfmom.py | 41 ++++++++++++++++++++++++++++++++++++-----
 1 file changed, 36 insertions(+), 5 deletions(-)

diff --git a/ngmix/prepsfmom.py b/ngmix/prepsfmom.py
index 2994e9bd..d365ffb2 100644
--- a/ngmix/prepsfmom.py
+++ b/ngmix/prepsfmom.py
@@ -37,11 +37,12 @@ class PrePSFMom(object):
         The apodization radius for the stamp in pixels. The default of 1.5 is likely
         fine for most ground based surveys.
     """
-    def __init__(self, fwhm, kernel, pad_factor=4, ap_rad=1.5):
+    def __init__(self, fwhm, kernel, pad_factor=4, ap_rad=1.5, fwhm_smooth=0):
         self.fwhm = fwhm
         self.pad_factor = pad_factor
         self.kernel = kernel
         self.ap_rad = ap_rad
+        self.fwhm_smooth = fwhm_smooth
         if self.kernel == "ksigma":
             self.kind = "ksigma"
         elif self.kernel in ["gauss", "pgauss"]:
@@ -132,6 +133,7 @@ def _meas(self, obs, psf_obs, return_kernels):
                 self.fwhm,
                 obs.jacobian.dvdrow, obs.jacobian.dvdcol,
                 obs.jacobian.dudrow, obs.jacobian.dudcol,
+                fwhm_smooth=self.fwhm_smooth,
             )
         elif self.kernel in ["gauss", "pgauss"]:
             kernels = _gauss_kernels(
@@ -139,6 +141,7 @@ def _meas(self, obs, psf_obs, return_kernels):
                 self.fwhm,
                 obs.jacobian.dvdrow, obs.jacobian.dvdcol,
                 obs.jacobian.dudrow, obs.jacobian.dudcol,
+                fwhm_smooth=self.fwhm_smooth,
             )
         else:
             raise ValueError(
@@ -193,8 +196,11 @@ class KSigmaMom(PrePSFMom):
         The apodization radius for the stamp in pixels. The default of 1.5 is likely
         fine for most ground based surveys.
     """
-    def __init__(self, fwhm, pad_factor=4, ap_rad=1.5):
-        super().__init__(fwhm, 'ksigma', pad_factor=pad_factor, ap_rad=ap_rad)
+    def __init__(self, fwhm, pad_factor=4, ap_rad=1.5, fwhm_smooth=0):
+        super().__init__(
+            fwhm, 'ksigma', pad_factor=pad_factor, ap_rad=ap_rad,
+            fwhm_smooth=fwhm_smooth,
+        )
 
 
 class PGaussMom(PrePSFMom):
@@ -217,8 +223,11 @@ class PGaussMom(PrePSFMom):
         The apodization radius for the stamp in pixels. The default of 1.5 is likely
         fine for most ground based surveys.
     """
-    def __init__(self, fwhm, pad_factor=4, ap_rad=1.5):
-        super().__init__(fwhm, 'pgauss', pad_factor=pad_factor, ap_rad=ap_rad)
+    def __init__(self, fwhm, pad_factor=4, ap_rad=1.5, fwhm_smooth=0):
+        super().__init__(
+            fwhm, 'pgauss', pad_factor=pad_factor, ap_rad=ap_rad,
+            fwhm_smooth=fwhm_smooth,
+        )
 
 
 # keep this here for API consistency
@@ -402,6 +411,7 @@ def _ksigma_kernels(
     dim,
     kernel_size,
     dvdrow, dvdcol, dudrow, dudcol,
+    fwhm_smooth=0,
 ):
     """This function builds a ksigma kernel in Fourier-space.
 
@@ -472,6 +482,16 @@ def _ksigma_kernels(
             "norm = %f (should be 1)!" % (kernel_size, nrm)
         )
 
+    # add smoothing after norm check above for kernel size
+    if fwhm_smooth > 0:
+        sigma_smooth = fwhm_to_sigma(fwhm_smooth)
+        chi2_2_smooth = sigma_smooth * sigma_smooth / 2 * fmag2
+        exp_val_smooth = np.zeros_like(karg)
+        msk_smooth = (chi2_2_smooth < FASTEXP_MAX_CHI2/2) & (chi2_2_smooth >= 0)
+        exp_val_smooth[msk_smooth] = fexp_arr(-chi2_2_smooth[msk_smooth])
+        fkf *= exp_val_smooth
+        knrm *= exp_val_smooth
+
     # the moment kernels take a bit more work
     # product by u^2 in real space is -dk^2/dku^2 in Fourier space
     # same holds for v and cross deriv is -dk^2/dkudkv
@@ -509,6 +529,7 @@ def _gauss_kernels(
     dim,
     kernel_size,
     dvdrow, dvdcol, dudrow, dudcol,
+    fwhm_smooth=0,
 ):
     """This function builds a Gaussian kernel in Fourier-space.
 
@@ -564,6 +585,16 @@ def _gauss_kernels(
             "norm = %f (should be 1)!" % (kernel_size, nrm)
         )
 
+    # add smoothing after norm check above for kernel size
+    if fwhm_smooth > 0:
+        sigma_smooth = fwhm_to_sigma(fwhm_smooth)
+        chi2_2_smooth = sigma_smooth * sigma_smooth / 2 * fmag2
+        exp_val_smooth = np.zeros_like(exp_val)
+        msk_smooth = (chi2_2_smooth < FASTEXP_MAX_CHI2/2) & (chi2_2_smooth >= 0)
+        exp_val_smooth[msk_smooth] = fexp_arr(-chi2_2_smooth[msk_smooth])
+        fkf *= exp_val_smooth
+        knrm *= exp_val_smooth
+
     # the moment kernels take a bit more work
     # product by u^2 in real space is -dk^2/dku^2 in Fourier space
     # same holds for v and cross deriv is -dk^2/dkudkv

From a17168d3a9060b0841ad91f429066cd9958cb190 Mon Sep 17 00:00:00 2001
From: conda-forge-admin <pelson.pub+conda-forge@gmail.com>
Date: Sun, 19 Jun 2022 14:15:12 -0500
Subject: [PATCH 02/15] TST add shear bands for mdet tests

---
 mdet_tests/test_mdet_regression.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/mdet_tests/test_mdet_regression.py b/mdet_tests/test_mdet_regression.py
index c4047a80..e48e565b 100644
--- a/mdet_tests/test_mdet_regression.py
+++ b/mdet_tests/test_mdet_regression.py
@@ -265,7 +265,7 @@ def test_mdet_regression(fname, write=False):
                     ),
                 }
             else:
-                assert col in ["shear"]
+                assert col in ["shear", "shear_bands"]
 
 
 if __name__ == "__main__":

From 3e19aaf6d44bb1b776a43245a03d72986c15f0ae Mon Sep 17 00:00:00 2001
From: beckermr <becker.mr@gmail.com>
Date: Wed, 29 Jun 2022 06:01:22 -0500
Subject: [PATCH 03/15] DOC added docs and inv pixel var weighting

---
 ngmix/prepsfmom.py            |  72 ++++++++++++++---
 ngmix/tests/test_prepsfmom.py | 148 +++++++++++++++++++---------------
 2 files changed, 146 insertions(+), 74 deletions(-)

diff --git a/ngmix/prepsfmom.py b/ngmix/prepsfmom.py
index b65faeb4..a058d55f 100644
--- a/ngmix/prepsfmom.py
+++ b/ngmix/prepsfmom.py
@@ -37,13 +37,26 @@ class PrePSFMom(object):
     ap_rad : float, optional
         The apodization radius for the stamp in pixels. The default of 1.5 is likely
         fine for most ground based surveys.
+    fwhm_smooth : float, optional
+        If non-zero, this optional applies additional Gaussian smoothing to the
+        object before computing the moments. Typically a non-zero value results
+        in less shape noise.
+    use_pix_weight : bool, optional
+        If `True`, this option applies inverse pixel variance weighting in Fourier
+        space when computing the moments. Typically this optional results in poorer
+        performance, but for round PSFs, it results in similar performance as
+        `fwhm_smooth`.
     """
-    def __init__(self, fwhm, kernel, pad_factor=4, ap_rad=1.5, fwhm_smooth=0):
+    def __init__(
+        self, fwhm, kernel, pad_factor=4, ap_rad=1.5,
+        fwhm_smooth=0, use_pix_weight=False
+    ):
         self.fwhm = fwhm
         self.pad_factor = pad_factor
         self.kernel = kernel
         self.ap_rad = ap_rad
         self.fwhm_smooth = fwhm_smooth
+        self.use_pix_weight = use_pix_weight
         if self.kernel == "ksigma":
             self.kind = "ksigma"
         elif self.kernel in ["gauss", "pgauss"]:
@@ -157,6 +170,7 @@ def _meas(self, obs, psf_obs, return_kernels):
         mom, mom_cov = _measure_moments_fft(
             kim, kpsf_im, tot_var, eff_pad_factor, kernels,
             im_row - psf_im_row, im_col - psf_im_col,
+            use_pix_weight=self.use_pix_weight,
         )
         res = make_mom_result(mom, mom_cov)
         if res['flags'] != 0:
@@ -196,11 +210,22 @@ class KSigmaMom(PrePSFMom):
     ap_rad : float, optional
         The apodization radius for the stamp in pixels. The default of 1.5 is likely
         fine for most ground based surveys.
+    fwhm_smooth : float, optional
+        If non-zero, this optional applies additional Gaussian smoothing to the
+        object before computing the moments. Typically a non-zero value results
+        in less shape noise.
+    use_pix_weight : bool, optional
+        If `True`, this option applies inverse pixel variance weighting in Fourier
+        space when computing the moments. Typically this optional results in poorer
+        performance, but for round PSFs, it results in similar performance as
+        `fwhm_smooth`.
     """
-    def __init__(self, fwhm, pad_factor=4, ap_rad=1.5, fwhm_smooth=0):
+    def __init__(
+        self, fwhm, pad_factor=4, ap_rad=1.5, fwhm_smooth=0, use_pix_weight=False,
+    ):
         super().__init__(
             fwhm, 'ksigma', pad_factor=pad_factor, ap_rad=ap_rad,
-            fwhm_smooth=fwhm_smooth,
+            fwhm_smooth=fwhm_smooth, use_pix_weight=use_pix_weight,
         )
 
 
@@ -223,11 +248,22 @@ class PGaussMom(PrePSFMom):
     ap_rad : float, optional
         The apodization radius for the stamp in pixels. The default of 1.5 is likely
         fine for most ground based surveys.
+    fwhm_smooth : float, optional
+        If non-zero, this optional applies additional Gaussian smoothing to the
+        object before computing the moments. Typically a non-zero value results
+        in less shape noise.
+    use_pix_weight : bool, optional
+        If `True`, this option applies inverse pixel variance weighting in Fourier
+        space when computing the moments. Typically this optional results in poorer
+        performance, but for round PSFs, it results in similar performance as
+        `fwhm_smooth`.
     """
-    def __init__(self, fwhm, pad_factor=4, ap_rad=1.5, fwhm_smooth=0):
+    def __init__(
+        self, fwhm, pad_factor=4, ap_rad=1.5, fwhm_smooth=0, use_pix_weight=False,
+    ):
         super().__init__(
             fwhm, 'pgauss', pad_factor=pad_factor, ap_rad=ap_rad,
-            fwhm_smooth=fwhm_smooth,
+            fwhm_smooth=fwhm_smooth, use_pix_weight=use_pix_weight,
         )
 
 
@@ -235,7 +271,9 @@ def __init__(self, fwhm, pad_factor=4, ap_rad=1.5, fwhm_smooth=0):
 PrePSFGaussMom = PGaussMom
 
 
-def _measure_moments_fft(kim, kpsf_im, tot_var, eff_pad_factor, kernels, drow, dcol):
+def _measure_moments_fft(
+    kim, kpsf_im, tot_var, eff_pad_factor, kernels, drow, dcol, use_pix_weight=False,
+):
     # we only need to do things where the kernel is non-zero
     # this saves a bunch of CPU cycles
     msk = kernels["msk"]
@@ -271,10 +309,15 @@ def _measure_moments_fft(kim, kpsf_im, tot_var, eff_pad_factor, kernels, drow, d
     fkp = kernels["fkp"]
     fkc = kernels["fkc"]
 
-    mf = np.sum((kim * fkf).real) * df2
-    mr = np.sum((kim * fkr).real) * df2
-    mp = np.sum((kim * fkp).real) * df2
-    mc = np.sum((kim * fkc).real) * df2
+    if use_pix_weight:
+        wgt = (kpsf_im * np.conj(kpsf_im)).real
+    else:
+        wgt = 1.0
+
+    mf = np.sum((kim * fkf * wgt).real) * df2
+    mr = np.sum((kim * fkr * wgt).real) * df2
+    mp = np.sum((kim * fkp * wgt).real) * df2
+    mc = np.sum((kim * fkc * wgt).real) * df2
 
     # build a covariance matrix of the moments
     # here we assume each Fourier mode is independent and sum the variances
@@ -288,10 +331,17 @@ def _measure_moments_fft(kim, kpsf_im, tot_var, eff_pad_factor, kernels, drow, d
     m_cov[1, 1] = 1
     tot_var *= eff_pad_factor**2
     tot_var_df4 = tot_var * df4
-    kerns = [fkp / kpsf_im, fkc / kpsf_im, fkr / kpsf_im, fkf / kpsf_im]
+    psf_kerns_fac = wgt / kpsf_im
+    kerns = [
+        fkp * psf_kerns_fac,
+        fkc * psf_kerns_fac,
+        fkr * psf_kerns_fac,
+        fkf * psf_kerns_fac,
+    ]
     conj_kerns = [np.conj(k) for k in kerns]
     for i in range(2, 6):
         for j in range(i, 6):
+            # subtract two since kernels start at second moments
             m_cov[i, j] = np.sum((kerns[i-2] * conj_kerns[j-2]).real) * tot_var_df4
             m_cov[j, i] = m_cov[i, j]
 
diff --git a/ngmix/tests/test_prepsfmom.py b/ngmix/tests/test_prepsfmom.py
index f82e6aec..c5bb5f0b 100644
--- a/ngmix/tests/test_prepsfmom.py
+++ b/ngmix/tests/test_prepsfmom.py
@@ -3,6 +3,7 @@
 import pytest
 import time
 from flaky import flaky
+from numpy.testing import assert_allclose
 
 from ngmix.prepsfmom import (
     KSigmaMom, PGaussMom,
@@ -23,6 +24,7 @@ def _report_info(s, arr, mn, err):
             "%s:" % s,
             np.mean(arr), mn, np.mean(arr)/mn - 1,
             np.std(arr), err, np.std(arr)/err - 1,
+            np.abs(np.mean(arr))/np.std(arr),
             flush=True,
         )
     else:
@@ -30,6 +32,7 @@ def _report_info(s, arr, mn, err):
             "%s:" % s,
             np.mean(arr), None, None,
             np.std(arr), None, None,
+            None,
             flush=True,
         )
 
@@ -370,34 +373,34 @@ def test_prepsfmom_gauss(
         _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)
-        assert np.allclose(
+        assert_allclose(
             np.abs(np.mean(res["flux"]) - flux_true)/np.mean(res["flux_err"]),
             0,
             atol=4,
             rtol=0,
         )
-        assert np.allclose(
+        assert_allclose(
             np.mean(res["flux"]), flux_true, atol=0, rtol=0.1)
-        assert np.allclose(
+        assert_allclose(
             np.std(res["flux"]), np.mean(res["flux_err"]), atol=0, rtol=0.2)
 
 
 @pytest.mark.parametrize("cls,mom_fwhm,snr", [
-    (KSigmaMom, 2.0, 1e2),
+    # (KSigmaMom, 2.0, 1e2),
     (PGaussMom, 2.0, 1e2),
 ])
 @pytest.mark.parametrize('pixel_scale', [0.25])
-@pytest.mark.parametrize('fwhm,psf_fwhm', [
-    (2.0, 1.0),
-])
-@pytest.mark.parametrize('image_size', [
-    53,
-])
-@pytest.mark.parametrize('pad_factor', [
-    1.5,
-])
-def test_prepsfmom_mn_cov(
+@pytest.mark.parametrize('fwhm,psf_fwhm', [(2.0, 1.0)])
+@pytest.mark.parametrize('image_size', [53])
+@pytest.mark.parametrize('pad_factor', [1.5])
+@pytest.mark.parametrize('round', [True, False])
+@pytest.mark.parametrize(
+    'fwhm_smooth,use_pix_weight',
+    [(0, False), (0, True), (1, False)],
+)
+def test_prepsfmom_mn_cov_psf(
     pad_factor, image_size, fwhm, psf_fwhm, pixel_scale, snr, mom_fwhm, cls,
+    use_pix_weight, fwhm_smooth, round,
 ):
     """Slower test to make sure means and errors are right
     w/ tons of monte carlo samples.
@@ -434,9 +437,11 @@ def test_prepsfmom_mn_cov(
     )
     psf = galsim.Gaussian(
         fwhm=psf_fwhm
-    ).shear(
-        g1=0.3, g2=-0.15
     )
+    if not round:
+        psf = psf.shear(
+            g1=0.3, g2=-0.15
+        )
     im = galsim.Convolve([gal, psf]).drawImage(
         nx=image_size,
         ny=image_size,
@@ -456,19 +461,35 @@ def test_prepsfmom_mn_cov(
     fitter = cls(
         fwhm=mom_fwhm,
         pad_factor=pad_factor,
+        use_pix_weight=use_pix_weight,
+        fwhm_smooth=fwhm_smooth,
     )
 
     # get true flux
-    im_true = gal.drawImage(
-        nx=image_size,
-        ny=image_size,
-        wcs=gs_wcs,
-        method='no_pixel').array
-    obs = Observation(
-        image=im_true,
-        jacobian=jac,
-    )
-    res = cls(fwhm=mom_fwhm, pad_factor=pad_factor).go(obs=obs, no_psf=True)
+    if use_pix_weight:
+        im_true = galsim.Convolve([gal, psf]).drawImage(
+            nx=image_size,
+            ny=image_size,
+            wcs=gs_wcs,
+            method='no_pixel').array
+        obs = Observation(
+            image=im_true,
+            weight=wgt,
+            jacobian=jac,
+            psf=Observation(image=psf_im, jacobian=psf_jac),
+        )
+        res = fitter.go(obs=obs)
+    else:
+        im_true = gal.drawImage(
+            nx=image_size,
+            ny=image_size,
+            wcs=gs_wcs,
+            method='no_pixel').array
+        obs = Observation(
+            image=im_true,
+            jacobian=jac,
+        )
+        res = fitter.go(obs=obs, no_psf=True)
     flux_true = res["flux"]
     T_true = res["T"]
     g1_true = res["e"][0]
@@ -501,26 +522,26 @@ def test_prepsfmom_mn_cov(
     print("mom cov meas:\n", mom_cov, flush=True)
     print("mom cov pred:\n", np.mean(res["mom_cov"], axis=0), flush=True)
 
-    assert np.allclose(np.mean(res["flux"]), flux_true, atol=0, rtol=1e-2)
-    assert np.allclose(np.mean(res["T"]), T_true, atol=0, rtol=1e-2)
-    assert np.allclose(np.mean(res["e"][:, 0]), g1_true, atol=0, rtol=1e-2)
-    assert np.allclose(np.mean(res["e"][:, 1]), g2_true, atol=0, rtol=1e-2)
+    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)
+    assert_allclose(np.mean(res["e"][:, 0]), g1_true, atol=0, rtol=2e-2)
+    assert_allclose(np.mean(res["e"][:, 1]), g2_true, atol=0, rtol=2e-2)
 
-    assert np.allclose(np.std(res["flux"]), np.mean(res["flux_err"]), atol=0, rtol=2e-2)
-    assert np.allclose(np.std(res["T"]), np.mean(res["T_err"]), atol=0, rtol=2e-2)
-    assert np.allclose(
+    assert_allclose(np.std(res["flux"]), np.mean(res["flux_err"]), atol=0, rtol=2e-2)
+    assert_allclose(np.std(res["T"]), np.mean(res["T_err"]), atol=0, rtol=2e-2)
+    assert_allclose(
         np.std(res["e"][:, 0]), np.mean(res["e_err"][:, 0]), atol=0, rtol=2e-2)
-    assert np.allclose(
+    assert_allclose(
         np.std(res["e"][:, 1]), np.mean(res["e_err"][:, 1]), atol=0, rtol=2e-2)
 
-    assert np.allclose(
+    assert_allclose(
         mom_cov[2:, 2:],
         np.mean(res["mom_cov"][:, 2:, 2:], axis=0),
         atol=2.5e-1,
         rtol=0,
     )
 
-    assert np.allclose(
+    assert_allclose(
         np.diagonal(mom_cov[2:, 2:]),
         np.diagonal(np.mean(res["mom_cov"][:, 2:, 2:], axis=0)),
         atol=0,
@@ -529,21 +550,16 @@ def test_prepsfmom_mn_cov(
 
 
 @pytest.mark.parametrize("cls,mom_fwhm,snr", [
-    (KSigmaMom, 2.0, 1e2),
     (PGaussMom, 2.0, 1e2),
+    (KSigmaMom, 2.0, 1e2),
 ])
 @pytest.mark.parametrize('pixel_scale', [0.25])
-@pytest.mark.parametrize('fwhm', [
-    2,
-])
-@pytest.mark.parametrize('image_size', [
-    53,
-])
-@pytest.mark.parametrize('pad_factor', [
-    1.5,
-])
+@pytest.mark.parametrize('fwhm', [2])
+@pytest.mark.parametrize('image_size', [53])
+@pytest.mark.parametrize('pad_factor', [1.5])
+@pytest.mark.parametrize('use_pix_weight', [True, False])
 def test_prepsfmom_mn_cov_nopsf(
-    pad_factor, image_size, fwhm, pixel_scale, snr, mom_fwhm, cls,
+    pad_factor, image_size, fwhm, pixel_scale, snr, mom_fwhm, cls, use_pix_weight,
 ):
     """Slower test to make sure means and errors are right
     w/ tons of monte carlo samples.
@@ -582,6 +598,7 @@ def test_prepsfmom_mn_cov_nopsf(
     fitter = cls(
         fwhm=mom_fwhm,
         pad_factor=pad_factor,
+        use_pix_weight=True,
     )
 
     # get true flux
@@ -594,7 +611,9 @@ def test_prepsfmom_mn_cov_nopsf(
         image=im_true,
         jacobian=jac,
     )
-    res = cls(fwhm=mom_fwhm, pad_factor=pad_factor).go(obs=obs, no_psf=True)
+    res = cls(
+        fwhm=mom_fwhm, pad_factor=pad_factor, use_pix_weight=True,
+    ).go(obs=obs, no_psf=True)
     flux_true = res["flux"]
     T_true = res["T"]
     g1_true = res["e"][0]
@@ -626,26 +645,26 @@ def test_prepsfmom_mn_cov_nopsf(
     print("mom cov meas:\n", mom_cov, flush=True)
     print("mom cov pred:\n", np.mean(res["mom_cov"], axis=0), flush=True)
 
-    assert np.allclose(np.mean(res["flux"]), flux_true, atol=0, rtol=1e-2)
-    assert np.allclose(np.mean(res["T"]), T_true, atol=0, rtol=1e-2)
-    assert np.allclose(np.mean(res["e"][:, 0]), g1_true, atol=0, rtol=1e-2)
-    assert np.allclose(np.mean(res["e"][:, 1]), g2_true, atol=0, rtol=1e-2)
+    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)
+    assert_allclose(np.mean(res["e"][:, 0]), g1_true, atol=0, rtol=1e-2)
+    assert_allclose(np.mean(res["e"][:, 1]), g2_true, atol=0, rtol=1e-2)
 
-    assert np.allclose(np.std(res["flux"]), np.mean(res["flux_err"]), atol=0, rtol=2e-2)
-    assert np.allclose(np.std(res["T"]), np.mean(res["T_err"]), atol=0, rtol=2e-2)
-    assert np.allclose(
+    assert_allclose(np.std(res["flux"]), np.mean(res["flux_err"]), atol=0, rtol=2e-2)
+    assert_allclose(np.std(res["T"]), np.mean(res["T_err"]), atol=0, rtol=2e-2)
+    assert_allclose(
         np.std(res["e"][:, 0]), np.mean(res["e_err"][:, 0]), atol=0, rtol=2e-2)
-    assert np.allclose(
+    assert_allclose(
         np.std(res["e"][:, 1]), np.mean(res["e_err"][:, 1]), atol=0, rtol=2e-2)
 
-    assert np.allclose(
+    assert_allclose(
         mom_cov[2:, 2:],
         np.mean(res["mom_cov"][:, 2:, 2:], axis=0),
         atol=2.5e-1,
         rtol=0,
     )
 
-    assert np.allclose(
+    assert_allclose(
         np.diagonal(mom_cov[2:, 2:]),
         np.diagonal(np.mean(res["mom_cov"][:, 2:, 2:], axis=0)),
         atol=0,
@@ -795,7 +814,7 @@ def test_prepsfmom_gauss_true_flux(
     )
     res = fitter.go(obs=obs, no_psf=True)
     flux_true = res["flux"]
-    assert np.allclose(flux_true, 400, atol=0, rtol=5e-3)
+    assert_allclose(flux_true, 400, atol=0, rtol=5e-3)
 
     obs = Observation(
         image=im,
@@ -805,7 +824,7 @@ def test_prepsfmom_gauss_true_flux(
     )
     res = fitter.go(obs=obs)
     flux_true = res["flux"]
-    assert np.allclose(flux_true, 400, atol=0, rtol=5e-3)
+    assert_allclose(flux_true, 400, atol=0, rtol=5e-3)
 
 
 @pytest.mark.parametrize('pixel_scale', [0.25, 0.125])
@@ -821,8 +840,9 @@ def test_prepsfmom_gauss_true_flux(
 @pytest.mark.parametrize('mom_fwhm', [
     2, 2.5,
 ])
+@pytest.mark.parametrize('use_pix_weight', [True, False])
 def test_prepsfmom_comp_to_gaussmom(
-    pad_factor, image_size, fwhm, pixel_scale, mom_fwhm,
+    pad_factor, image_size, fwhm, pixel_scale, mom_fwhm, use_pix_weight,
 ):
     rng = np.random.RandomState(seed=100)
 
@@ -858,7 +878,9 @@ def test_prepsfmom_comp_to_gaussmom(
         image=im_true,
         jacobian=jac,
     )
-    res = PGaussMom(fwhm=mom_fwhm, pad_factor=pad_factor).go(
+    res = PGaussMom(
+        fwhm=mom_fwhm, pad_factor=pad_factor, use_pix_weight=use_pix_weight
+    ).go(
         obs=obs, no_psf=True, return_kernels=True,
     )
 
@@ -870,7 +892,7 @@ def test_prepsfmom_comp_to_gaussmom(
             print("%s:" % k, res[k], res_gmom[k])
 
     for k in ["flux", "flux_err", "T", "T_err", "e", "e_cov"]:
-        assert np.allclose(res[k], res_gmom[k], atol=0, rtol=1e-2)
+        assert_allclose(res[k], res_gmom[k], atol=0, rtol=1e-2)
 
 
 def _sim_apodize(flux_factor, ap_rad):

From 3e7cfcab5c60218b42f8902b26e40852c46154d5 Mon Sep 17 00:00:00 2001
From: beckermr <becker.mr@gmail.com>
Date: Wed, 29 Jun 2022 06:48:33 -0500
Subject: [PATCH 04/15] PROD bump version

---
 ngmix/_version.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/ngmix/_version.py b/ngmix/_version.py
index df7fede3..d1ba3084 100644
--- a/ngmix/_version.py
+++ b/ngmix/_version.py
@@ -1 +1 @@
-__version__ = '2.0.6'  # noqa
+__version__ = '2.1.0'  # noqa

From a3f06f80c9f35cada26535c85d68da56ed98b286 Mon Sep 17 00:00:00 2001
From: beckermr <becker.mr@gmail.com>
Date: Wed, 29 Jun 2022 06:50:10 -0500
Subject: [PATCH 05/15] DOC update the changelog

---
 CHANGES.md | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/CHANGES.md b/CHANGES.md
index 896b5c18..98811992 100644
--- a/CHANGES.md
+++ b/CHANGES.md
@@ -1,3 +1,13 @@
+## v2.1.0
+
+### 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 `use_pix_weight` keyword to pre-PSF moments routines to enable inverse
+      variance pixel weighting.
+
+
 ## v2.0.6
 
 ### bug fixes

From 71866c74ea2a7bccbb5364e45e490638926ee039 Mon Sep 17 00:00:00 2001
From: beckermr <becker.mr@gmail.com>
Date: Wed, 29 Jun 2022 07:04:00 -0500
Subject: [PATCH 06/15] DOC & TST added change log, try more

---
 CHANGES.md                    | 1 +
 ngmix/tests/test_prepsfmom.py | 2 +-
 2 files changed, 2 insertions(+), 1 deletion(-)

diff --git a/CHANGES.md b/CHANGES.md
index 98811992..749ec314 100644
--- a/CHANGES.md
+++ b/CHANGES.md
@@ -6,6 +6,7 @@
       smoothing of the profile before the moments are measured.
     - Added `use_pix_weight` keyword to pre-PSF moments routines to enable inverse
       variance pixel weighting.
+    - Added caching of FFTs in metacal and pre-PSF moment rountines.
 
 
 ## v2.0.6
diff --git a/ngmix/tests/test_prepsfmom.py b/ngmix/tests/test_prepsfmom.py
index c5bb5f0b..cb7724cc 100644
--- a/ngmix/tests/test_prepsfmom.py
+++ b/ngmix/tests/test_prepsfmom.py
@@ -101,7 +101,7 @@ def test_prepsfmom_raises_badjacob(cls):
     assert "same WCS Jacobia" in str(e.value)
 
 
-@flaky
+@flaky(max_runs=10)
 def test_prepsfmom_speed_and_cache():
     image_size = 48
     psf_image_size = 53

From bd14129c2c6cd4348e6e7ca1bc544db4e59e5283 Mon Sep 17 00:00:00 2001
From: beckermr <becker.mr@gmail.com>
Date: Wed, 29 Jun 2022 07:15:41 -0500
Subject: [PATCH 07/15] TST retry mcal cache too

---
 ngmix/tests/test_metacal_cache.py | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/ngmix/tests/test_metacal_cache.py b/ngmix/tests/test_metacal_cache.py
index 43ab99df..4d35b791 100644
--- a/ngmix/tests/test_metacal_cache.py
+++ b/ngmix/tests/test_metacal_cache.py
@@ -1,11 +1,15 @@
 import time
 import numpy as np
+
+from flaky import flaky
+
 import ngmix
 import ngmix.metacal.metacal
 from ._galsim_sims import _get_obs
 from ..metacal.metacal import _cached_galsim_stuff
 
 
+@flaky(max_runs=10)
 def test_metacal_cache():
     # first warm up numba
     rng = np.random.RandomState(seed=100)

From a8e3982c69344ee335f7e622f4f1581df6873d04 Mon Sep 17 00:00:00 2001
From: beckermr <becker.mr@gmail.com>
Date: Wed, 29 Jun 2022 08:32:17 -0500
Subject: [PATCH 08/15] TST more tests

---
 ngmix/tests/test_prepsfmom.py | 107 ++++++++++++++++++++++++++++------
 1 file changed, 88 insertions(+), 19 deletions(-)

diff --git a/ngmix/tests/test_prepsfmom.py b/ngmix/tests/test_prepsfmom.py
index cb7724cc..dae1b0ba 100644
--- a/ngmix/tests/test_prepsfmom.py
+++ b/ngmix/tests/test_prepsfmom.py
@@ -557,9 +557,9 @@ def test_prepsfmom_mn_cov_psf(
 @pytest.mark.parametrize('fwhm', [2])
 @pytest.mark.parametrize('image_size', [53])
 @pytest.mark.parametrize('pad_factor', [1.5])
-@pytest.mark.parametrize('use_pix_weight', [True, False])
+@pytest.mark.parametrize('fwhm_smooth', [0, 1])
 def test_prepsfmom_mn_cov_nopsf(
-    pad_factor, image_size, fwhm, pixel_scale, snr, mom_fwhm, cls, use_pix_weight,
+    pad_factor, image_size, fwhm, pixel_scale, snr, mom_fwhm, cls, fwhm_smooth,
 ):
     """Slower test to make sure means and errors are right
     w/ tons of monte carlo samples.
@@ -598,7 +598,7 @@ def test_prepsfmom_mn_cov_nopsf(
     fitter = cls(
         fwhm=mom_fwhm,
         pad_factor=pad_factor,
-        use_pix_weight=True,
+        fwhm_smooth=fwhm_smooth,
     )
 
     # get true flux
@@ -612,7 +612,8 @@ def test_prepsfmom_mn_cov_nopsf(
         jacobian=jac,
     )
     res = cls(
-        fwhm=mom_fwhm, pad_factor=pad_factor, use_pix_weight=True,
+        fwhm=mom_fwhm, pad_factor=pad_factor,
+        fwhm_smooth=fwhm_smooth,
     ).go(obs=obs, no_psf=True)
     flux_true = res["flux"]
     T_true = res["T"]
@@ -828,21 +829,12 @@ def test_prepsfmom_gauss_true_flux(
 
 
 @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,
-])
-@pytest.mark.parametrize('use_pix_weight', [True, False])
+@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(
-    pad_factor, image_size, fwhm, pixel_scale, mom_fwhm, use_pix_weight,
+    pad_factor, image_size, fwhm, pixel_scale, mom_fwhm
 ):
     rng = np.random.RandomState(seed=100)
 
@@ -879,7 +871,7 @@ def test_prepsfmom_comp_to_gaussmom(
         jacobian=jac,
     )
     res = PGaussMom(
-        fwhm=mom_fwhm, pad_factor=pad_factor, use_pix_weight=use_pix_weight
+        fwhm=mom_fwhm, pad_factor=pad_factor,
     ).go(
         obs=obs, no_psf=True, return_kernels=True,
     )
@@ -895,6 +887,83 @@ def test_prepsfmom_comp_to_gaussmom(
         assert_allclose(res[k], res_gmom[k], atol=0, rtol=1e-2)
 
 
+@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])
+@pytest.mark.parametrize('fwhm_smooth', [0, 1.5])
+def test_prepsfmom_comp_to_gaussmom_fwhm_smooth(
+    pad_factor, image_size, fwhm, pixel_scale, mom_fwhm, fwhm_smooth
+):
+    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)
+
+    gal = galsim.Gaussian(
+        fwhm=fwhm
+    ).shear(
+        g1=-0.1, g2=0.2
+    ).withFlux(
+        400
+    ).shift(
+        dx=shift[0], dy=shift[1]
+    )
+
+    # get true flux
+    im_true = gal.drawImage(
+        nx=image_size,
+        ny=image_size,
+        wcs=gs_wcs,
+    ).array
+    obs = Observation(
+        image=im_true,
+        jacobian=jac,
+    )
+    res = PGaussMom(
+        fwhm=mom_fwhm, pad_factor=pad_factor, fwhm_smooth=fwhm_smooth,
+    ).go(
+        obs=obs, no_psf=True, return_kernels=True,
+    )
+
+    from ngmix.gaussmom import GaussMom
+    if fwhm_smooth > 0:
+        im_true_smooth = galsim.Convolve(
+            [gal, galsim.Gaussian(fwhm=fwhm_smooth)]
+        ).drawImage(
+            nx=image_size,
+            ny=image_size,
+            wcs=gs_wcs,
+        ).array
+    else:
+        im_true_smooth = im_true
+    obs_smooth = Observation(
+        image=im_true_smooth,
+        jacobian=jac,
+    )
+    res_gmom = GaussMom(fwhm=mom_fwhm).go(obs=obs_smooth)
+
+    for k in sorted(res):
+        if k in res_gmom:
+            print("%s:" % k, res[k], res_gmom[k])
+
+    for k in ["flux", "T", "e"]:
+        assert_allclose(res[k], res_gmom[k], atol=0, rtol=1e-2)
+    # the error do not match - not sure why
+    # for k in ["flux_err", "T_err", "e_cov"]:
+    #     assert_allclose(res[k], res_gmom[k], atol=0, rtol=5e-1)
+
+
 def _sim_apodize(flux_factor, ap_rad):
     """
     we are simulating an object at the center with a bright object right on the

From d54011cba13e58e4f25daee1d5432eb5ac170858 Mon Sep 17 00:00:00 2001
From: beckermr <becker.mr@gmail.com>
Date: Wed, 29 Jun 2022 08:41:42 -0500
Subject: [PATCH 09/15] DOC add comment on tyest

---
 ngmix/tests/test_prepsfmom.py | 15 +++++++++------
 1 file changed, 9 insertions(+), 6 deletions(-)

diff --git a/ngmix/tests/test_prepsfmom.py b/ngmix/tests/test_prepsfmom.py
index dae1b0ba..435d370f 100644
--- a/ngmix/tests/test_prepsfmom.py
+++ b/ngmix/tests/test_prepsfmom.py
@@ -933,7 +933,7 @@ def test_prepsfmom_comp_to_gaussmom_fwhm_smooth(
     res = PGaussMom(
         fwhm=mom_fwhm, pad_factor=pad_factor, fwhm_smooth=fwhm_smooth,
     ).go(
-        obs=obs, no_psf=True, return_kernels=True,
+        obs=obs, no_psf=True,
     )
 
     from ngmix.gaussmom import GaussMom
@@ -957,11 +957,14 @@ def test_prepsfmom_comp_to_gaussmom_fwhm_smooth(
         if k in res_gmom:
             print("%s:" % k, res[k], res_gmom[k])
 
-    for k in ["flux", "T", "e"]:
-        assert_allclose(res[k], res_gmom[k], atol=0, rtol=1e-2)
-    # the error do not match - not sure why
-    # for k in ["flux_err", "T_err", "e_cov"]:
-    #     assert_allclose(res[k], res_gmom[k], atol=0, rtol=5e-1)
+    assert_allclose(res["flux"], res_gmom["flux"], atol=0, rtol=5e-4)
+    assert_allclose(res["T"], res_gmom["T"], atol=0, rtol=1e-3)
+    assert_allclose(res["e"], res_gmom["e"], atol=0, rtol=1e-3)
+    # the errors do not match - this is because the underlying noise model is
+    # different - the pure gaussian moments weight map is an error on the convolved
+    # profile whereas the pre-PSF case uses error propagation through the
+    # smoothing kernel treating the weight map as applying to the unconvolved profile
+    # thus we do not test the errors
 
 
 def _sim_apodize(flux_factor, ap_rad):

From 5a25fe886fef2057255123192565c39ac9d47824 Mon Sep 17 00:00:00 2001
From: beckermr <becker.mr@gmail.com>
Date: Thu, 30 Jun 2022 06:11:32 -0500
Subject: [PATCH 10/15] REF remove pixel weighting for now

---
 CHANGES.md                    |  2 --
 ngmix/prepsfmom.py            | 45 +++++++++--------------------------
 ngmix/tests/test_prepsfmom.py | 42 ++++++++++----------------------
 3 files changed, 23 insertions(+), 66 deletions(-)

diff --git a/CHANGES.md b/CHANGES.md
index 749ec314..2f136baa 100644
--- a/CHANGES.md
+++ b/CHANGES.md
@@ -4,8 +4,6 @@
 
     - Added `fwhm_smooth` keyword to pre-PSF moments routines to allow for extra
       smoothing of the profile before the moments are measured.
-    - Added `use_pix_weight` keyword to pre-PSF moments routines to enable inverse
-      variance pixel weighting.
     - Added caching of FFTs in metacal and pre-PSF moment rountines.
 
 
diff --git a/ngmix/prepsfmom.py b/ngmix/prepsfmom.py
index a058d55f..547d8b19 100644
--- a/ngmix/prepsfmom.py
+++ b/ngmix/prepsfmom.py
@@ -41,22 +41,15 @@ class PrePSFMom(object):
         If non-zero, this optional applies additional Gaussian smoothing to the
         object before computing the moments. Typically a non-zero value results
         in less shape noise.
-    use_pix_weight : bool, optional
-        If `True`, this option applies inverse pixel variance weighting in Fourier
-        space when computing the moments. Typically this optional results in poorer
-        performance, but for round PSFs, it results in similar performance as
-        `fwhm_smooth`.
     """
     def __init__(
-        self, fwhm, kernel, pad_factor=4, ap_rad=1.5,
-        fwhm_smooth=0, use_pix_weight=False
+        self, fwhm, kernel, pad_factor=4, ap_rad=1.5, fwhm_smooth=0
     ):
         self.fwhm = fwhm
         self.pad_factor = pad_factor
         self.kernel = kernel
         self.ap_rad = ap_rad
         self.fwhm_smooth = fwhm_smooth
-        self.use_pix_weight = use_pix_weight
         if self.kernel == "ksigma":
             self.kind = "ksigma"
         elif self.kernel in ["gauss", "pgauss"]:
@@ -170,7 +163,6 @@ def _meas(self, obs, psf_obs, return_kernels):
         mom, mom_cov = _measure_moments_fft(
             kim, kpsf_im, tot_var, eff_pad_factor, kernels,
             im_row - psf_im_row, im_col - psf_im_col,
-            use_pix_weight=self.use_pix_weight,
         )
         res = make_mom_result(mom, mom_cov)
         if res['flags'] != 0:
@@ -214,18 +206,13 @@ class KSigmaMom(PrePSFMom):
         If non-zero, this optional applies additional Gaussian smoothing to the
         object before computing the moments. Typically a non-zero value results
         in less shape noise.
-    use_pix_weight : bool, optional
-        If `True`, this option applies inverse pixel variance weighting in Fourier
-        space when computing the moments. Typically this optional results in poorer
-        performance, but for round PSFs, it results in similar performance as
-        `fwhm_smooth`.
     """
     def __init__(
-        self, fwhm, pad_factor=4, ap_rad=1.5, fwhm_smooth=0, use_pix_weight=False,
+        self, fwhm, pad_factor=4, ap_rad=1.5, fwhm_smooth=0
     ):
         super().__init__(
             fwhm, 'ksigma', pad_factor=pad_factor, ap_rad=ap_rad,
-            fwhm_smooth=fwhm_smooth, use_pix_weight=use_pix_weight,
+            fwhm_smooth=fwhm_smooth,
         )
 
 
@@ -252,18 +239,13 @@ class PGaussMom(PrePSFMom):
         If non-zero, this optional applies additional Gaussian smoothing to the
         object before computing the moments. Typically a non-zero value results
         in less shape noise.
-    use_pix_weight : bool, optional
-        If `True`, this option applies inverse pixel variance weighting in Fourier
-        space when computing the moments. Typically this optional results in poorer
-        performance, but for round PSFs, it results in similar performance as
-        `fwhm_smooth`.
     """
     def __init__(
-        self, fwhm, pad_factor=4, ap_rad=1.5, fwhm_smooth=0, use_pix_weight=False,
+        self, fwhm, pad_factor=4, ap_rad=1.5, fwhm_smooth=0,
     ):
         super().__init__(
             fwhm, 'pgauss', pad_factor=pad_factor, ap_rad=ap_rad,
-            fwhm_smooth=fwhm_smooth, use_pix_weight=use_pix_weight,
+            fwhm_smooth=fwhm_smooth,
         )
 
 
@@ -272,7 +254,7 @@ def __init__(
 
 
 def _measure_moments_fft(
-    kim, kpsf_im, tot_var, eff_pad_factor, kernels, drow, dcol, use_pix_weight=False,
+    kim, kpsf_im, tot_var, eff_pad_factor, kernels, drow, dcol,
 ):
     # we only need to do things where the kernel is non-zero
     # this saves a bunch of CPU cycles
@@ -309,15 +291,10 @@ def _measure_moments_fft(
     fkp = kernels["fkp"]
     fkc = kernels["fkc"]
 
-    if use_pix_weight:
-        wgt = (kpsf_im * np.conj(kpsf_im)).real
-    else:
-        wgt = 1.0
-
-    mf = np.sum((kim * fkf * wgt).real) * df2
-    mr = np.sum((kim * fkr * wgt).real) * df2
-    mp = np.sum((kim * fkp * wgt).real) * df2
-    mc = np.sum((kim * fkc * wgt).real) * df2
+    mf = np.sum((kim * fkf).real) * df2
+    mr = np.sum((kim * fkr).real) * df2
+    mp = np.sum((kim * fkp).real) * df2
+    mc = np.sum((kim * fkc).real) * df2
 
     # build a covariance matrix of the moments
     # here we assume each Fourier mode is independent and sum the variances
@@ -331,7 +308,7 @@ def _measure_moments_fft(
     m_cov[1, 1] = 1
     tot_var *= eff_pad_factor**2
     tot_var_df4 = tot_var * df4
-    psf_kerns_fac = wgt / kpsf_im
+    psf_kerns_fac = 1 / kpsf_im
     kerns = [
         fkp * psf_kerns_fac,
         fkc * psf_kerns_fac,
diff --git a/ngmix/tests/test_prepsfmom.py b/ngmix/tests/test_prepsfmom.py
index 435d370f..84b1f2aa 100644
--- a/ngmix/tests/test_prepsfmom.py
+++ b/ngmix/tests/test_prepsfmom.py
@@ -394,13 +394,10 @@ def test_prepsfmom_gauss(
 @pytest.mark.parametrize('image_size', [53])
 @pytest.mark.parametrize('pad_factor', [1.5])
 @pytest.mark.parametrize('round', [True, False])
-@pytest.mark.parametrize(
-    'fwhm_smooth,use_pix_weight',
-    [(0, False), (0, True), (1, False)],
-)
+@pytest.mark.parametrize('fwhm_smooth', [0, 1])
 def test_prepsfmom_mn_cov_psf(
     pad_factor, image_size, fwhm, psf_fwhm, pixel_scale, snr, mom_fwhm, cls,
-    use_pix_weight, fwhm_smooth, round,
+    fwhm_smooth, round,
 ):
     """Slower test to make sure means and errors are right
     w/ tons of monte carlo samples.
@@ -461,35 +458,20 @@ def test_prepsfmom_mn_cov_psf(
     fitter = cls(
         fwhm=mom_fwhm,
         pad_factor=pad_factor,
-        use_pix_weight=use_pix_weight,
         fwhm_smooth=fwhm_smooth,
     )
 
     # get true flux
-    if use_pix_weight:
-        im_true = galsim.Convolve([gal, psf]).drawImage(
-            nx=image_size,
-            ny=image_size,
-            wcs=gs_wcs,
-            method='no_pixel').array
-        obs = Observation(
-            image=im_true,
-            weight=wgt,
-            jacobian=jac,
-            psf=Observation(image=psf_im, jacobian=psf_jac),
-        )
-        res = fitter.go(obs=obs)
-    else:
-        im_true = gal.drawImage(
-            nx=image_size,
-            ny=image_size,
-            wcs=gs_wcs,
-            method='no_pixel').array
-        obs = Observation(
-            image=im_true,
-            jacobian=jac,
-        )
-        res = fitter.go(obs=obs, no_psf=True)
+    im_true = gal.drawImage(
+        nx=image_size,
+        ny=image_size,
+        wcs=gs_wcs,
+        method='no_pixel').array
+    obs = Observation(
+        image=im_true,
+        jacobian=jac,
+    )
+    res = fitter.go(obs=obs, no_psf=True)
     flux_true = res["flux"]
     T_true = res["T"]
     g1_true = res["e"][0]

From 205e1ca55cb8d8349c98db3572b2682c82da2a3d Mon Sep 17 00:00:00 2001
From: beckermr <becker.mr@gmail.com>
Date: Thu, 30 Jun 2022 06:14:02 -0500
Subject: [PATCH 11/15] REF remove duplicated code

---
 ngmix/prepsfmom.py | 21 +++++++++++----------
 1 file changed, 11 insertions(+), 10 deletions(-)

diff --git a/ngmix/prepsfmom.py b/ngmix/prepsfmom.py
index 547d8b19..9ca4361d 100644
--- a/ngmix/prepsfmom.py
+++ b/ngmix/prepsfmom.py
@@ -459,6 +459,15 @@ def _deconvolve_im_psf_inplace(kim, kpsf_im, max_amp, min_psf_frac=1e-5):
     return kim, kpsf_im, msk
 
 
+def _get_fwhm_smooth_profile(fwhm_smooth, fmag2):
+    sigma_smooth = fwhm_to_sigma(fwhm_smooth)
+    chi2_2_smooth = sigma_smooth * sigma_smooth / 2 * fmag2
+    exp_val_smooth = np.zeros_like(fmag2)
+    msk_smooth = (chi2_2_smooth < FASTEXP_MAX_CHI2/2) & (chi2_2_smooth >= 0)
+    exp_val_smooth[msk_smooth] = fexp_arr(-chi2_2_smooth[msk_smooth])
+    return exp_val_smooth
+
+
 @functools.lru_cache(maxsize=128)
 def _ksigma_kernels(
     dim,
@@ -537,11 +546,7 @@ def _ksigma_kernels(
 
     # add smoothing after norm check above for kernel size
     if fwhm_smooth > 0:
-        sigma_smooth = fwhm_to_sigma(fwhm_smooth)
-        chi2_2_smooth = sigma_smooth * sigma_smooth / 2 * fmag2
-        exp_val_smooth = np.zeros_like(karg)
-        msk_smooth = (chi2_2_smooth < FASTEXP_MAX_CHI2/2) & (chi2_2_smooth >= 0)
-        exp_val_smooth[msk_smooth] = fexp_arr(-chi2_2_smooth[msk_smooth])
+        exp_val_smooth = _get_fwhm_smooth_profile(fwhm_smooth, fmag2)
         fkf *= exp_val_smooth
         knrm *= exp_val_smooth
 
@@ -641,11 +646,7 @@ def _gauss_kernels(
 
     # add smoothing after norm check above for kernel size
     if fwhm_smooth > 0:
-        sigma_smooth = fwhm_to_sigma(fwhm_smooth)
-        chi2_2_smooth = sigma_smooth * sigma_smooth / 2 * fmag2
-        exp_val_smooth = np.zeros_like(exp_val)
-        msk_smooth = (chi2_2_smooth < FASTEXP_MAX_CHI2/2) & (chi2_2_smooth >= 0)
-        exp_val_smooth[msk_smooth] = fexp_arr(-chi2_2_smooth[msk_smooth])
+        exp_val_smooth = _get_fwhm_smooth_profile(fwhm_smooth, fmag2)
         fkf *= exp_val_smooth
         knrm *= exp_val_smooth
 

From fbd59eacf63cd43435285b078d9657f8fde78f76 Mon Sep 17 00:00:00 2001
From: beckermr <becker.mr@gmail.com>
Date: Thu, 30 Jun 2022 06:27:30 -0500
Subject: [PATCH 12/15] ENH always clear the cache

---
 ngmix/tests/test_metacal_cache.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/ngmix/tests/test_metacal_cache.py b/ngmix/tests/test_metacal_cache.py
index 4d35b791..b50a0d69 100644
--- a/ngmix/tests/test_metacal_cache.py
+++ b/ngmix/tests/test_metacal_cache.py
@@ -11,6 +11,8 @@
 
 @flaky(max_runs=10)
 def test_metacal_cache():
+    _cached_galsim_stuff.cache_clear()
+
     # first warm up numba
     rng = np.random.RandomState(seed=100)
     obs = _get_obs(rng, noise=0.005, set_noise_image=True, psf_fwhm=0.8, n=300)

From 6772a4286000ad4fa12952a18cc86806b24c8d37 Mon Sep 17 00:00:00 2001
From: beckermr <becker.mr@gmail.com>
Date: Thu, 30 Jun 2022 06:28:30 -0500
Subject: [PATCH 13/15] BUG if numba already warm then this assert is not true

---
 ngmix/tests/test_metacal_cache.py | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/ngmix/tests/test_metacal_cache.py b/ngmix/tests/test_metacal_cache.py
index b50a0d69..9e2e641d 100644
--- a/ngmix/tests/test_metacal_cache.py
+++ b/ngmix/tests/test_metacal_cache.py
@@ -40,8 +40,5 @@ def test_metacal_cache():
     print("third time: %r seconds (< %r?)" % (t2, t1*0.7), flush=True)
     print(_cached_galsim_stuff.cache_info(), flush=True)
 
-    # numba should be slower always but we do not care how much
-    assert t1 < t0
-
     # we expect roughly 30% gains
     assert t2 < t1*0.7

From abd66589d5e3ff4d699698825f0193c81664ede8 Mon Sep 17 00:00:00 2001
From: beckermr <becker.mr@gmail.com>
Date: Thu, 30 Jun 2022 06:49:17 -0500
Subject: [PATCH 14/15] TST remove round test

---
 ngmix/tests/test_prepsfmom.py | 9 +++------
 1 file changed, 3 insertions(+), 6 deletions(-)

diff --git a/ngmix/tests/test_prepsfmom.py b/ngmix/tests/test_prepsfmom.py
index 84b1f2aa..bfd061f0 100644
--- a/ngmix/tests/test_prepsfmom.py
+++ b/ngmix/tests/test_prepsfmom.py
@@ -393,11 +393,10 @@ def test_prepsfmom_gauss(
 @pytest.mark.parametrize('fwhm,psf_fwhm', [(2.0, 1.0)])
 @pytest.mark.parametrize('image_size', [53])
 @pytest.mark.parametrize('pad_factor', [1.5])
-@pytest.mark.parametrize('round', [True, False])
 @pytest.mark.parametrize('fwhm_smooth', [0, 1])
 def test_prepsfmom_mn_cov_psf(
     pad_factor, image_size, fwhm, psf_fwhm, pixel_scale, snr, mom_fwhm, cls,
-    fwhm_smooth, round,
+    fwhm_smooth,
 ):
     """Slower test to make sure means and errors are right
     w/ tons of monte carlo samples.
@@ -434,11 +433,9 @@ def test_prepsfmom_mn_cov_psf(
     )
     psf = galsim.Gaussian(
         fwhm=psf_fwhm
+    ).shear(
+        g1=0.3, g2=-0.15
     )
-    if not round:
-        psf = psf.shear(
-            g1=0.3, g2=-0.15
-        )
     im = galsim.Convolve([gal, psf]).drawImage(
         nx=image_size,
         ny=image_size,

From 8851f4db7065a9f4982f6f8a7a9b45c546da3077 Mon Sep 17 00:00:00 2001
From: beckermr <becker.mr@gmail.com>
Date: Thu, 30 Jun 2022 06:54:45 -0500
Subject: [PATCH 15/15] TST added test for shape s/n with smoothing

---
 ngmix/tests/test_prepsfmom.py | 104 +++++++++++++++++++++++++++++++++-
 1 file changed, 103 insertions(+), 1 deletion(-)

diff --git a/ngmix/tests/test_prepsfmom.py b/ngmix/tests/test_prepsfmom.py
index bfd061f0..99718cfd 100644
--- a/ngmix/tests/test_prepsfmom.py
+++ b/ngmix/tests/test_prepsfmom.py
@@ -386,7 +386,7 @@ def test_prepsfmom_gauss(
 
 
 @pytest.mark.parametrize("cls,mom_fwhm,snr", [
-    # (KSigmaMom, 2.0, 1e2),
+    (KSigmaMom, 2.0, 1e2),
     (PGaussMom, 2.0, 1e2),
 ])
 @pytest.mark.parametrize('pixel_scale', [0.25])
@@ -528,6 +528,108 @@ def test_prepsfmom_mn_cov_psf(
     )
 
 
+@pytest.mark.parametrize("cls,mom_fwhm,snr", [(PGaussMom, 2.0, 1e2)])
+@pytest.mark.parametrize('pixel_scale', [0.25])
+@pytest.mark.parametrize('fwhm,psf_fwhm', [(2.0, 1.0)])
+@pytest.mark.parametrize('image_size', [53])
+@pytest.mark.parametrize('pad_factor', [1.5])
+def test_prepsfmom_fwhm_smooth_snr(
+    pad_factor, image_size, fwhm, psf_fwhm, pixel_scale, snr, mom_fwhm, cls,
+):
+    def _run_sim_fwhm_smooth(fwhm_smooth):
+        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)
+        psf_shift = rng.uniform(low=-scale/2, high=scale/2, size=2)
+        xy = gs_wcs.toImage(galsim.PositionD(shift))
+        psf_xy = gs_wcs.toImage(galsim.PositionD(psf_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)
+
+        psf_jac = Jacobian(
+            y=26 + psf_xy.y, x=26 + psf_xy.x,
+            dudx=gs_wcs.dudx, dudy=gs_wcs.dudy,
+            dvdx=gs_wcs.dvdx, dvdy=gs_wcs.dvdy)
+
+        gal = galsim.Gaussian(
+            fwhm=fwhm
+        ).shear(
+            g1=-0.1, g2=0.2
+        ).withFlux(
+            400
+        ).shift(
+            dx=shift[0], dy=shift[1]
+        )
+        psf = galsim.Gaussian(
+            fwhm=psf_fwhm
+        ).shear(
+            g1=0.3, g2=-0.15
+        )
+        im = galsim.Convolve([gal, psf]).drawImage(
+            nx=image_size,
+            ny=image_size,
+            wcs=gs_wcs
+        ).array
+        noise = np.sqrt(np.sum(im**2)) / snr
+        wgt = np.ones_like(im) / noise**2
+
+        psf_im = psf.shift(
+            dx=psf_shift[0], dy=psf_shift[1]
+        ).drawImage(
+            nx=53,
+            ny=53,
+            wcs=gs_wcs
+        ).array
+
+        fitter = cls(
+            fwhm=mom_fwhm,
+            pad_factor=pad_factor,
+            fwhm_smooth=fwhm_smooth,
+        )
+
+        # get true flux
+        im_true = gal.drawImage(
+            nx=image_size,
+            ny=image_size,
+            wcs=gs_wcs,
+            method='no_pixel').array
+        obs = Observation(
+            image=im_true,
+            jacobian=jac,
+        )
+        res = fitter.go(obs=obs, no_psf=True)
+
+        res = []
+        for _ in range(1_000):
+            _im = im + rng.normal(size=im.shape, scale=noise)
+            obs = Observation(
+                image=_im,
+                weight=wgt,
+                jacobian=jac,
+                psf=Observation(image=psf_im, jacobian=psf_jac),
+            )
+
+            _res = fitter.go(obs=obs)
+            if _res['flags'] == 0:
+                res.append(_res)
+
+        res = _stack_list_of_dicts(res)
+
+        return np.abs(np.mean(res["e"], axis=0))/np.std(res["e"], axis=0)
+
+    e_snr = _run_sim_fwhm_smooth(0)
+    e_snr_smooth = _run_sim_fwhm_smooth(1)
+
+    assert np.all(e_snr_smooth > e_snr)
+
+
 @pytest.mark.parametrize("cls,mom_fwhm,snr", [
     (PGaussMom, 2.0, 1e2),
     (KSigmaMom, 2.0, 1e2),