Coverage for metadetect / metadetect.py: 6%

1import logging 

2import time 

3 

4import numpy as np 

5import ngmix 

6from ngmix.gexceptions import BootPSFFailure 

7import esutil as eu 

8 

9from . import fitting 

10from . import procflags 

11from . import shearpos 

12from .util import Namer 

13from .mfrac import measure_mfrac 

14from .fitting import ( 

15 fit_mbobs_list_wavg, 

16 combine_fit_res, 

17 fit_mbobs_list_joint, 

18 MAX_NUM_SHEAR_BANDS, 

19) 

20 

21logger = logging.getLogger(__name__) 

22 

23 

24def do_metadetect( 

25 config, mbobs, rng, shear_band_combs=None, 

26 color_key_func=None, color_dep_mbobs=None, 

27 det_band_combs=None, 

28): 

29 """Run metadetect on the multi-band observations. 

30 

31 Parameters 

32 ---------- 

33 config: dict 

34 Configuration dictionary. Possible entries are 

35 

36 metacal 

37 weight 

38 model 

39 

40 mbobs: ngmix.MultiBandObsList 

41 We will do detection and measurements on these images 

42 rng: numpy.random.RandomState 

43 Random number generator 

44 shear_band_combs: list of list of int, optional 

45 If given, each element of the outer list is a list of indices into mbobs to use 

46 for shear measurement. Shear measurements will be made for each element of the 

47 outer list. If None, then shear measurements will be made for all entries in 

48 mbobs. 

49 det_band_combs: list of list of int or str, optional 

50 If given, the set of bands to use for detection. The default of None uses all 

51 of the bands. If the string "shear_bands" is passed, the code uses the bands 

52 used for shear. 

53 color_key_func: function, optional 

54 If given, a function that computes a color or tuple of colors to key the 

55 `color_dep_mbobs` dictionary given an input set of fluxes from the mbobs. 

56 color_dep_mbobs: dict of mbobs, optional 

57 A dictionary of color-dependently rendered observations of the mbobs for use 

58 in color-dependent metadetect. 

59 

60 Returns 

61 ------- 

62 res: dict 

63 The fitting data keyed on the shear component. 

64 """ 

65 md = Metadetect( 

66 config, mbobs, rng, 

67 shear_band_combs=shear_band_combs, 

68 color_key_func=color_key_func, 

69 color_dep_mbobs=color_dep_mbobs, 

70 det_band_combs=det_band_combs, 

71 ) 

72 md.go() 

73 return md.result 

74 

75 

76class Metadetect(dict): 

77 """Metadetect fitter for multi-band observations. 

78 

79 This class runs the metadetect shear fitting algorithm for observations 

80 consisting of one image per band. It can use several different fitters for the 

81 actual shape measurement. 

82 

83 wmom - Perform a post-PSF weighted moments measurement. The shear measurement 

84 is a moment computed from the inverse variance weighted average across 

85 the bands. 

86 

87 ksigma - Perform a pre-PSF weighted moments measurement. The shear measurement 

88 is a moment computed from the inverse variance weighted average across 

89 the bands. 

90 

91 pgauss - Perform a pre-PSF weighted moments measurement with a Gaussian filter. 

92 The shear measurement is a moment computed from the inverse variance 

93 weighted average across the bands. 

94 

95 am or admom - Use adaptive moments. The shear measurement is compute from fitting 

96 adaptive moments on a coadd of the bands used for shear. 

97 

98 Parameters 

99 ---------- 

100 config: dict 

101 Configuration dictionary. Possible entries are 

102 

103 metacal 

104 weight 

105 model 

106 

107 mbobs: ngmix.MultiBandObsList 

108 We will do detection and measurements on these images 

109 rng: numpy.random.RandomState 

110 Random number generator 

111 show: bool, optional 

112 If True, show the images using descwl_coadd.vis. 

113 shear_band_combs: list of list of int, optional 

114 If given, each element of the outer list is a list of indices into mbobs to use 

115 for shear measurement. Shear measurements will be made for each element of the 

116 outer list. If None, then shear measurements will be made for all entries in 

117 mbobs. 

118 det_band_combs: list of list of int or str, optional 

119 If given, the set of bands to use for detection. The default of None uses all 

120 of the bands. If the string "shear_bands" is passed, the code uses the bands 

121 used for shear. 

122 color_key_func: function, optional 

123 If given, a function that computes a color or tuple of colors to key the 

124 `color_dep_mbobs` dictionary given an input set of fluxes from the mbobs. If 

125 this function returns None, the PSF from the original observation is used. 

126 color_dep_mbobs: dict of mbobs, optional 

127 A dictionary of color-dependently rendered observations of the mbobs for use 

128 in color-dependent metadetect. 

129 """ 

