Merge pull request #241 from esheldon/rho4

add rho4 to adaptive moments

Author: esheldon

.github/workflows/mdet.yml

@@ -21,13 +21,12 @@ jobs:
           channel-priority: strict
           show-channel-urls: true
           miniforge-version: latest
-          miniforge-variant: Mambaforge
 
       - name: configure conda and install code
         shell: bash -l {0}
         run: |
-          mamba config --set always_yes yes
-          mamba install --quiet \
+          conda config --set always_yes yes
+          conda install --quiet \
             pip \
             setuptools \
             numpy \
@@ -60,7 +59,7 @@ jobs:
           popd
 
           pip uninstall ngmix -y
-          mamba install ngmix==2.0.3 -y
+          conda install ngmix==2.0.3 -y
           pushd mdet_tests
           python test_mdet_regression.py v2.0.3
           popd

.github/workflows/test.yml

@@ -11,7 +11,7 @@ jobs:
     name: tests
     strategy:
       matrix:
-        pyver: ["3.8", "3.9", "3.10"]
+        pyver: ["3.9", "3.10", "3.11"]
 
     runs-on: "ubuntu-latest"
 
@@ -25,13 +25,12 @@ jobs:
           channel-priority: strict
           show-channel-urls: true
           miniforge-version: latest
-          miniforge-variant: Mambaforge
 
       - name: configure conda and install code
         shell: bash -l {0}
         run: |
-          mamba config --set always_yes yes
-          mamba install --quiet \
+          conda config --set always_yes yes
+          conda install --quiet \
             pip \
             setuptools \
             numpy \

CHANGES.md

@@ -1,5 +1,9 @@
 ## v2.3.1
 
+### new features
+
+    - Add calculation of "rho4" to the adaptive moments code.
+
 ## Bug Fixes
 
     - The metacal psf reconvolution method 'gauss' (which uses MetacalGaussPSF)

ngmix/admom/admom.py

