Coverage for metadetect / lsst / tests / test_lsst_metadetect.py: 12%

1""" 

2test using lsst simple sim 

3""" 

4import sys 

5import numpy as np 

6import pytest 

7 

8import logging 

9import ngmix 

10import metadetect 

11from metadetect import procflags 

12from metadetect.lsst.metadetect import run_metadetect 

13from metadetect.lsst.measure import get_pgauss_fitter 

14from metadetect.lsst.configs import get_config 

15from metadetect.lsst import util 

16import lsst.afw.image as afw_image 

17from lsst.utils import getPackageDir 

18 

19try: 

20 getPackageDir('descwl_shear_sims') 

21 getPackageDir('descwl_coadd') 

22 skip_tests_on_simulations = False 

23except LookupError: 

24 skip_tests_on_simulations = True 

25 

26 

27logging.basicConfig( 

28 stream=sys.stdout, 

29 level=logging.INFO, 

30) 

31 

32 

33def make_lsst_sim(seed, mag=20, hlr=0.5, bands=None): 

34 import descwl_shear_sims 

35 

36 rng = np.random.RandomState(seed=seed) 

37 coadd_dim = 251 

38 

39 if bands is None: 

40 bands = ['i'] 

41 

42 galaxy_catalog = descwl_shear_sims.galaxies.FixedGalaxyCatalog( 

43 rng=rng, 

44 coadd_dim=coadd_dim, 

45 buff=20, 

46 layout='grid', 

47 mag=mag, 

48 hlr=hlr, 

49 ) 

50 

51 psf = descwl_shear_sims.psfs.make_fixed_psf(psf_type='gauss') 

52 

53 sim_data = descwl_shear_sims.make_sim( 

54 rng=rng, 

55 galaxy_catalog=galaxy_catalog, 

56 coadd_dim=coadd_dim, 

57 g1=0.02, 

58 g2=0.00, 

59 psf=psf, 

60 bands=bands, 

61 ) 

62 return sim_data 

63 

64 

65def do_coadding(rng, sim_data, nowarp): 

66 from descwl_coadd.coadd import make_coadd 

67 from descwl_coadd.coadd_nowarp import make_coadd_nowarp 

68 

69 bands = list(sim_data['band_data'].keys()) 

70 

71 if nowarp: 

72 coadd_data_list = [ 

73 make_coadd_nowarp( 

74 exp=sim_data['band_data'][band][0], 

75 psf_dims=sim_data['psf_dims'], 

76 rng=rng, 

77 remove_poisson=False, 

78 ) 

79 for band in bands 

80 ] 

81 else: 

82 coadd_data_list = [ 

83 make_coadd( 

84 exps=sim_data['band_data'][band], 

85 psf_dims=sim_data['psf_dims'], 

86 rng=rng, 

87 coadd_wcs=sim_data['coadd_wcs'], 

88 coadd_bbox=sim_data['coadd_bbox'], 

89 remove_poisson=False, 

90 ) 

91 for band in bands 

92 ] 

93 

94 return util.extract_multiband_coadd_data(coadd_data_list) 

95 

96 

97@pytest.mark.skipif( 

98 skip_tests_on_simulations, 

99 reason='descwl_shear_sims and/or descwl_coadd not available' 

100) 

101@pytest.mark.parametrize('subtract_sky', [None, False, True]) 

102@pytest.mark.parametrize("metacal_types_option", [None, "1p1m", "full"]) 

103def test_lsst_metadetect_smoke(subtract_sky, metacal_types_option): 

104 rng = np.random.RandomState(seed=116) 

105 

106 bands = ['r', 'i'] 

107 sim_data = make_lsst_sim(116, bands=bands) 

108 data = do_coadding(rng=rng, sim_data=sim_data, nowarp=True) 

109 

110 config = {} 

111 

112 if subtract_sky is not None: 

113 config['subtract_sky'] = subtract_sky 

114 

115 if metacal_types_option is not None: 

116 if metacal_types_option == "1p1m": 

117 metacal_types = ['noshear', '1p', '1m'] 

118 config['metacal'] = {} 

119 elif metacal_types_option == "full": 

120 metacal_types = ['noshear', '1p', '1m', '2p', '2m'] 

121 config['metacal'] = {} 

122 config['metacal']['types'] = metacal_types 

123 else: 

124 metacal_types = ['noshear', '1p', '1m'] 

125 

126 detected = afw_image.Mask.getPlaneBitMask('DETECTED') 

127 res = run_metadetect(rng=rng, config=config, **data) 

128 

129 # we remove the DETECTED bit 