130 def __init__( 

131 self, config, mbobs, rng, show=False, 

132 shear_band_combs=None, 

133 color_key_func=None, 

134 color_dep_mbobs=None, 

135 det_band_combs=None, 

136 ): 

137 self._show = show 

138 

139 self._set_config(config) 

140 self.mbobs = mbobs 

141 self.nband = len(mbobs) 

142 self.rng = rng 

143 

144 self.color_key_func = color_key_func 

145 self.color_dep_mbobs = color_dep_mbobs 

146 if ( 

147 (color_dep_mbobs is not None and color_key_func is None) 

148 or 

149 (color_dep_mbobs is None and color_key_func is not None) 

150 ): 

151 raise RuntimeError( 

152 "You must both `color_dep_mbobs` and `color_key_func`!" 

153 ) 

154 

155 self._set_fitter() 

156 

157 if shear_band_combs is None: 

158 shear_band_combs = [ 

159 list(range(self.nband)), 

160 ] 

161 

162 self._shear_band_combs = shear_band_combs 

163 

164 if det_band_combs is None: 

165 det_band_combs = ( 

166 [list(range(self.nband))] 

167 * len(self._shear_band_combs) 

168 ) 

169 elif det_band_combs == "shear_bands": 

170 det_band_combs = self._shear_band_combs 

171 

172 self._det_band_combs = det_band_combs 

173 

174 def _set_config(self, config): 

175 """ 

176 set the config, dealing with defaults 

177 """ 

178 

179 self.update(config) 

180 assert 'metacal' in self, \ 

181 'metacal setting must be present in config' 

182 assert 'meds' in self, \ 

183 'meds setting must be present in config' 

184 self['nodet_flags'] = self.get('nodet_flags', 0) 

185 

186 def _get_ormask_and_bmask(self, mbobs): 

187 """ 

188 set the ormask and bmask, ored from all epochs 

189 """ 

190 

191 for band, obslist in enumerate(mbobs): 

192 nepoch = len(obslist) 

193 assert nepoch == 1, 'expected 1 epoch, got %d' % nepoch 

194 

195 obs = obslist[0] 

196 

197 if band == 0: 

198 ormask = obs.ormask.copy() 

199 bmask = obs.bmask.copy() 

200 else: 

201 ormask |= obs.ormask 

202 bmask |= obs.bmask 

203 

204 return ormask, bmask 

205 

206 def _get_mfrac(self, mbobs): 

207 """ 

208 get the masked fraction image, averaged over all bands 

209 """ 

210 wgts = [] 

211 mfrac = np.zeros_like(mbobs[0][0].image) 

212 for band, obslist in enumerate(mbobs): 

213 nepoch = len(obslist) 

214 assert nepoch == 1, 'expected 1 epoch, got %d' % nepoch 

215 

216 obs = obslist[0] 

217 msk = obs.weight > 0 

218 if not np.any(msk): 

219 wgt = 0 

220 else: 

221 wgt = np.median(obs.weight[msk]) 

222 if hasattr(obs, "mfrac"): 

223 mfrac += (obs.mfrac * wgt) 