@@ -167,14 +167,66 @@ def find_cen_admom(
 class AdmomResult(dict):
     """
     Represent a fit using adaptive moments, and generate images and mixtures
-    for the best fit
+    for the best fit.  This inherits from dict and has entries for
+    the results of adaptive moments.
 
     Parameters
     ----------
     obs: observation(s)
         Observation, ObsList, or MultiBandObsList
     result: dict
-        the basic fit result, to bad added to this object's keys
+        the basic fit result, to be added to this object's keys
+
+    The entries will include, but not be limited to
+    -----------------------------------------------
+    flags: int
+        flags for processing
+    flagstr: str
+        Explanation of flags
+    numiter: int
+        number of iterations in adaptive moments algorithm
+    npix: int
+        Number of pixels used
+    flux_flags: int
+        Flags for flux
+    flux_flagstr: int
+        Explanation of flux flags
+    flux: float
+        Flux estimate
+    flux_err: float
+        Error on flux estimate
+    T_flags: int
+        Flags for T
+    T_flagstr: int
+        Explanation of T flags
+    T: float
+        T for gaussian
+    T_err: float
+        Error on T
+    e1: float
+        First ellipticity parameter
+    e1_err: float
+        Error on first ellipticity parameter
+    e2: float
+        Second ellipticity parameter
+    e2_err: float
+        Error on second ellipticity parameter
+    rho4_flags: int
+        Flags for rho4
+    rho4_flagstr: int
+        Explanation of rho4 flags
+    rho4: float
+        rho4 for gaussian
+    rho4_err: float
+        Error on rho4
+    pars: array
+        Array of gaussian pars, size 6 with [v, u, e1, e2, T, flux]
+    wsum: float
+        Sum of weights over pixels
+    sums: array
+        Array of size 7, holding sums for
+    sums_cov: array
+        Covariance sums, shape (7, 7)
     """
 
     def __init__(self, obs, result):
@@ -368,7 +420,7 @@ def get_result(ares, jac_area, wgt_norm):
         if n == 'sums':
             res[n] = ares[n].copy()
         elif n == 'sums_cov':
-            res[n] = ares[n].reshape((6, 6)).copy()
+            res[n] = ares[n].reshape((7, 7)).copy()
         else:
             res[n] = ares[n]
     res["sums_norm"] = ares["wsum"]
@@ -378,12 +430,16 @@ def get_result(ares, jac_area, wgt_norm):
     res["flux_flagstr"] = ""
     res["T_flags"] = 0
     res["T_flagstr"] = ""
+    res["rho4_flags"] = 0
+    res["rho4_flagstr"] = ""
 
     res['flux'] = np.nan
     res['flux_mean'] = np.nan
     res["flux_err"] = np.nan
     res["T"] = np.nan
     res["T_err"] = np.nan
+    res["rho4"] = np.nan
+    res["rho4_err"] = np.nan
     res["s2n"] = np.nan
     res["e1"] = np.nan
     res["e2"] = np.nan
@@ -396,6 +452,7 @@ def get_result(ares, jac_area, wgt_norm):
     # set things we always set if flags are ok
     if res['flags'] == 0:
         res['T'] = res['pars'][4]
+        res['rho4'] = ares['rho4']
         flux_sum = res['sums'][5]
         res['flux_mean'] = flux_sum/res['wsum']
         res['pars'][5] = res['flux_mean']
@@ -420,8 +477,8 @@ def get_result(ares, jac_area, wgt_norm):
     else:
         res['flux_flags'] |= res['flags']
 
-    # handle flux+T only
     if res['flags'] == 0:
+        # T
         if res['sums_cov'][4, 4] > 0 and res['sums_cov'][5, 5] > 0:
             if res['sums'][5] > 0:
                 # the sums include the weight, so need factor of two to correct
@@ -437,8 +494,28 @@ def get_result(ares, jac_area, wgt_norm):
                 res["T_flags"] |= ngmix.flags.NONPOS_FLUX
         else:
             res["T_flags"] |= ngmix.flags.NONPOS_VAR
+
+        # rho4
+        if res['sums_cov'][6, 6] > 0 and res['sums_cov'][5, 5] > 0:
+            if res['sums'][5] > 0:
+                res['rho4'] = res['sums'][6] / res['sums'][5]
+                # the sums include the weight, so need factor of two to correct
+                res['rho4_err'] = 4*get_ratio_error(
+                    res['sums'][6],
+                    res['sums'][5],
+                    res['sums_cov'][6, 6],
+                    res['sums_cov'][5, 5],
+                    res['sums_cov'][6, 5],
+                )
+            else:
+                # flux <= 0.0
+                res["rho4_flags"] |= ngmix.flags.NONPOS_FLUX
+        else:
+            res["rho4_flags"] |= ngmix.flags.NONPOS_VAR
+
     else:
         res['T_flags'] |= res['flags']
+        res['rho4_flags'] |= res['flags']
 
     # now handle full flags
     if not np.all(np.diagonal(res['sums_cov'][2:, 2:]) > 0):
@@ -486,22 +563,23 @@ def get_result(ares, jac_area, wgt_norm):
     res['flagstr'] = ngmix.flags.get_flags_str(res['flags'])
     res['flux_flagstr'] = ngmix.flags.get_flags_str(res['flux_flags'])
     res['T_flagstr'] = ngmix.flags.get_flags_str(res['T_flags'])
+    res['rho4_flagstr'] = ngmix.flags.get_flags_str(res['rho4_flags'])
 
     return res
 
 
 _admom_result_dtype = [
     ('flags', 'i4'),
     ('numiter', 'i4'),
-    ('nimage', 'i4'),
     ('npix', 'i4'),
     ('wsum', 'f8'),
 
-    ('sums', 'f8', 6),
-    ('sums_cov', 'f8', (6, 6)),
+    ('sums', 'f8', 7),
+    ('sums_cov', 'f8', (7, 7)),
     ('pars', 'f8', 6),
+    ('rho4', 'f8'),
     # temporary
-    ('F', 'f8', 6),
+    ('F', 'f8', 7),
 ]
 
 _admom_conf_dtype = [

ngmix/admom/admom_nb.py

@@ -86,6 +86,7 @@ def admom(confarray, wt, pixels, resarray):
             res['pars'][3] = 2.0*wt['irc'][0]
             res['pars'][4] = wt['icc'][0] + wt['irr'][0]
             res['pars'][5] = 1.0
+            res['rho4'] = res['sums'][6] / res['sums'][5]
 
             break
 
@@ -151,19 +152,26 @@ def admom_momsums(wt, pixels, res):
         wdata = weight*pixel['val']
         w2 = weight*weight
 
+        chi2 = (
+            wt["dcc"][0] * vmod * vmod
+            + wt["drr"][0] * umod * umod
+            - 2.0 * wt["drc"][0] * vmod * umod
+        )
+
         F[0] = pixel['v']
         F[1] = pixel['u']
         F[2] = umod*umod - vmod*vmod
         F[3] = 2*vmod*umod
         F[4] = umod*umod + vmod*vmod
         F[5] = 1.0
+        F[6] = chi2 * chi2
 
         res['wsum'] += weight
         res['npix'] += 1
 
-        for i in range(6):
+        for i in range(7):
             res['sums'][i] += wdata*F[i]
-            for j in range(6):
+            for j in range(7):
                 res['sums_cov'][i, j] += w2*var*F[i]*F[j]
 
 
@@ -228,8 +236,8 @@ def clear_result(res):
     res['sums'][:] = 0.0
     res['sums_cov'][:, :] = 0.0
     res['pars'][:] = np.nan
+    res['rho4'] = np.nan
 
     # res['flags']=0
     # res['numiter']=0
-    # res['nimage']=0
     # res['F'][:]=0.0

ngmix/fitting/leastsqbound.py

@@ -538,7 +538,7 @@ def wDfun(x, *args):  # wrapped Dfun
                              retval[1]['ipvt'] - 1)).T
         cov_x = None
         if info in [1, 2, 3, 4]:
