From b9243ae51563ac6055129f05fe5119e3676cd266 Mon Sep 17 00:00:00 2001 From: Erin Sheldon Date: Fri, 14 Jun 2024 11:56:04 -0400 Subject: [PATCH 1/6] add rho4 to adaptive moments Add some docs as well and remove unused entry in admom dtype --- ngmix/admom/admom.py | 94 +++++++++++++++++++++++++++++++++++---- ngmix/admom/admom_nb.py | 14 ++++-- ngmix/gmix/gmix_nb.py | 13 +++++- ngmix/tests/test_admom.py | 11 ++++- 4 files changed, 118 insertions(+), 14 deletions(-) diff --git a/ngmix/admom/admom.py b/ngmix/admom/admom.py index ba2164fc..771d6726 100644 --- a/ngmix/admom/admom.py +++ b/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,6 +563,7 @@ 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 @@ -493,15 +571,15 @@ def get_result(ares, jac_area, wgt_norm): _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 = [ diff --git a/ngmix/admom/admom_nb.py b/ngmix/admom/admom_nb.py index febea928..79b4dae7 100644 --- a/ngmix/admom/admom_nb.py +++ b/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 diff --git a/ngmix/gmix/gmix_nb.py b/ngmix/gmix/gmix_nb.py index 25624832..0fb7700f 100644 --- a/ngmix/gmix/gmix_nb.py +++ b/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 diff --git a/ngmix/tests/test_admom.py b/ngmix/tests/test_admom.py index 4d77a95c..3abd244e 100644 --- a/ngmix/tests/test_admom.py +++ b/ngmix/tests/test_admom.py @@ -27,12 +27,14 @@ 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()) @@ -40,6 +42,7 @@ def test_admom(g1_true, g2_true, wcs_g1, wcs_g2): 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) From 929029e3b39ba8ded19855a1827f1c0508e68ece Mon Sep 17 00:00:00 2001 From: Erin Sheldon Date: Fri, 14 Jun 2024 11:58:47 -0400 Subject: [PATCH 2/6] update CI --- .github/workflows/test.yml | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index c66eb668..0f9c1277 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -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 \ From dda2bd318bc9f34767c8ca2f81ac5873f14e0587 Mon Sep 17 00:00:00 2001 From: Erin Sheldon Date: Fri, 14 Jun 2024 12:22:38 -0400 Subject: [PATCH 3/6] get inv from linalg rather than dual dual has been removed from numpy Update CI --- .coveragerc | 1 + .github/workflows/mdet.yml | 7 +++---- ngmix/fitting/leastsqbound.py | 2 +- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/.coveragerc b/.coveragerc index d0190664..047d473d 100644 --- a/.coveragerc +++ b/.coveragerc @@ -6,6 +6,7 @@ omit = ngmix/admom/*_nb.py ngmix/gmix/*_nb.py ngmix/pixels/*_nb.py + ngmix/tests/* [report] exclude_lines = diff --git a/.github/workflows/mdet.yml b/.github/workflows/mdet.yml index ef097360..dd37eae8 100644 --- a/.github/workflows/mdet.yml +++ b/.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 diff --git a/ngmix/fitting/leastsqbound.py b/ngmix/fitting/leastsqbound.py index 09ecb09b..cfee3831 100644 --- a/ngmix/fitting/leastsqbound.py +++ b/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, :]) From 2386770384e2a615009e3bdfab75f1e9f5ca6b68 Mon Sep 17 00:00:00 2001 From: Erin Sheldon Date: Fri, 14 Jun 2024 12:49:39 -0400 Subject: [PATCH 4/6] loosen up speed test also - use newer pythons - put tests back into coverage --- .coveragerc | 1 - .github/workflows/test.yml | 2 +- ngmix/tests/test_prepsfmom.py | 4 ++-- 3 files changed, 3 insertions(+), 4 deletions(-) diff --git a/.coveragerc b/.coveragerc index 047d473d..d0190664 100644 --- a/.coveragerc +++ b/.coveragerc @@ -6,7 +6,6 @@ omit = ngmix/admom/*_nb.py ngmix/gmix/*_nb.py ngmix/pixels/*_nb.py - ngmix/tests/* [report] exclude_lines = diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 0f9c1277..220e50f2 100644 --- a/.github/workflows/test.yml +++ b/.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" diff --git a/ngmix/tests/test_prepsfmom.py b/ngmix/tests/test_prepsfmom.py index 46c756f9..526faf4d 100644 --- a/ngmix/tests/test_prepsfmom.py +++ b/ngmix/tests/test_prepsfmom.py @@ -270,9 +270,9 @@ def test_prepsfmom_speed_and_cache(use_cache): # if numba stuff is cached this does not work so commented out # assert dt2 < dt1 if use_cache: - assert dt3/nfit < dt2*0.6 + assert dt3/nfit < dt2*0.8 else: - assert dt3/nfit >= dt2*0.6 + assert dt3/nfit >= dt2*0.8 def _stack_list_of_dicts(res): From 23b24076fd2ca698a2e44f2234ac0e01bb38f9b7 Mon Sep 17 00:00:00 2001 From: Erin Sheldon Date: Fri, 14 Jun 2024 13:01:00 -0400 Subject: [PATCH 5/6] add some "pragma: no cover" to tests --- ngmix/tests/test_em.py | 8 ++++---- ngmix/tests/test_leastsqbound.py | 2 +- ngmix/tests/test_metacal_accuracy.py | 2 +- ngmix/tests/test_prepsfmom.py | 3 ++- ngmix/tests/test_util.py | 4 ++-- 5 files changed, 10 insertions(+), 9 deletions(-) diff --git a/ngmix/tests/test_em.py b/ngmix/tests/test_em.py index 28f41142..6525151e 100644 --- a/ngmix/tests/test_em.py +++ b/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: diff --git a/ngmix/tests/test_leastsqbound.py b/ngmix/tests/test_leastsqbound.py index 7ddbe79a..4abfa661 100644 --- a/ngmix/tests/test_leastsqbound.py +++ b/ngmix/tests/test_leastsqbound.py @@ -105,7 +105,7 @@ def test_leastsqbound_smoke(use_prior): try: res = bootstrap(obs=obs, runner=runner, psf_runner=psf_runner) allflags[i] = res['flags'] - except ngmix.BootPSFFailure: + except ngmix.BootPSFFailure: # pragma: no cover allflags[i] = 1 assert np.any(allflags == 0) diff --git a/ngmix/tests/test_metacal_accuracy.py b/ngmix/tests/test_metacal_accuracy.py index 2cb30b6e..fc2b092b 100644 --- a/ngmix/tests/test_metacal_accuracy.py +++ b/ngmix/tests/test_metacal_accuracy.py @@ -77,7 +77,7 @@ def test_metacal_accuracy(psf): assert np.abs(m - 0.00034) < 1.0e-4 -def _make_struct(res, obs, shear_type): +def _make_struct(res, obs, shear_type): # pragma: no cover """ make the data structure Parameters diff --git a/ngmix/tests/test_prepsfmom.py b/ngmix/tests/test_prepsfmom.py index 526faf4d..ce81965c 100644 --- a/ngmix/tests/test_prepsfmom.py +++ b/ngmix/tests/test_prepsfmom.py @@ -275,7 +275,8 @@ def test_prepsfmom_speed_and_cache(use_cache): assert dt3/nfit >= dt2*0.8 -def _stack_list_of_dicts(res): +def _stack_list_of_dicts(res): # pragma: no cover + def _get_dtype(v): if isinstance(v, float): return ('f8',) diff --git a/ngmix/tests/test_util.py b/ngmix/tests/test_util.py index b6627aec..00d9e3df 100644 --- a/ngmix/tests/test_util.py +++ b/ngmix/tests/test_util.py @@ -46,5 +46,5 @@ def test_get_ratio_err_array(): if __name__ == '__main__': - test_get_ratio_err_scalar() - test_get_ratio_err_array() + test_get_ratio_err_scalar() # pragma: no cover + test_get_ratio_err_array() # pragma: no cover From 9f8f53d94303039c3baf6e151a52a4ea1d9bf9ee Mon Sep 17 00:00:00 2001 From: Erin Sheldon Date: Fri, 14 Jun 2024 13:23:40 -0400 Subject: [PATCH 6/6] update change log --- CHANGES.md | 4 ++++ ngmix/tests/_fakemeds.py | 2 +- 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/CHANGES.md b/CHANGES.md index 9930ac26..edab49c7 100644 --- a/CHANGES.md +++ b/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) diff --git a/ngmix/tests/_fakemeds.py b/ngmix/tests/_fakemeds.py index 5097b866..a1b4f7f0 100644 --- a/ngmix/tests/_fakemeds.py +++ b/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