130 assert np.all(res['noshear']['bmask'] & detected == 0) 

131 

132 for metacal_type in metacal_types: 

133 assert ( 

134 metacal_type in res.keys() 

135 ), f"metacal_type={metacal_type} not in res.keys()" 

136 

137 for front in ['gauss', 'pgauss']: 

138 if front == 'gauss': 

139 gname = f'{front}_g' 

140 assert gname in res['noshear'].dtype.names 

141 

142 flux_name = f'{front}_band_flux' 

143 

144 for shear in metacal_types: 

145 # 5x5 grid 

146 assert res[shear].size == 25 

147 

148 assert np.any(res[shear][f"{front}_flags"] == 0) 

149 assert np.all(res[shear]["mfrac"] == 0) 

150 

151 assert len(res[shear][flux_name].shape) == len(bands) 

152 assert len(res[shear][flux_name][0]) == len(bands) 

153 

154 

155@pytest.mark.skipif( 

156 skip_tests_on_simulations, 

157 reason='descwl_shear_sims and/or descwl_coadd not available' 

158) 

159def test_lsst_metadetect_shear_bands_missing(): 

160 rng = np.random.RandomState(seed=116) 

161 

162 bands = ['g', 'r', 'i', 'z'] 

163 sim_data = make_lsst_sim(116, bands=bands) 

164 data = do_coadding(rng=rng, sim_data=sim_data, nowarp=True) 

165 config = {"shear_bands": ["r", "Y"]} 

166 with pytest.raises(RuntimeError) as e: 

167 run_metadetect(rng=rng, config=config, **data) 

168 

169 assert "'r', 'Y'" in str(e.value) 

170 

171 

172@pytest.mark.skipif( 

173 skip_tests_on_simulations, 

174 reason='descwl_shear_sims and/or descwl_coadd not available' 

175) 

176def test_lsst_metadetect_shear_bands(): 

177 rng = np.random.RandomState(seed=116) 

178 

179 bands = ['g', 'r', 'i', 'z'] 

180 nband = len(bands) 

181 sim_data = make_lsst_sim(116, bands=bands) 

182 data = do_coadding(rng=rng, sim_data=sim_data, nowarp=True) 

183 

184 config = {"shear_bands": ["r", "z"]} 

185 metacal_types = ['noshear', '1p', '1m'] 

186 

187 detected = afw_image.Mask.getPlaneBitMask('DETECTED') 

188 res = run_metadetect(rng=rng, config=config, **data) 

189 

190 # we remove the DETECTED bit 

191 assert np.all(res['noshear']['bmask'] & detected == 0) 

192 

193 for metacal_type in metacal_types: 

194 assert ( 

195 metacal_type in res.keys() 

196 ), f"metacal_type={metacal_type} not in res.keys()" 

197 

198 for front in ['gauss', 'pgauss']: 

199 if front == 'gauss': 

200 gname = f'{front}_g' 

201 assert gname in res['noshear'].dtype.names 

202 

203 flux_name = f'{front}_band_flux' 

204 

205 for shear in metacal_types: 

206 # 5x5 grid 

207 assert res[shear].size == 25 

208 

209 assert np.any(res[shear][f"{front}_flags"] == 0) 

210 assert np.all(res[shear]["mfrac"] == 0) 

211 assert res[shear][flux_name].shape == (25, nband) 

212 

213 for shear in metacal_types: 

214 assert np.all(res[shear]["shear_bands"] == np.array([["13"]])) 

215 # g and i band should be all NaNs for gauss 

216 assert np.all(np.isnan(res[shear]["gauss_band_flux"][:, 0])) 

217 assert np.all(np.isnan(res[shear]["gauss_band_flux"][:, 2])) 

218 # rest should be finite 

219 assert np.all(np.isfinite(res[shear]["gauss_band_flux"][:, 1])) 

220 assert np.all(np.isfinite(res[shear]["gauss_band_flux"][:, 3])) 

221 assert np.all(np.isfinite(res[shear]["pgauss_band_flux"])) 

222 

223 

224@pytest.mark.skipif( 

225 skip_tests_on_simulations, 

226 reason='descwl_shear_sims and/or descwl_coadd not available' 

227) 

228def test_lsst_metadetect_pgauss(): 

229 rng = np.random.RandomState(seed=882) 

230 

231 bands = ['r', 'i'] 

232 sim_data = make_lsst_sim(116, bands=bands) 

233 data = do_coadding(rng=rng, sim_data=sim_data, nowarp=True) 

234 

235 fwhm = 2.0 

236 config = { 

237 'pgauss': { 

238 'fwhm': fwhm, 

239 } 

240 } 