-            from numpy.dual import inv
+            from numpy.linalg import inv
             from numpy.linalg import LinAlgError
             perm = take(eye(n), retval[1]['ipvt'] - 1, 0)
             r = triu(transpose(retval[1]['fjac'])[:n, :])

ngmix/gmix/gmix_nb.py

@@ -672,7 +672,18 @@ def g1g2_to_e1e2(g1, g2):
 @njit
 def get_weighted_sums(wt, pixels, res, maxrad):
     """
-    do sums for calculating the weighted moments
+    Do sums for calculating the weighted moments.
+
+    Parameters
+    ----------
+    wt: array
+        The gaussian mixture with dtype ngmix.gmix.gmix._gauss2d_dtype
+    pixels: array
+        Array of pixels
+    res: array
+        The result array
+    maxrad: float
+        Maximum radius in u, v coordinates
     """
 
     maxrad2 = maxrad ** 2

ngmix/tests/_fakemeds.py

@@ -34,7 +34,7 @@ def make_fake_meds(
     import fitsio
 
     if cutout_types is None:
-        cutout_types = copy.deepcopy(CUTOUT_TYPES)
+        cutout_types = copy.deepcopy(CUTOUT_TYPES)  # pragma: no cover
 
     with_psf = 'psf' in cutout_types
 

ngmix/tests/test_admom.py

@@ -27,19 +27,22 @@ def test_admom(g1_true, g2_true, wcs_g1, wcs_g2):
     ).shear(
         g1=g1_true, g2=g2_true
     ).withFlux(
-        400)
+        400,
+    )
     im = obj.drawImage(
         nx=image_size,
         ny=image_size,
         wcs=gs_wcs,
-        method='no_pixel').array
+        method='no_pixel',
+    ).array
     noise = np.sqrt(np.sum(im**2)) / 1e18
     wgt = np.ones_like(im) / noise**2
     scale = np.sqrt(gs_wcs.pixelArea())
 
     g1arr = []
     g2arr = []
     Tarr = []
+    rho4arr = []
     for _ in range(50):
         shift = rng.uniform(low=-scale/2, high=scale/2, size=2)
         xy = gs_wcs.toImage(galsim.PositionD(shift))
@@ -75,6 +78,7 @@ def test_admom(g1_true, g2_true, wcs_g1, wcs_g2):
             g1arr.append(_g1)
             g2arr.append(_g2)
             Tarr.append(_T)
+            rho4arr.append(res['rho4'])
 
             fim = res.make_image()
             assert fim.shape == im.shape
@@ -94,6 +98,9 @@ def test_admom(g1_true, g2_true, wcs_g1, wcs_g2):
     if g1_true == 0 and g2_true == 0:
         T = np.mean(Tarr)
         assert np.abs(T - fwhm_to_T(fwhm)) < 1e-6
+        rho4_mean = np.mean(rho4arr)
+        # gaussians have rho4 of 2.0
+        assert np.abs(rho4_mean - 2) < 1e-5
 
     with pytest.raises(ValueError):
         ngmix.admom.run_admom(None, None)

ngmix/tests/test_em.py

@@ -145,9 +145,9 @@ def test_em_2gauss(noise):
     f1 = pars[0]
     f2 = pars[6]
     if f1 > f2:
-        indices = [1, 0]
+        indices = [1, 0]  # pragma: no cover
     else:
-        indices = [0, 1]
+        indices = [0, 1]  # pragma: no cover
 
     # only check pars for no noise
     if noise == 0.0:
@@ -213,9 +213,9 @@ def test_em_2gauss_withpsf(noise):
     f1 = pars[0]
     f2 = pars[6]
     if f1 > f2:
-        indices = [1, 0]
+        indices = [1, 0]  # pragma: no cover
     else:
-        indices = [0, 1]
+        indices = [0, 1]  # pragma: no cover
 
     # only check pars for no noise
     if noise == 0.0: