Coverage for metadetect / lsst / tests / test_lsst_metadetect.py: 14%
293 statements
« prev ^ index » next coverage.py v7.12.0, created at 2025-12-17 16:57 -0800
1"""
2test using lsst simple sim
3"""
4import sys
5import numpy as np
6import pytest
8import logging
9import ngmix
10import metadetect
11from metadetect import procflags
12from metadetect.lsst.metadetect import run_metadetect, get_fitter
13from metadetect.lsst.configs import get_config
14from metadetect.lsst import util
15import lsst.afw.image as afw_image
16from lsst.utils import getPackageDir
18try:
19 getPackageDir('descwl_shear_sims')
20 getPackageDir('descwl_coadd')
21 skip_tests_on_simulations = False
22except LookupError:
23 skip_tests_on_simulations = True
26ngmix_v = float(ngmix.__version__[:3])
28logging.basicConfig(
29 stream=sys.stdout,
30 level=logging.INFO,
31)
34def make_lsst_sim(seed, mag=20, hlr=0.5, bands=None):
35 import descwl_shear_sims
37 rng = np.random.RandomState(seed=seed)
38 coadd_dim = 251
40 if bands is None:
41 bands = ['i']
43 galaxy_catalog = descwl_shear_sims.galaxies.FixedGalaxyCatalog(
44 rng=rng,
45 coadd_dim=coadd_dim,
46 buff=20,
47 layout='grid',
48 mag=mag,
49 hlr=hlr,
50 )
52 psf = descwl_shear_sims.psfs.make_fixed_psf(psf_type='gauss')
54 sim_data = descwl_shear_sims.make_sim(
55 rng=rng,
56 galaxy_catalog=galaxy_catalog,
57 coadd_dim=coadd_dim,
58 g1=0.02,
59 g2=0.00,
60 psf=psf,
61 bands=bands,
62 )
63 return sim_data
66def do_coadding(rng, sim_data, nowarp):
67 from descwl_coadd.coadd import make_coadd
68 from descwl_coadd.coadd_nowarp import make_coadd_nowarp
70 bands = list(sim_data['band_data'].keys())
72 if nowarp:
73 coadd_data_list = [
74 make_coadd_nowarp(
75 exp=sim_data['band_data'][band][0],
76 psf_dims=sim_data['psf_dims'],
77 rng=rng,
78 remove_poisson=False,
79 )
80 for band in bands
81 ]
82 else:
83 coadd_data_list = [
84 make_coadd(
85 exps=sim_data['band_data'][band],
86 psf_dims=sim_data['psf_dims'],
87 rng=rng,
88 coadd_wcs=sim_data['coadd_wcs'],
89 coadd_bbox=sim_data['coadd_bbox'],
90 remove_poisson=False,
91 )
92 for band in bands
93 ]
95 return util.extract_multiband_coadd_data(coadd_data_list)
98@pytest.mark.skipif(
99 skip_tests_on_simulations,
100 reason='descwl_shear_sims and/or descwl_coadd not available'
101)
102@pytest.mark.parametrize('meas_type', [None, 'wmom', 'ksigma', 'pgauss'])
103@pytest.mark.parametrize('subtract_sky', [None, False, True])
104@pytest.mark.parametrize("metacal_types_option", [None, "1p1m", "full"])
105def test_lsst_metadetect_smoke(meas_type, subtract_sky, metacal_types_option):
106 rng = np.random.RandomState(seed=116)
108 bands = ['r', 'i']
109 sim_data = make_lsst_sim(116, bands=bands)
110 data = do_coadding(rng=rng, sim_data=sim_data, nowarp=True)
112 config = {}
114 if subtract_sky is not None:
115 config['subtract_sky'] = subtract_sky
117 if meas_type is not None:
118 config['meas_type'] = meas_type
120 if metacal_types_option is not None:
121 if metacal_types_option == "1p1m":
122 metacal_types = ['noshear', '1p', '1m']
123 config['metacal'] = {}
124 elif metacal_types_option == "full":
125 metacal_types = ['noshear', '1p', '1m', '2p', '2m']
126 config['metacal'] = {}
127 config['metacal']['types'] = metacal_types
128 else:
129 metacal_types = ['noshear', '1p', '1m']
131 detected = afw_image.Mask.getPlaneBitMask('DETECTED')
132 res = run_metadetect(rng=rng, config=config, **data)
134 # we remove the DETECTED bit
135 assert np.all(res['noshear']['bmask'] & detected == 0)
137 for metacal_type in metacal_types:
138 assert (
139 metacal_type in res.keys()
140 ), f"metacal_type={metacal_type} not in res.keys()"
142 if meas_type is None:
143 front = 'wmom'
144 else:
145 front = meas_type
147 gname = f'{front}_g'
148 flux_name = f'{front}_band_flux'
149 assert gname in res['noshear'].dtype.names
151 for shear in metacal_types:
152 # 5x5 grid
153 assert res[shear].size == 25
155 assert np.any(res[shear][f"{front}_flags"] == 0)
156 assert np.all(res[shear]["mfrac"] == 0)
158 assert len(res[shear][flux_name].shape) == len(bands)
159 assert len(res[shear][flux_name][0]) == len(bands)
162@pytest.mark.skipif(
163 skip_tests_on_simulations,
164 reason='descwl_shear_sims and/or descwl_coadd not available'
165)
166@pytest.mark.skipif(ngmix_v < 2.1, reason="requires ngmix 2.1 or higher")
167@pytest.mark.parametrize('meas_type', ['wmom', 'ksigma', 'pgauss'])
168@pytest.mark.parametrize('fwhm_smooth', [None, 1.2])
169def test_lsst_metadetect_weight(meas_type, fwhm_smooth):
170 rng = np.random.RandomState(seed=882)
172 bands = ['r', 'i']
173 sim_data = make_lsst_sim(116, bands=bands)
174 data = do_coadding(rng=rng, sim_data=sim_data, nowarp=True)
176 fwhm = 2.0
177 config = {
178 'meas_type': meas_type,
179 'weight': {
180 'fwhm': fwhm,
181 }
182 }
183 if meas_type != 'wmom' and fwhm_smooth is not None:
184 config['weight']['fwhm_smooth'] = fwhm_smooth
186 fitter = get_fitter(config=get_config(config), rng=rng)
187 assert fitter.fwhm == fwhm
188 if meas_type != 'wmom' and fwhm_smooth is not None:
189 assert fitter.fwhm_smooth == fwhm_smooth
191 res = run_metadetect(rng=rng, config=config, **data)
193 gname = f'{meas_type}_g'
194 flux_name = f'{meas_type}_band_flux'
195 assert gname in res['noshear'].dtype.names
197 for shear in ('noshear', '1p', '1m'):
198 # 5x5 grid
199 assert res[shear].size == 25
201 assert np.any(res[shear][f"{meas_type}_flags"] == 0)
202 assert np.all(res[shear]["mfrac"] == 0)
204 assert len(res[shear][flux_name].shape) == len(bands)
205 assert len(res[shear][flux_name][0]) == len(bands)
208@pytest.mark.skipif(
209 skip_tests_on_simulations,
210 reason='descwl_shear_sims and/or descwl_coadd not available'
211)
212@pytest.mark.parametrize('fwhm_reg', [0, 0.8])
213def test_lsst_metadetect_fwhm_reg(fwhm_reg):
214 rng = np.random.RandomState(seed=882)
216 meas_type = 'pgauss'
217 bands = ['r', 'i']
218 sim_data = make_lsst_sim(35, bands=bands)
219 data = do_coadding(rng=rng, sim_data=sim_data, nowarp=True)
221 fwhm = 2.0
222 config = {
223 'meas_type': meas_type,
224 'weight': {
225 'fwhm': fwhm,
226 'fwhm_reg': fwhm_reg,
227 }
228 }
230 fitter = get_fitter(config=get_config(config), rng=rng)
231 assert fitter.fwhm == fwhm
233 res = run_metadetect(rng=rng, config=config, **data)
235 gname = f'{meas_type}_g'
236 flux_name = f'{meas_type}_band_flux'
237 assert gname in res['noshear'].dtype.names
239 for shear in ('noshear', '1p', '1m'):
240 # 5x5 grid
241 assert res[shear].size == 25
243 assert np.any(res[shear][f"{meas_type}_flags"] == 0)
244 assert np.all(res[shear]["mfrac"] == 0)
246 assert len(res[shear][flux_name].shape) == len(bands)
247 assert len(res[shear][flux_name][0]) == len(bands)
250@pytest.mark.skipif(
251 skip_tests_on_simulations,
252 reason='descwl_shear_sims and/or descwl_coadd not available'
253)
254def test_lsst_metadetect_fwhm_reg_shapenoise():
256 sdevs = {}
257 for fwhm_reg in [0, 0.8]:
258 meas_type = 'pgauss'
259 bands = ['r', 'i']
260 sim_data = make_lsst_sim(35, bands=bands)
261 rng = np.random.RandomState(seed=882)
262 data = do_coadding(rng=rng, sim_data=sim_data, nowarp=True)
264 fwhm = 2.0
265 config = {
266 'meas_type': meas_type,
267 'weight': {
268 'fwhm': fwhm,
269 'fwhm_reg': fwhm_reg,
270 }
271 }
273 fitter = get_fitter(config=get_config(config), rng=rng)
274 assert fitter.fwhm == fwhm
276 res = run_metadetect(rng=rng, config=config, **data)
278 w, = np.where(res['noshear'][f'{meas_type}_flags'] == 0)
280 sdevs[fwhm_reg] = res['noshear'][f'{meas_type}_g'][w, 0].std()
282 assert sdevs[0.8] < 0.6 * sdevs[0]
285@pytest.mark.skipif(
286 skip_tests_on_simulations,
287 reason='descwl_shear_sims and/or descwl_coadd not available'
288)
289def test_lsst_metadetect_am():
290 rng = np.random.RandomState(seed=882)
292 # only single band for am currently
293 bands = ['i']
294 sim_data = make_lsst_sim(116, bands=bands)
295 data = do_coadding(rng=rng, sim_data=sim_data, nowarp=True)
297 meas_type = 'am'
298 config = {'meas_type': meas_type}
300 res = run_metadetect(rng=rng, config=config, **data)
302 gname = f'{meas_type}_g'
303 flux_name = f'{meas_type}_band_flux'
304 assert gname in res['noshear'].dtype.names
306 for shear in ('noshear', '1p', '1m'):
307 # 5x5 grid
308 assert res[shear].size == 25
310 assert np.any(res[shear]["am_flags"] == 0)
311 assert np.all(res[shear]["mfrac"] == 0)
313 assert len(res[shear][flux_name].shape) == len(bands)
314 with pytest.raises(TypeError):
315 len(res[shear][flux_name][0])
318@pytest.mark.skipif(
319 skip_tests_on_simulations,
320 reason='descwl_shear_sims and/or descwl_coadd not available'
321)
322def test_lsst_metadetect_gauss():
323 rng = np.random.RandomState(seed=918)
325 bands = ['r', 'i', 'z']
326 sim_data = make_lsst_sim(116, bands=bands)
327 data = do_coadding(rng=rng, sim_data=sim_data, nowarp=True)
329 meas_type = 'gauss'
330 config = {'meas_type': meas_type}
332 res = run_metadetect(rng=rng, config=config, **data)
334 gname = f'{meas_type}_g'
335 flux_name = f'{meas_type}_band_flux'
336 assert gname in res['noshear'].dtype.names
338 for shear in ('noshear', '1p', '1m'):
339 # 5x5 grid of sources
340 assert res[shear].size == 25
342 assert np.all(res[shear][f"{meas_type}_flags"] == 0)
343 assert np.all(res[shear]["mfrac"] == 0)
345 assert res[shear][flux_name].shape[1] == len(bands)
346 assert np.all(np.isfinite(res[shear][flux_name]))
349@pytest.mark.skipif(
350 skip_tests_on_simulations,
351 reason='descwl_shear_sims and/or descwl_coadd not available'
352)
353def test_lsst_metadetect_fullcoadd_smoke():
354 rng = np.random.RandomState(seed=116)
356 bands = ['r', 'i']
357 sim_data = make_lsst_sim(882, bands=bands)
358 data = do_coadding(rng=rng, sim_data=sim_data, nowarp=False)
360 config = {'meas_type': 'pgauss'}
361 res = run_metadetect(config=config, rng=rng, **data)
363 front = 'pgauss'
364 gname = f'{front}_g'
365 flux_name = f'{front}_band_flux'
366 assert gname in res['noshear'].dtype.names
368 for shear in ('noshear', '1p', '1m'):
369 # 5x5 grid
370 assert res[shear].size == 25
372 assert np.any(res[shear]["pgauss_flags"] == 0)
373 assert np.all(res[shear]["mfrac"] == 0)
375 assert len(res[shear][flux_name].shape) == len(bands)
376 assert len(res[shear][flux_name][0]) == len(bands)
379@pytest.mark.skipif(
380 skip_tests_on_simulations,
381 reason='descwl_shear_sims and/or descwl_coadd not available'
382)
383def test_lsst_zero_weights(show=False):
384 """
385 At time of writing, DM stack will still detect in regions with inf
386 variance. Test this continues to be true.
388 However, we don't have detections in BRIGHT, see test
389 test_lsst_masked_as_bright
390 """
391 nobj = []
392 seed = 55
393 for do_zero in [False, True]:
394 rng = np.random.RandomState(seed)
395 sim_data = make_lsst_sim(seed, mag=23)
396 data = do_coadding(rng=rng, sim_data=sim_data, nowarp=False)
398 if do_zero:
399 data['mbexp']['i'].variance.array[50:100, 50:100] = np.inf
400 data['noise_mbexp']['i'].variance.array[50:100, 50:100] = np.inf
402 if show:
403 import matplotlib.pyplot as mplt
404 fig, axs = mplt.subplots(ncols=2)
405 axs[0].imshow(data['mbexp']['i'].image.array)
406 axs[1].imshow(data['mbexp']['i'].variance.array)
407 mplt.show()
409 resdict = run_metadetect(rng=rng, config=None, **data)
411 if do_zero:
412 for shear_type, tres in resdict.items():
413 assert np.any(tres['wmom_flags'] & procflags.ZERO_WEIGHTS != 0)
414 assert np.any(tres['wmom_psf_flags'] & procflags.NO_ATTEMPT != 0)
415 else:
416 for shear_type, tres in resdict.items():
417 # 5x5 grid
418 assert tres.size == 25
420 nobj.append(resdict['noshear'].size)
422 assert nobj[0] == nobj[1]
425@pytest.mark.skipif(
426 skip_tests_on_simulations,
427 reason='descwl_shear_sims and/or descwl_coadd not available'
428)
429def test_lsst_masked_as_bright(show=False):
430 """
431 Make sure we don't detect in areas marked BRIGHT
432 """
433 seed = 55
434 afw_image.Mask.addMaskPlane('BRIGHT')
435 bright = afw_image.Mask.getPlaneBitMask('BRIGHT')
436 for do_zero in [False, True]:
437 rng = np.random.RandomState(seed)
438 sim_data = make_lsst_sim(seed, mag=23)
439 data = do_coadding(rng=rng, sim_data=sim_data, nowarp=False)
441 if do_zero:
442 data['mbexp']['i'].variance.array[50:100, 50:100] = np.inf
443 data['mbexp']['i'].mask.array[50:100, 50:100] |= bright
444 data['noise_mbexp']['i'].variance.array[50:100, 50:100] = np.inf
445 data['noise_mbexp']['i'].mask.array[50:100, 50:100] |= bright
447 resdict = run_metadetect(rng=rng, config=None, **data)
449 if show:
450 import matplotlib.pyplot as mplt
451 fig, axs = mplt.subplots(ncols=2)
452 axs[0].imshow(data['mbexp']['i'].image.array)
453 axs[1].imshow(data['mbexp']['i'].variance.array)
455 axs[0].scatter(
456 resdict['noshear']['col'] - resdict['noshear']['col0'],
457 resdict['noshear']['row'] - resdict['noshear']['row0'],
458 s=4,
459 c='red',
460 )
461 mplt.show()
463 if do_zero:
464 for shear_type, tres in resdict.items():
465 assert tres.size == 24
466 else:
467 for shear_type, tres in resdict.items():
468 # 5x5 grid
469 assert tres.size == 25
472@pytest.mark.skipif(
473 skip_tests_on_simulations,
474 reason='descwl_shear_sims and/or descwl_coadd not available'
475)
476@pytest.mark.parametrize('meas_type', ['ksigma', 'pgauss'])
477def test_lsst_metadetect_prepsf_stars(meas_type):
478 seed = 55
479 rng = np.random.RandomState(seed=seed)
481 sim_data = make_lsst_sim(seed, hlr=1.0e-4, mag=23)
482 data = do_coadding(rng=rng, sim_data=sim_data, nowarp=True)
484 config = {'meas_type': meas_type}
486 res = run_metadetect(rng=rng, config=config, **data)
488 n = metadetect.util.Namer(front=meas_type)
490 data = res['noshear']
492 wlowT, = np.where(data[f'{meas_type}_flags'] != 0)
493 wgood, = np.where(data[f'{meas_type}_flags'] == 0)
495 # some will have T < 0 due to noise. Expect some with flags set
496 assert wlowT.size > 0
498 assert np.any((data[n('flags')][wlowT] & ngmix.flags.NONPOS_SIZE) != 0)
500 assert np.any(np.isnan(data[n('g')][wlowT]))
501 for field in data.dtype.names:
502 if field != "shear_bands":
503 assert np.all(np.isfinite(data[field][wgood])), field
506@pytest.mark.skipif(
507 skip_tests_on_simulations,
508 reason='descwl_shear_sims and/or descwl_coadd not available'
509)
510def test_lsst_metadetect_mfrac_ormask(show=False):
511 rng = np.random.RandomState(seed=116)
513 ntrial = 1
514 flag = 2**30
516 for trial in range(ntrial):
517 sim_data = make_lsst_sim(rng.randint(0, 2**30))
518 data = do_coadding(rng=rng, sim_data=sim_data, nowarp=False)
520 data['mfrac_mbexp']['i'].image.array[:, :] = rng.uniform(
521 size=data['mbexp']['i'].image.array.shape, low=0.2, high=0.8
522 )
524 for ormask in data['ormasks']:
525 ormask[30:150, 30:150] = flag
526 if show:
527 import matplotlib.pyplot as mplt
528 fig, axs = mplt.subplots(ncols=2)
529 axs[0].imshow(data['mbexp']['i'].image.array)
530 axs[1].imshow(ormask)
531 mplt.show()
533 res = run_metadetect(config=None, rng=rng, **data)
535 for shear in ('noshear', '1p', '1m'):
536 assert np.any(res[shear]["wmom_flags"] == 0)
537 assert np.any(
538 (res[shear]["mfrac"] > 0.40)
539 & (res[shear]["mfrac"] < 0.60)
540 )
541 assert np.any(res[shear]["ormask"] & flag != 0)
544if __name__ == '__main__':
545 # test_lsst_metadetect_am()
546 test_lsst_metadetect_gauss()
547 # test_lsst_masked_as_bright(show=True)
548 # test_lsst_metadetect_smoke('wmom', 'False')
549 # test_lsst_metadetect_mfrac_ormask(show=True)