241 

242 fitter = get_pgauss_fitter(config=get_config(config)) 

243 assert fitter.fwhm == fwhm 

244 

245 res = run_metadetect(rng=rng, config=config, **data) 

246 

247 for front in ['gauss', 'pgauss']: 

248 if front == 'gauss': 

249 gname = f'{front}_g' 

250 assert gname in res['noshear'].dtype.names 

251 

252 flux_name = f'{front}_band_flux' 

253 

254 for shear in ('noshear', '1p', '1m'): 

255 # 5x5 grid 

256 assert res[shear].size == 25 

257 

258 assert np.any(res[shear][f"{front}_flags"] == 0) 

259 assert np.all(res[shear]["mfrac"] == 0) 

260 

261 assert len(res[shear][flux_name].shape) == len(bands) 

262 assert len(res[shear][flux_name][0]) == len(bands) 

263 

264 

265@pytest.mark.skipif( 

266 skip_tests_on_simulations, 

267 reason='descwl_shear_sims and/or descwl_coadd not available' 

268) 

269def test_lsst_metadetect_fullcoadd_smoke(): 

270 rng = np.random.RandomState(seed=116) 

271 

272 bands = ['r', 'i'] 

273 sim_data = make_lsst_sim(882, bands=bands) 

274 data = do_coadding(rng=rng, sim_data=sim_data, nowarp=False) 

275 

276 config = {} 

277 res = run_metadetect(config=config, rng=rng, **data) 

278 

279 front = 'gauss' 

280 gname = f'{front}_g' 

281 flux_name = f'{front}_band_flux' 

282 assert gname in res['noshear'].dtype.names 

283 

284 for shear in ('noshear', '1p', '1m'): 

285 # 5x5 grid 

286 assert res[shear].size == 25 

287 

288 assert np.any(res[shear][f"{front}_flags"] == 0) 

289 assert np.all(res[shear]["mfrac"] == 0) 

290 

291 assert len(res[shear][flux_name].shape) == len(bands) 

292 assert len(res[shear][flux_name][0]) == len(bands) 

293 

294 

295@pytest.mark.skipif( 

296 skip_tests_on_simulations, 

297 reason='descwl_shear_sims and/or descwl_coadd not available' 

298) 

299def test_lsst_zero_weights(show=False): 

300 """ 

301 At time of writing, DM stack will still detect in regions with inf 

302 variance. Test this continues to be true. 

303 

304 However, we don't have detections in BRIGHT, see test 

305 test_lsst_masked_as_bright 

306 """ 

307 nobj = [] 

308 seed = 55 

309 for do_zero in [False, True]: 

310 rng = np.random.RandomState(seed) 

311 sim_data = make_lsst_sim(seed, mag=23) 

312 data = do_coadding(rng=rng, sim_data=sim_data, nowarp=False) 

313 

314 if do_zero: 

315 data['mbexp']['i'].variance.array[50:100, 50:100] = np.inf 

316 data['noise_mbexp']['i'].variance.array[50:100, 50:100] = np.inf 

317 

318 if show: 

319 import matplotlib.pyplot as mplt 

320 fig, axs = mplt.subplots(ncols=2) 

321 axs[0].imshow(data['mbexp']['i'].image.array) 

322 axs[1].imshow(data['mbexp']['i'].variance.array) 

323 mplt.show() 

324 

325 resdict = run_metadetect(rng=rng, config=None, **data) 

326 

327 if do_zero: 

328 for shear_type, tres in resdict.items(): 

329 w, = np.where( 

330 tres['stamp_flags'] & procflags.ZERO_WEIGHTS != 0 

331 ) 

332 assert w.size > 0, 'expected some stamp_flags set' 

333 assert np.all(tres['gauss_flags'][w] == procflags.NO_ATTEMPT) 

334 assert np.all(tres['pgauss_flags'][w] == procflags.NO_ATTEMPT) 

335 

336 else: 

337 for shear_type, tres in resdict.items(): 

338 # 5x5 grid 

339 assert tres.size == 25 

340 

341 nobj.append(resdict['noshear'].size) 

342 

343 assert nobj[0] == nobj[1] 

344 

345 

346@pytest.mark.skipif( 

347 skip_tests_on_simulations, 

348 reason='descwl_shear_sims and/or descwl_coadd not available' 

349) 

350def test_lsst_masked_as_bright(show=False): 

351 """ 

352 Make sure we don't detect in areas marked BRIGHT 

353 """ 

354 seed = 55 

355 afw_image.Mask.addMaskPlane('BRIGHT') 