224 wgts.append(wgt) 

225 

226 if np.sum(wgts) > 0: 

227 mfrac = mfrac / np.sum(wgts) 

228 else: 

229 mfrac[:, :] = 1.0 

230 

231 return mfrac 

232 

233 def _set_fitter(self): 

234 """ 

235 set the fitter to be used 

236 """ 

237 

238 def _get_fitter(cfg): 

239 model = cfg.get('model', 'wmom') 

240 symmetrize = cfg.get("symmetrize", True) 

241 

242 if "fwhm_smooth" in cfg.get("weight", {}): 

243 kwargs = {"fwhm_smooth": cfg["weight"]["fwhm_smooth"]} 

244 else: 

245 kwargs = {} 

246 

247 if model == 'wmom': 

248 fitter = ngmix.gaussmom.GaussMom(fwhm=cfg["weight"]["fwhm"]) 

249 is_wavg = True 

250 coadd = False 

251 elif model == 'ksigma': 

252 fitter = ngmix.prepsfmom.KSigmaMom( 

253 fwhm=cfg["weight"]["fwhm"], 

254 **kwargs, 

255 ) 

256 is_wavg = True 

257 coadd = False 

258 elif model == "pgauss": 

259 fitter = ngmix.prepsfmom.PGaussMom( 

260 fwhm=cfg["weight"]["fwhm"], 

261 **kwargs, 

262 ) 

263 is_wavg = True 

264 coadd = False 

265 elif model in ["admom", "am", "gauss"]: 

266 # we pass the name to our codes 

267 fitter = model 

268 is_wavg = False 

269 

270 # we set this defualt 

271 # it may be used to set the masked fraction measurement 

272 # aperture 

273 if "weight" not in cfg: 

274 cfg["weight"] = {} 

275 if "fwhm" not in cfg["weight"]: 

276 cfg["weight"]["fwhm"] = 1.2 

277 

278 coadd = cfg.get("coadd", False) 

279 else: 

280 raise ValueError("bad model: '%s'" % model) 

281 

282 if "fwhm_reg" in cfg.get("weight", {}): 

283 fwhm_reg = cfg["weight"]["fwhm_reg"] 

284 fitter.kind = fitter.kind + "_reg%0.2f" % cfg["weight"]["fwhm_reg"] 

285 else: 

286 fwhm_reg = 0 

287 

288 return ( 

289 model, fitter, cfg["weight"]["fwhm"], fwhm_reg, 

290 is_wavg, symmetrize, coadd, 

291 ) 

292 

293 if "fitters" in self and ( 

294 "model" in self 

295 or "weight" in self 

296 or "symmetrize" in self 

297 or "coadd" in self 

298 ): 

299 raise RuntimeError( 

300 "You can only specify one of fitters or " 

301 "model+weight+symmetrize+coadd!" 

302 ) 

303 

304 if "fitters" in self: 

305 fitters = [] 

306 fwhms = [] 

307 fwhm_regs = [] 

308 fitter_is_wavg = [] 

309 fitter_symmetrize = [] 

310 fitter_coadd = [] 

311 for fitter_cfg in self["fitters"]: 

312 _, fitter, fwhm, fwhm_reg, is_wavg, symmetrize, coadd = _get_fitter( 

313 fitter_cfg 

314 ) 

315 fitters.append(fitter) 

316 fwhms.append(fwhm) 

317 fwhm_regs.append(fwhm_reg) 

318 fitter_is_wavg.append(is_wavg) 

319 fitter_symmetrize.append(symmetrize) 

320 fitter_coadd.append(coadd) 

321 self._fitters = fitters 

322 self._fwhms = fwhms 

323 self._fwhm_regs = fwhm_regs 

324 self._fitter_is_wavg = fitter_is_wavg 

325 self._fitter_symmetrize = fitter_symmetrize 

326 self._fitter_coadd = fitter_coadd 

327 else: 

