From fb1cd68a7e9ce2371e647802c9cb31cc8417cd9a Mon Sep 17 00:00:00 2001 From: Erin Sheldon Date: Tue, 3 Oct 2023 14:30:33 -0400 Subject: [PATCH 1/4] fix bug with metacal psf method "gauss" The psf determination was being done on the pre pixel rather than post pixel psf. This required a Convolve with the pixel afterward, which does not work as expected when the wcs is complex Added a metacal accuracy test as well --- CHANGES.md | 7 + ngmix/_version.py | 2 +- ngmix/metacal/metacal.py | 9 +- ngmix/tests/test_metacal_accuracy.py | 218 +++++++++++++++++++++++++++ 4 files changed, 231 insertions(+), 5 deletions(-) create mode 100644 ngmix/tests/test_metacal_accuracy.py diff --git a/CHANGES.md b/CHANGES.md index b8ee7777..9930ac26 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -1,3 +1,10 @@ +## v2.3.1 + +## Bug Fixes + + - The metacal psf reconvolution method 'gauss' (which uses MetacalGaussPSF) + now works with a complex WCS + ## v2.3.0 ### new features diff --git a/ngmix/_version.py b/ngmix/_version.py index 89ddae6d..2531adcc 100644 --- a/ngmix/_version.py +++ b/ngmix/_version.py @@ -1 +1 @@ -__version__ = '2.3.0' # noqa +__version__ = '2.3.1' # noqa diff --git a/ngmix/metacal/metacal.py b/ngmix/metacal/metacal.py index e892105a..dbf3a3dd 100644 --- a/ngmix/metacal/metacal.py +++ b/ngmix/metacal/metacal.py @@ -419,6 +419,7 @@ def _set_data(self): # interpolated psf deconvolved from pixel. This is what # we dilate, shear, etc and reconvolve the image by self.psf_int_nopix = galsim.Convolve([psf_int, self.pixel_inv]) + self.psf_int = psf_int def get_wcs(self): """ @@ -572,15 +573,15 @@ def _get_dilated_psf(self, shear, doshear=False): _do_dilate for the algorithm """ - import galsim + # import galsim assert doshear is False, 'no shearing gauss psf' gauss_psf = _get_gauss_target_psf( - self.psf_int_nopix, flux=self.psf_flux, + self.psf_int, + flux=self.psf_flux, ) - psf_grown_nopix = _do_dilate(gauss_psf, shear) - psf_grown = galsim.Convolve(psf_grown_nopix, self.pixel) + psf_grown = _do_dilate(gauss_psf, shear) return psf_grown def _make_psf_obs(self, gsim): diff --git a/ngmix/tests/test_metacal_accuracy.py b/ngmix/tests/test_metacal_accuracy.py new file mode 100644 index 00000000..2cb30b6e --- /dev/null +++ b/ngmix/tests/test_metacal_accuracy.py @@ -0,0 +1,218 @@ +import numpy as np +import ngmix +import galsim +import pytest + + +@pytest.mark.parametrize('psf', ['gauss', 'fitgauss', 'galsim_obj']) +def test_metacal_accuracy(psf): + + ntrial = 100 + seed = 99 + + if psf == 'galsim_obj': + psf = galsim.Gaussian(fwhm=1.05) + + # Wacky WCS + wcs_g1 = 0.1 + wcs_g2 = 0.0 + wcs = galsim.ShearWCS(0.263, galsim.Shear(g1=wcs_g1, g2=wcs_g2)).jacobian() + + print("WCS:", wcs.getDecomposition()) + print() + print() + + shear_true = 0.02 + rng = np.random.RandomState(seed) + + # We will measure moments with a fixed gaussian weight function + weight_fwhm = 1.2 + fitter = ngmix.gaussmom.GaussMom(fwhm=weight_fwhm) + psf_fitter = ngmix.gaussmom.GaussMom(fwhm=weight_fwhm) + + # these "runners" run the measurement code on observations + psf_runner = ngmix.runners.PSFRunner(fitter=psf_fitter) + runner = ngmix.runners.Runner(fitter=fitter) + + # this "bootstrapper" runs the metacal image shearing as well as both psf + # and object measurements + boot = ngmix.metacal.MetacalBootstrapper( + runner=runner, psf_runner=psf_runner, + rng=rng, + psf=psf, + types=['noshear', '1p', '1m'], + ) + dlist = [] + for i in range(ntrial): + im, psf_im, obs = _make_data(rng=rng, shear=shear_true, wcs=wcs) + resdict, obsdict = boot.go(obs) + for stype, sres in resdict.items(): + st = _make_struct(res=sres, obs=obsdict[stype], shear_type=stype) + dlist.append(st) + + print() + data = np.hstack(dlist) + + w = _select(data=data, shear_type='noshear') + w_1p = _select(data=data, shear_type='1p') + w_1m = _select(data=data, shear_type='1m') + + g1 = data['g'][w, 0].mean() + g1err = data['g'][w, 0].std() / np.sqrt(w.size) + g1_1p = data['g'][w_1p, 0].mean() + g1_1m = data['g'][w_1m, 0].mean() + + R11 = (g1_1p - g1_1m)/0.02 + + shear = g1 / R11 + shear_err = g1err / R11 + m = shear / shear_true - 1 + merr = shear_err / shear_true + + s2n = data['s2n'][w].mean() + print('S/N: %g' % s2n) + print('R11: %g' % R11) + print('m: %g +/- %g (99.7%% conf)' % (m, merr*3)) + + assert np.abs(m - 0.00034) < 1.0e-4 + + +def _make_struct(res, obs, shear_type): + """ + make the data structure + Parameters + ---------- + res: dict + With keys 's2n', 'e', and 'T' + obs: ngmix.Observation + The observation for this shear type + shear_type: str + The shear type + Returns + ------- + 1-element array with fields + """ + dt = [ + ('flags', 'i4'), + ('shear_type', 'U7'), + ('s2n', 'f8'), + ('g', 'f8', 2), + ('T', 'f8'), + ('Tpsf', 'f8'), + ] + data = np.zeros(1, dtype=dt) + data['shear_type'] = shear_type + data['flags'] = res['flags'] + if res['flags'] == 0: + data['s2n'] = res['s2n'] + # for moments we are actually measureing e, the elliptity + data['g'] = res['e'] + data['T'] = res['T'] + else: + data['s2n'] = np.nan + data['g'] = np.nan + data['T'] = np.nan + data['Tpsf'] = np.nan + # we only have one epoch and band, so we can get the psf T from the + # observation rather than averaging over epochs/bands + data['Tpsf'] = obs.psf.meta['result']['T'] + return data + + +def _select(data, shear_type): + """ + select the data by shear type and size + Parameters + ---------- + data: array + The array with fields shear_type and T + shear_type: str + e.g. 'noshear', '1p', etc. + Returns + ------- + array of indices + """ + # raw moments, so the T is the post-psf T. This the + # selection is > 1.2 rather than something smaller like 0.5 + # for pre-psf T from one of the maximum likelihood fitters + wtype, = np.where( + (data['shear_type'] == shear_type) & + (data['flags'] == 0) + ) + w, = np.where(data['T'][wtype]/data['Tpsf'][wtype] > 1.2) + print('%s kept: %d/%d' % (shear_type, w.size, wtype.size)) + w = wtype[w] + return w + + +def _make_data(rng, shear, wcs): + """ + simulate an exponential object with moffat psf + the hlr of the exponential is drawn from a gaussian + with mean 0.4 arcseconds and sigma 0.2 + Parameters + ---------- + rng: np.random.RandomState + The random number generator + noise: float + Noise for the image + shear: (g1, g2) + The shear in each component + Returns + ------- + ngmix.Observation + """ + noise = 1.0e-8 + psf_noise = 1.0e-8 + stamp_size = 91 + + psf_fwhm = 0.9 + gal_hlr = 0.5 + psf = galsim.Moffat( + beta=2.5, + fwhm=psf_fwhm, + ).shear( + g1=0.02, + g2=-0.01, + ) + obj0 = galsim.Exponential( + half_light_radius=gal_hlr, + ).shear( + g1=shear, + ) + obj = galsim.Convolve(psf, obj0) + + psf_im = psf.drawImage(nx=stamp_size, ny=stamp_size, wcs=wcs).array + im = obj.drawImage(nx=stamp_size, ny=stamp_size, wcs=wcs).array + + psf_im += rng.normal(scale=psf_noise, size=psf_im.shape) + im += rng.normal(scale=noise, size=im.shape) + + cen = (np.array(im.shape)-1.0)/2.0 + psf_cen = (np.array(psf_im.shape)-1.0)/2.0 + + jacobian = ngmix.Jacobian( + x=cen[1], y=cen[0], wcs=wcs.jacobian( + image_pos=galsim.PositionD(cen[1], cen[0]) + ), + ) + psf_jacobian = ngmix.Jacobian( + x=psf_cen[1], y=psf_cen[0], wcs=wcs.jacobian( + image_pos=galsim.PositionD(psf_cen[1], psf_cen[0]) + ), + ) + + wt = im*0 + 1.0/noise**2 + psf_wt = psf_im*0 + 1.0/psf_noise**2 + psf_obs = ngmix.Observation( + psf_im, + weight=psf_wt, + jacobian=psf_jacobian, + ) + obs = ngmix.Observation( + im, + weight=wt, + jacobian=jacobian, + psf=psf_obs, + ) + return im, psf_im, obs From 4dbca02abe9e1dc1f983a45355c8e81110edd1e8 Mon Sep 17 00:00:00 2001 From: Erin Sheldon Date: Wed, 4 Oct 2023 08:25:19 -0400 Subject: [PATCH 2/4] update test to include det_bands column --- 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 d202b523..7f2a298f 100644 --- a/mdet_tests/test_mdet_regression.py +++ b/mdet_tests/test_mdet_regression.py @@ -268,7 +268,7 @@ def test_mdet_regression(fname, write=False): ), } else: - assert col in ["shear", "shear_bands"] + assert col in ["det_bands", "shear", "shear_bands"] if __name__ == "__main__": From 8f97869817d1efc4592513f9ffda7f1a520c1ba0 Mon Sep 17 00:00:00 2001 From: Erin Sheldon Date: Wed, 4 Oct 2023 08:30:56 -0400 Subject: [PATCH 3/4] satisfy the flake8 Note sure why this wasn't failing before --- ngmix/observation.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/ngmix/observation.py b/ngmix/observation.py index f4653fae..e21b5806 100644 --- a/ngmix/observation.py +++ b/ngmix/observation.py @@ -463,8 +463,9 @@ def set_ormask(self, ormask): ormask = np.ascontiguousarray(ormask) assert len(ormask.shape) == 2, "ormask must be 2d" - assert (ormask.shape == image.shape),\ + assert (ormask.shape == image.shape), ( "image and ormask must be same shape" + ) self._ormask = ormask @@ -1006,8 +1007,9 @@ def append(self, obs_list): An ObsList. An AssertionError will be raised if `obs_list` is not an `ngmix.ObsList`. """ - assert isinstance(obs_list, ObsList),\ + assert isinstance(obs_list, ObsList), ( 'obs_list should be of type ObsList' + ) super(MultiBandObsList, self).append(obs_list) def get_s2n(self): @@ -1092,8 +1094,9 @@ def __setitem__(self, index, obs_list): """ over-riding this for type safety """ - assert isinstance(obs_list, ObsList),\ + assert isinstance(obs_list, ObsList), ( 'obs_list should be of type ObsList' + ) super(MultiBandObsList, self).__setitem__(index, obs_list) @@ -1306,8 +1309,9 @@ def __setitem__(self, index, kobs): """ over-riding this for type safety """ - assert isinstance(kobs, KObservation),\ + assert isinstance(kobs, KObservation), ( 'kobs should be of type KObservation' + ) super(KObsList, self).__setitem__(index, kobs) From 45c2dc36a183ab47924cc75c7e506dd914ec680e Mon Sep 17 00:00:00 2001 From: Erin Sheldon Date: Wed, 4 Oct 2023 09:03:38 -0400 Subject: [PATCH 4/4] update expected response test for fixed psf='gauss' metacal reconv method --- ngmix/tests/test_metacal_bootstrap.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/ngmix/tests/test_metacal_bootstrap.py b/ngmix/tests/test_metacal_bootstrap.py index a28d041b..d0fa6911 100644 --- a/ngmix/tests/test_metacal_bootstrap.py +++ b/ngmix/tests/test_metacal_bootstrap.py @@ -178,6 +178,4 @@ def test_metacal_bootstrap_gaussmom_response(metacal_caching): Rvals[i] = (res1p['e'][0] - res1m['e'][0])/0.02 Rmean = Rvals.mean() - # this response value comes from 2.0.0, and had an error less than of - # 1.0e-5 on it. Allow for differences in rng/arch etc. by taking 1.0e-4 - assert abs(Rmean - 0.28535) < 1.0e-4 + assert abs(Rmean - 0.28159) < 1.0e-4