356 bright = afw_image.Mask.getPlaneBitMask('BRIGHT') 

357 for do_zero in [False, True]: 

358 rng = np.random.RandomState(seed) 

359 sim_data = make_lsst_sim(seed, mag=23) 

360 data = do_coadding(rng=rng, sim_data=sim_data, nowarp=False) 

361 

362 if do_zero: 

363 data['mbexp']['i'].variance.array[50:100, 50:100] = np.inf 

364 data['mbexp']['i'].mask.array[50:100, 50:100] |= bright 

365 data['noise_mbexp']['i'].variance.array[50:100, 50:100] = np.inf 

366 data['noise_mbexp']['i'].mask.array[50:100, 50:100] |= bright 

367 

368 resdict = run_metadetect(rng=rng, config=None, **data) 

369 

370 if show: 

371 import matplotlib.pyplot as mplt 

372 fig, axs = mplt.subplots(ncols=2) 

373 axs[0].imshow(data['mbexp']['i'].image.array) 

374 axs[1].imshow(data['mbexp']['i'].variance.array) 

375 

376 axs[0].scatter( 

377 resdict['noshear']['col'] - resdict['noshear']['col0'], 

378 resdict['noshear']['row'] - resdict['noshear']['row0'], 

379 s=4, 

380 c='red', 

381 ) 

382 mplt.show() 

383 

384 if do_zero: 

385 for shear_type, tres in resdict.items(): 

386 assert tres.size == 24 

387 else: 

388 for shear_type, tres in resdict.items(): 

389 # 5x5 grid 

390 assert tres.size == 25 

391 

392 

393@pytest.mark.skipif( 

394 skip_tests_on_simulations, 

395 reason='descwl_shear_sims and/or descwl_coadd not available' 

396) 

397def test_lsst_metadetect_prepsf_stars(): 

398 seed = 55 

399 rng = np.random.RandomState(seed=seed) 

400 

401 sim_data = make_lsst_sim(seed, hlr=1.0e-4, mag=23) 

402 data = do_coadding(rng=rng, sim_data=sim_data, nowarp=True) 

403 

404 config = {} 

405 

406 res = run_metadetect(rng=rng, config=config, **data) 

407 

408 n = metadetect.util.Namer(front='pgauss') 

409 

410 data = res['noshear'] 

411 

412 wlowT, = np.where(data[n('flags')] != 0) 

413 wgood, = np.where(data[n('flags')] == 0) 

414 

415 # some will have T < 0 due to noise. Expect some with flags set 

416 assert wlowT.size > 0 

417 

418 assert np.any((data[n('flags')][wlowT] & ngmix.flags.NONPOS_SIZE) != 0) 

419 

420 for field in data.dtype.names: 

421 if field != "shear_bands": 

422 assert np.all(np.isfinite(data[field][wgood])), field 

423 

424 

425@pytest.mark.skipif( 

426 skip_tests_on_simulations, 

427 reason='descwl_shear_sims and/or descwl_coadd not available' 

428) 

429def test_lsst_metadetect_mfrac_ormask(show=False): 

430 rng = np.random.RandomState(seed=116) 

431 

432 ntrial = 1 

433 flag = 2**30 

434 

435 for trial in range(ntrial): 

436 sim_data = make_lsst_sim(rng.randint(0, 2**30)) 

437 data = do_coadding(rng=rng, sim_data=sim_data, nowarp=False) 

438 

439 data['mfrac_mbexp']['i'].image.array[:, :] = rng.uniform( 

440 size=data['mbexp']['i'].image.array.shape, low=0.2, high=0.8 

441 ) 

442 

443 for ormask in data['ormasks']: 

444 ormask[30:150, 30:150] = flag 

445 if show: 

446 import matplotlib.pyplot as mplt 

447 fig, axs = mplt.subplots(ncols=2) 

448 axs[0].imshow(data['mbexp']['i'].image.array) 

449 axs[1].imshow(ormask) 

450 mplt.show() 

451 

452 res = run_metadetect(config=None, rng=rng, **data) 

453 

454 for shear in ('noshear', '1p', '1m'): 

455 assert np.any(res[shear]["gauss_flags"] == 0) 

456 assert np.any( 

457 (res[shear]["mfrac"] > 0.40) 

458 & (res[shear]["mfrac"] < 0.60) 

459 ) 

460 assert np.any(res[shear]["ormask"] & flag != 0) 

461 

462 

463if __name__ == '__main__': 

464 test_lsst_masked_as_bright(show=True) 

465 # test_lsst_metadetect_smoke('wmom', 'False') 

466 # test_lsst_metadetect_mfrac_ormask(show=True)