328 _, fitter, fwhm, fwhm_reg, is_wavg, symmetrize, coadd = _get_fitter(self) 

329 self._fitters = [fitter] 

330 self._fwhms = [fwhm] 

331 self._fwhm_regs = [fwhm_reg] 

332 self._fitter_is_wavg = [is_wavg] 

333 self._fitter_symmetrize = [symmetrize] 

334 self._fitter_coadd = [coadd] 

335 

336 @property 

337 def result(self): 

338 """ 

339 get the result dict, keyed by the metacal type such 

340 as 'noshear', '1p', '1m', '2p', '2m' 

341 """ 

342 if not hasattr(self, '_result'): 

343 raise RuntimeError('run go() first') 

344 

345 return self._result 

346 

347 def go(self): 

348 """Run metadetect and set the result.""" 

349 

350 mfrac = self._get_mfrac(self.mbobs) 

351 any_all_zero_weight = False 

352 any_all_masked = False 

353 for obsl in self.mbobs: 

354 for obs in obsl: 

355 if np.all(obs.weight == 0): 

356 any_all_zero_weight = True 

357 

358 if np.all((obs.bmask & self['nodet_flags']) != 0): 

359 any_all_masked = True 

360 

361 # if the there are no pixels with mfrac < 1 or it is all zero weight 

362 # or it is all masked, we cannot measure anything so set result to None 

363 if (not np.any(mfrac < 1)) or any_all_zero_weight or any_all_masked: 

364 self._result = None 

365 return 

366 

367 # we do metacal on everything so we can get fluxes for non-shear bands later 

368 mcal_res = self._get_mbobs_data(None, list(range(self.nband)))["mcal_res"] 

369 

370 if mcal_res is None: 

371 self._result = None 

372 return 

373 

374 # past this point, the code should always return a dictionary with the minimal 

375 # metacal types 

376 # this indicates that a measurement should have been possible 

377 # we may find nothing, but that is a different thing 

378 all_res = {} 

379 for shear_bands, det_bands in zip( 

380 self._shear_band_combs, self._det_band_combs 

381 ): 

382 if self.color_key_func is not None and self.color_dep_mbobs is not None: 

383 res = self._go_bands_with_color(shear_bands, mcal_res, det_bands) 

384 else: 

385 res = self._go_bands(shear_bands, mcal_res, det_bands) 

386 if res is not None: 

387 for k, v in res.items(): 

388 if v is None: 

389 continue 

390 if k not in all_res: 

391 all_res[k] = [] 

392 all_res[k].append(v) 

393 

394 for k in all_res: 

395 all_res[k] = np.hstack(all_res[k]) 

396 

397 for mcal_type in self['metacal'].get( 

398 "types", ngmix.metacal.METACAL_MINIMAL_TYPES 

399 ): 

400 if mcal_type not in all_res: 

401 # this can happen if we do not detect sources for one of the 

402 # metacal images 

403 all_res[mcal_type] = None 

404 

405 self._result = all_res 

406 

407 def _go_bands(self, shear_bands, mcal_res, det_bands): 

408 kdata = self._get_mbobs_data(None, shear_bands) 

409 

410 _result = {} 

411 for shear_str, shear_mbobs in mcal_res.items(): 

412 cat, mbobs_list = self._do_detect( 

413 shear_mbobs, 

414 det_bands, 

415 ) 

416 _result[shear_str] = self._measure( 

417 mbobs_list=mbobs_list, 

418 shear_bands=shear_bands, 

419 det_bands=det_bands, 

420 cat=cat, 

421 shear_str=shear_str, 

422 mfrac=kdata["mfrac"], 

423 bmask=kdata["bmask"], 

424 ormask=kdata["ormask"], 

425 psf_stats=kdata["psf_stats"], 

426 ) 

427 

428 return _result 

429 

430 def _go_bands_with_color(self, shear_bands, mcal_res, det_bands): 

431 from . import detect 

432 

433 _result = {} 

434 for shear_str, shear_mbobs in mcal_res.items(): 

435 if not self._fitter_is_wavg[0]: 

436 raise RuntimeError( 

437 "Color-dependent metadetect can only run if first" 

438 " fitter is one of wmom, pgauss, or ksigma!" 

439 ) 

440 

441 # we first detect and get color of each detection 

442 cat, mbobs_list = self._do_detect(shear_mbobs, det_bands) 

443 nocolor_data = fit_mbobs_list_wavg( 

444 mbobs_list=mbobs_list, 

445 fitter=self._fitters[0], 

446 shear_bands=shear_bands, 

447 bmask_flags=self.get("bmask_flags", 0), 

448 ) 

449 if nocolor_data is None: 

450 _result[shear_str] = None 

451 continue 

452 

453 # now we map color to the mbobs for that color 

454 n = Namer(self._fitters[0].kind) 

455 col = n("band_flux") 

456 color_keys = [ 

457 self.color_key_func(nocolor_data[col][i]) 

458 for i in range(nocolor_data.shape[0]) 

459 ] 

460 

461 # now we remeasure the object at that mbobs 

462 color_data = [] 

463 for i, color_key in enumerate(color_keys): 

464 kdata = self._get_mbobs_data(color_key, shear_bands) 

465 if kdata["mcal_res"] is None or kdata["mcal_res"][shear_str] is None: 

466 continue 

467 

468 _medsifier = detect.CatalogMEDSifier( 

469 kdata["mcal_res"][shear_str], 

470 cat['x'][i:i+1], 

471 cat['y'][i:i+1], 

472 cat['box_size'][i:i+1], 

473 ) 

474 mbm = _medsifier.get_multiband_meds() 

475 mbobs_list = mbm.get_mbobs_list( 

476 weight_type=self["meds"].get("weight_type", "weight"), 

477 ) 

478 

479 _data = self._measure( 

480 mbobs_list=mbobs_list, 

481 shear_bands=shear_bands, 

482 det_bands=det_bands, 

483 cat=cat[i:i+1], 

484 shear_str=shear_str, 

485 mfrac=kdata["mfrac"], 

486 bmask=kdata["bmask"], 

487 ormask=kdata["ormask"], 

488 psf_stats=kdata["psf_stats"], 

489 ) 

490 if _data is not None: 

491 color_data.append(_data) 

492 

493 if len(color_data) > 0: 

494 _result[shear_str] = np.hstack(color_data) 

495 else: 

496 _result[shear_str] = None 

497 

498 return _result 

499 

500 def _get_mbobs_data(self, key, shear_bands): 

501 logger.info("computing mbobs data: %s %s", key, shear_bands) 

502 

503 if key is None: 

504 mbobs = self.mbobs 

505 else: 

506 mbobs = self.color_dep_mbobs[key] 

507 # the caching will miss cases where the input mbobs is passed as one 

508 # of the color dependent ones as well 

509 # thus we test if they are the same object and set the key to None if 

510 # this happens 

511 if mbobs is self.mbobs: 

512 key = None 

513 

514 if not hasattr(self, "_mbobs_data_cache"): 

515 logger.info("set mbobs data caches") 

516 self._mbobs_data_cache = {} 

517 self._mcalpsf_data_cache = {} 

518 

519 if key not in self._mbobs_data_cache: 

520 logger.info("key %s not in mbobs data cache", key) 

521 self._mbobs_data_cache[key] = {} 

522 self._mcalpsf_data_cache[key] = {} 

523 

524 t0 = time.time() 

525 try: 

526 fitting.fit_all_psfs(mbobs, self.rng) 

527 _psf_fit_flags = 0 

528 except BootPSFFailure: 

529 _psf_fit_flags = procflags.PSF_FAILURE 

530 self._mcalpsf_data_cache[key]["psf_fit_flags"] = _psf_fit_flags 

531 logger.info("PSF fits took %s seconds", time.time() - t0) 

532 

533 mcal_res = self._get_all_metacal(mbobs) 

534 self._mcalpsf_data_cache[key]["mcal_res"] = mcal_res 

535 

536 sbkey = tuple(sorted(shear_bands)) 

537 if sbkey not in self._mbobs_data_cache[key]: 

538 logger.info("shear_bands key %s not in mbobs data cache", sbkey) 

539 self._mbobs_data_cache[key][sbkey] = {} 

540 self._mbobs_data_cache[key][sbkey].update( 

541 self._mcalpsf_data_cache[key] 

542 ) 

543 

544 _mbobs = ngmix.MultiBandObsList() 

545 for band in shear_bands: 

546 _mbobs.append(mbobs[band]) 

547 mfrac = self._get_mfrac(_mbobs) 

548 ormask, bmask = self._get_ormask_and_bmask(_mbobs) 

549 psf_stats = _get_psf_stats( 

550 _mbobs, 

551 self._mbobs_data_cache[key][sbkey]["psf_fit_flags"], 

552 ) 

553 self._mbobs_data_cache[key][sbkey]["mfrac"] = mfrac 

554 self._mbobs_data_cache[key][sbkey]["bmask"] = bmask 

555 self._mbobs_data_cache[key][sbkey]["ormask"] = ormask 

556 self._mbobs_data_cache[key][sbkey]["psf_stats"] = psf_stats 

557 

558 return self._mbobs_data_cache[key][sbkey] 

559 

560 def _measure( 

561 self, *, mbobs_list, shear_bands, cat, shear_str, mfrac, bmask, 

562 ormask, psf_stats, det_bands, 

563 ): 

564 

565 t0 = time.time() 

566 all_res = [] 

567 for fitter, fwhm_reg, is_wavg, symm, coadd in zip( 

568 self._fitters, self._fwhm_regs, 

569 self._fitter_is_wavg, self._fitter_symmetrize, 

570 self._fitter_coadd, 

571 ): 

572 ft0 = time.time() 

573 if is_wavg: 

574 res = fit_mbobs_list_wavg( 

575 mbobs_list=mbobs_list, 

576 fitter=fitter, 

577 shear_bands=shear_bands, 

578 bmask_flags=self.get("bmask_flags", 0), 

579 fwhm_reg=fwhm_reg, 

580 symmetrize=symm, 

581 ) 

582 else: 

583 res = fit_mbobs_list_joint( 

584 mbobs_list=mbobs_list, 

585 fitter_name=fitter, 

586 shear_bands=shear_bands, 

587 bmask_flags=self.get("bmask_flags", 0), 

588 rng=self.rng, 

589 symmetrize=symm, 

590 coadd=coadd, 

591 ) 

592 ft0 = time.time() - ft0 

593 logger.info( 

594 "fitter %s took %s seconds", 

595 fitter.kind 

596 if hasattr(fitter, "kind") 

597 else fitter, 

598 ft0, 

599 ) 

600 all_res.append(res) 

601 

602 res = combine_fit_res(all_res) 

603 

604 if res is not None: 

605 res = self._add_positions_and_psf( 

606 cat=cat, 

607 res=res, 

608 shear_str=shear_str, 

609 mfrac=mfrac, 

610 bmask=bmask, 

611 ormask=ormask, 

612 psf_stats=psf_stats, 

613 det_bands=det_bands, 

614 ) 

615 logger.info("src measurements took %s seconds", time.time() - t0) 

616 

617 return res 

618 

619 def _add_positions_and_psf( 

620 self, *, cat, res, shear_str, mfrac, bmask, ormask, psf_stats, det_bands, 

621 ): 

622 """ 

623 add catalog positions to the result 

624 """ 

625 

626 new_dt = [ 

627 ('sx_row', 'f4'), 

628 ('sx_col', 'f4'), 

629 ('sx_row_noshear', 'f4'), 

630 ('sx_col_noshear', 'f4'), 

631 ('ormask', 'i4'), 

632 ('mfrac', 'f4'), 

633 ('bmask', 'i4'), 

634 ('mfrac_img', 'f4'), 

635 ('ormask_noshear', 'i4'), 

636 ('mfrac_noshear', 'f4'), 

637 ('bmask_noshear', 'i4'), 

638 ("det_bands", "U%d" % MAX_NUM_SHEAR_BANDS), 

639 ] 

640 if 'psfrec_flags' not in res.dtype.names: 

641 new_dt += [ 

642 ('psfrec_flags', 'i4'), # psfrec is the original psf 

643 ('psfrec_g', 'f8', 2), 

644 ('psfrec_T', 'f8'), 

645 ] 

646 newres = eu.numpy_util.add_fields( 

647 res, 

648 new_dt, 

649 ) 

650 

651 assert len(det_bands) <= MAX_NUM_SHEAR_BANDS 

652 newres["det_bands"] = "".join("%s" % b for b in sorted(det_bands)) 

653 

654 if 'psfrec_flags' not in res.dtype.names: 

655 newres['psfrec_flags'] = procflags.NO_ATTEMPT 

656 

657 newres['psfrec_flags'][:] = psf_stats['flags'] 

658 newres['psfrec_g'][:, 0] = psf_stats['g1'] 

659 newres['psfrec_g'][:, 1] = psf_stats['g2'] 

660 newres['psfrec_T'][:] = psf_stats['T'] 

661 newres['mfrac_img'][:] = np.mean(mfrac) 

662 

663 if cat.size > 0: 

664 obs = self.mbobs[0][0] 

665 

666 newres['sx_col'] = cat['x'] 

667 newres['sx_row'] = cat['y'] 

668 

669 rows_noshear, cols_noshear = shearpos.unshear_positions_obs( 

670 newres['sx_row'], 

671 newres['sx_col'], 

672 shear_str, 

673 obs, # an example for jacobian and image shape 

674 # default is 0.01 but make sure to use the passed in default 

675 # if needed 

676 step=self['metacal'].get("step", shearpos.DEFAULT_STEP), 

677 ) 

678 

679 newres['sx_row_noshear'] = rows_noshear 

680 newres['sx_col_noshear'] = cols_noshear 

681 

682 if 'ormask_region' in self and self['ormask_region'] > 1: 

683 ormask_region = self['ormask_region'] 

684 elif 'mask_region' in self and self['mask_region'] > 1: 

685 ormask_region = self['mask_region'] 

686 else: 

687 ormask_region = 1 

688 

689 if 'mask_region' in self and self['mask_region'] > 1: 

690 bmask_region = self['mask_region'] 

691 else: 

692 bmask_region = 1 

693 

694 logger.debug( 

695 'ormask|bmask region: %s|%s', 

696 ormask_region, 

697 bmask_region, 

698 ) 

699 

700 newres["ormask"] = _fill_in_mask_col( 

701 mask_region=ormask_region, 

702 rows=newres['sx_row'], 

703 cols=newres['sx_col'], 

704 mask=ormask, 

705 ) 

706 newres["ormask_noshear"] = _fill_in_mask_col( 

707 mask_region=ormask_region, 

708 rows=newres['sx_row_noshear'], 

709 cols=newres['sx_col_noshear'], 

710 mask=ormask, 

711 ) 

712 

713 newres["bmask"] = _fill_in_mask_col( 

714 mask_region=bmask_region, 

715 rows=newres['sx_row'], 

716 cols=newres['sx_col'], 

717 mask=bmask, 

718 ) 

719 newres["bmask_noshear"] = _fill_in_mask_col( 

720 mask_region=bmask_region, 

721 rows=newres['sx_row_noshear'], 

722 cols=newres['sx_col_noshear'], 

723 mask=bmask, 

724 ) 

725 

726 if np.any(mfrac > 0): 

727 newres["mfrac"] = measure_mfrac( 

728 mfrac=mfrac, 

729 x=newres["sx_col"], 

730 y=newres["sx_row"], 

731 box_sizes=cat["box_size"], 

732 obs=obs, 

733 fwhm=self.get("mfrac_fwhm", self._fwhms[0]), 

734 ) 

735 

736 newres["mfrac_noshear"] = measure_mfrac( 

737 mfrac=mfrac, 

738 x=newres["sx_col_noshear"], 

739 y=newres["sx_row_noshear"], 

740 box_sizes=cat["box_size"], 

741 obs=obs, 

742 fwhm=self.get("mfrac_fwhm", self._fwhms[0]), 

743 ) 

744 else: 

745 newres["mfrac"] = 0 

746 newres["mfrac_noshear"] = 0 

747 

748 return newres 

749 

750 def _do_detect(self, mbobs, det_bands): 

751 """ 

752 use a MEDSifier to run detection 

753 """ 

754 from . import detect 

755 

756 t0 = time.time() 

757 det_mbobs = ngmix.MultiBandObsList() 

758 for band in det_bands: 

759 det_mbobs.append(mbobs[band]) 

760 

761 medsifier = detect.MEDSifier( 

762 mbobs=det_mbobs, 

763 sx_config=self.get('sx', None), 

764 meds_config=self['meds'], 

765 nodet_flags=self['nodet_flags'], 

766 ) 

767 

768 if self._show: 

769 import descwl_coadd.vis 

770 descwl_coadd.vis.show_image(medsifier.seg) 

771 

772 all_medsifier = detect.CatalogMEDSifier( 

773 mbobs, 

774 medsifier.cat['x'], 

775 medsifier.cat['y'], 

776 medsifier.cat['box_size'], 

777 seg=medsifier.seg, 

778 number=medsifier.cat['number'], 

779 ) 

780 mbm = all_medsifier.get_multiband_meds() 

781 mbobs_list = mbm.get_mbobs_list( 

782 weight_type=self["meds"].get("weight_type", "weight"), 

783 ) 

784 logger.info("detect took %s seconds", time.time() - t0) 

785 

786 return medsifier.cat, mbobs_list 

787 

788 def _get_all_metacal(self, mbobs): 

789 """ 

790 get the sheared versions of the observations 

791 """ 

792 t0 = time.time() 

793 try: 

794 odict = ngmix.metacal.get_all_metacal( 

795 mbobs, 

796 rng=self.rng, 

797 **self['metacal'] 

798 ) 

799 except BootPSFFailure: 

800 odict = None 

801 logger.info("metacal took %s seconds", time.time() - t0) 

802 

803 if self._show and odict is not None: 

804 import descwl_coadd.vis 

805 

806 orig_mbobs = mbobs 

807 

808 for mtype, mbobs in odict.items(): 

809 for band in range(len(mbobs)): 

810 

811 obslist = mbobs[band] 

812 orig_obslist = orig_mbobs[band] 

813 

814 for iobs in range(len(obslist)): 

815 obs = obslist[iobs] 

816 orig_obs = orig_obslist[iobs] 

817 

818 descwl_coadd.vis.show_images( 

819 [ 

820 obs.image, 

821 obs.bmask, 

822 obs.weight, 

823 orig_obs.noise, 

824 ], 

825 title=mtype, 

826 ) 

827 

828 return odict 

829 

830 

831def _get_psf_stats(mbobs, global_flags): 

832 if global_flags != 0: 

833 flags = procflags.PSF_FAILURE | global_flags 

834 g1 = np.nan 

835 g2 = np.nan 

836 T = np.nan 

837 else: 

838 try: 

839 g1sum = 0.0 

840 g2sum = 0.0 

841 Tsum = 0.0 

842 wsum = 0.0 

843 

844 for obslist in mbobs: 

845 for obs in obslist: 

846 wt = obs.weight.max() 

847 res = obs.psf.meta['result'] 

848 T = res['T'] 

849 if 'e' in res: