Coverage for tests/test_deblend_task.py: 99%

1# This file is part of meas_extensions_scarlet. 

2# 

3# Developed for the LSST Data Management System. 

4# This product includes software developed by the LSST Project 

5# (https://www.lsst.org). 

6# See the COPYRIGHT file at the top-level directory of this distribution 

7# for details of code ownership. 

8# 

9# This program is free software: you can redistribute it and/or modify 

10# it under the terms of the GNU General Public License as published by 

11# the Free Software Foundation, either version 3 of the License, or 

12# (at your option) any later version. 

13# 

14# This program is distributed in the hope that it will be useful, 

15# but WITHOUT ANY WARRANTY; without even the implied warranty of 

16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 

17# GNU General Public License for more details. 

18# 

19# You should have received a copy of the GNU General Public License 

20# along with this program. If not, see <https://www.gnu.org/licenses/>. 

21 

22"""Tests for ``ScarletDeblendTask``. 

23 

24Exercises the deblend task on the cached ``multi-blend`` scene from 

25``pipeline.py``: catalog structure, heavy-footprint attachment and 

26model recovery, and skip / failure semantics. The skip tests run 

27against targeted single-blend scenes. 

28""" 

29 

30import unittest 

31from unittest.mock import patch 

32 

33import lsst.afw.detection as afwDet 

34import lsst.afw.image as afwImage 

35import lsst.afw.table as afwTable 

36import lsst.meas.extensions.scarlet as mes 

37import lsst.scarlet.lite as scl 

38import lsst.utils.tests 

39import numpy as np 

40from lsst.afw.detection import PeakTable 

41from lsst.afw.geom import SpanSet 

42from lsst.afw.table import Schema 

43from lsst.geom import Box2I, Point2I 

44from lsst.meas.extensions.scarlet.scarletDeblendTask import ( 

45 ScarletDeblendContext, 

46 ScarletDeblendTask, 

47 deblend, 

48) 

49from lsst.meas.extensions.scarlet.utils import nonzeroBandSupport 

50 

51import pipeline 

52from scenes import SCENES 

53 

54 

55class TestDeblendTask(lsst.utils.tests.TestCase): 

56 """Tests for ``ScarletDeblendTask`` in 

57 ``lsst.meas.extensions.scarlet.scarletDeblendTask``. 

58 

59 Exercises the deblend task's catalog structure, heavy-footprint 

60 attachment / model recovery, and skip / failure semantics. The 

61 ``multi-blend`` scene from ``pipeline.py`` is shared across the 

62 catalog-structure and heavy-footprint tests; the skip tests run 

63 on targeted single-blend scenes. 

64 """ 

65 

66 def _deblend(self, scene, config=None): 

67 # Run the cached pipeline through deblend on a given scene 

68 # (with an optional non-default deblend config). 

69 # ``pipeline.deblend`` memoizes by (scene, configs), so repeated 

70 # calls across tests with the same arguments reuse the bundle. 

71 image = pipeline.build_image(scene) 

72 detection = pipeline.detect(image) 

73 deconv = pipeline.deconvolve(detection) 

74 return pipeline.deblend(deconv, config=config) 

75 

76 def _attach_band_footprints(self, bundle, band, useFlux): 

77 # Hydrate every child row in the bundle's catalog with the 

78 # HeavyFootprint for one (band, useFlux) combination. 

79 image = bundle.image 

80 imageForRedistribution = image.mCoadd[band] if useFlux else None 

81 mes.io.updateCatalogFootprints( 

82 bundle.result.scarletModelData, 

83 bundle.result.deblendedCatalog, 

84 band=band, 

85 imageForRedistribution=imageForRedistribution, 

86 removeScarletData=False, 

87 updateFluxColumns=True, 

88 ) 

89 

90 def _iter_multipeak_children(self, bundle): 

91 # Yield ``(parent, child)`` for every child row whose parent 

92 # is a top-level multi-peak source — i.e. the rows that came 

93 # from the scarlet deblend rather than from an isolated parent. 

94 catalog = bundle.result.deblendedCatalog 

95 objectParents = bundle.result.objectParents 

96 parents = objectParents[ 

97 (objectParents["parent"] == 0) & (objectParents["deblend_nPeaks"] > 1) 

98 ] 

99 for parent in parents: 

100 for child in catalog[catalog["parent"] == parent.get("id")]: 

101 yield parent, child 

102 

103 def _scarlet_blend_for_child(self, bundle, parent, child, band): 

104 # Reconstruct the per-band scarlet blend and pick out the 

105 # source matching ``child.getId()``. Returns the blend (whose 

106 # ``observation`` can be rebound for flux redistribution), the 

107 # source, and the parent's afw footprint. 

108 modelData = bundle.result.scarletModelData 

109 parentBlendData = modelData.blends[parent.getId()] 

110 parentFootprint = parent.getFootprint() 

111 

112 blendData = parentBlendData.children[child["deblend_blendId"]] 

113 full_blend = blendData.minimal_data_to_blend( 

114 model_psf=modelData.model_psf[None, :, :], 

115 psf=modelData.psf, 

116 bands=modelData.bands, 

117 ) 

118 blend = full_blend[band] 

119 source = next( 

120 src for src in blend.sources if src.metadata["id"] == child.getId() 

121 ) 

122 return blend, source, parentFootprint 

123 

124 def _neighbor_overlap(self, blend, source): 

125 # Neighbor flux on ``source``'s support, summed one neighbor at a 

126 # time rather than via ``blend.get_model()`` as in 

127 # ``setDeblenderMetrics``. Returns per-band (max, sum). 

128 model = source.get_model() 

129 support = nonzeroBandSupport(model.data) 

130 neighbors = np.zeros(model.data.shape, dtype=float) 

131 for other in blend.sources: 

132 if other is not source: 

133 neighbors += other.get_model().project(bbox=model.bbox).data 

134 onSupport = neighbors[:, support] 

135 return onSupport.max(axis=1, initial=0.0), onSupport.sum(axis=1) 

136 

137 def test_skip_too_big(self): 

138 """A parent footprint exceeding ``maxFootprintArea`` is skipped 

139 with the ``deblend_skipped_parentTooBig`` flag set. 

140 

141 Uses the ``large_two_sersic`` scene (single parent, two large 

142 overlapping Sersics) with ``maxFootprintArea=2000``; the 

143 deconvolved footprint exceeds that limit so the parent is 

144 skipped. 

145 """ 

146 config = ScarletDeblendTask.ConfigClass() 

147 config.maxFootprintArea = 2000 

148 config.catchFailures = False 

149 

150 bundle = self._deblend(SCENES["large_two_sersic"], config=config) 

151 

152 catalog = bundle.result.objectParents 

153 parents = catalog[catalog["parent"] == 0] 

154 self.assertEqual(len(parents), 1) 

155 parent = parents[0] 

156 self.assertTrue(parent.get("deblend_skipped")) 

157 self.assertTrue(parent.get("deblend_skipped_parentTooBig")) 

158 self.assertFalse(parent.get("deblend_skipped_tooManyPeaks")) 

159 

160 def test_skip_too_many_peaks(self): 

161 """A parent with more than ``maxNumberOfPeaks`` peaks is 

162 skipped with the ``deblend_skipped_tooManyPeaks`` flag set. 

163 

164 Uses the ``three_source_blend`` scene (single parent, three 

165 peaks) with ``maxNumberOfPeaks=2``. 

166 

167 Also pins that a skipped blend leaves 

168 ``deblend_blendConvergenceFailedFlag`` unset: it was never fit, 

169 so convergence does not apply (finding C-1 of the 

170 ``audits/audit-2026-05-05.md`` audit). 

171 """ 

172 config = ScarletDeblendTask.ConfigClass() 

173 config.maxNumberOfPeaks = 2 

174 config.catchFailures = False 

175 

176 bundle = self._deblend(SCENES["three_source_blend"], config=config) 

177 

178 catalog = bundle.result.objectParents 

179 parents = catalog[catalog["parent"] == 0] 

180 self.assertEqual(len(parents), 1) 

181 parent = parents[0] 

182 self.assertTrue(parent.get("deblend_skipped")) 

183 self.assertTrue(parent.get("deblend_skipped_tooManyPeaks")) 

184 self.assertFalse(parent.get("deblend_skipped_parentTooBig")) 

185 self.assertFalse(parent.get("deblend_blendConvergenceFailedFlag")) 

186 

187 def test_all_subblends_failed_updates_parent_record(self): 

188 """When every sub-blend of a multi-peak parent fails, the 

189 aggregate ``deblend_*`` columns must land on the *parent* 

190 record, not on the last sub-blend. 

191 

192 Triggers the all-sub-blends-failed branch by capping 

193 ``maxNumberOfPeaks=2`` on the three-peak ``three_source_blend``, 

194 so every sub-blend is skipped as ``deblend_skipped_tooManyPeaks`` 

195 and the aggregate-update branch runs. Asserts the parent record 

196 carries the summary (``deblend_nPeaks`` matches the parent 

197 footprint's peak count). Under the bug the aggregate landed on 

198 the trailing sub-blend's record and the parent stayed at default. 

199 Regression test for finding C-5 of the 

200 ``audits/audit-2026-05-05.md`` audit. 

201 """ 

202 config = ScarletDeblendTask.ConfigClass() 

203 config.maxNumberOfPeaks = 2 

204 config.catchFailures = False 

205 bundle = self._deblend(SCENES["three_source_blend"], config=config) 

206 catalog = bundle.result.objectParents 

207 parents = catalog[catalog["parent"] == 0] 

208 self.assertEqual(len(parents), 1) 

209 parent = parents[0] 

210 parentNPeaks = len(parent.getFootprint().peaks) 

211 self.assertEqual(parent.get("deblend_nPeaks"), parentNPeaks) 

212 

213 def test_skip_doesnt_affect_other_parents(self): 

214 """One parent being skipped does not affect deblending of 

215 other parents. 

216 

217 The ``multi-blend`` scene detects four parents; with 

218 ``maxNumberOfPeaks=2`` exactly one (the three-peak blend) is 

219 skipped. The remaining three parents are not skipped, and 

220 each non-isolated one produces ``nChild == nPeaks`` children 

221 — the deblend invariant. 

222 """ 

223 config = ScarletDeblendTask.ConfigClass() 

224 config.maxNumberOfPeaks = 2 

225 config.catchFailures = False 

226 

227 bundle = self._deblend(SCENES["multi-blend"], config=config) 

228 

229 catalog = bundle.result.objectParents 

230 parents = catalog[catalog["parent"] == 0] 

231 self.assertEqual(np.sum(parents["deblend_skipped"]), 1) 

232 

233 skipped = parents[parents["deblend_skipped"]] 

234 self.assertTrue(skipped[0].get("deblend_skipped_tooManyPeaks")) 

235 

236 non_skipped = parents[~parents["deblend_skipped"]] 

237 self.assertEqual(len(non_skipped), 3) 

238 blend_parents = non_skipped[~non_skipped["deblend_skipped_isolatedParent"]] 

239 self.assertEqual(len(blend_parents), 2) 

240 for p in blend_parents: 

241 self.assertEqual(p.get("deblend_nChild"), p.get("deblend_nPeaks")) 

242 

243 def test_convergence_flag_false_when_converged(self): 

244 """``deblend_blendConvergenceFailedFlag`` is unset for a blend 

245 that reaches convergence. 

246 

247 The ``three_source_blend`` parent converges in well under 

248 ``maxIter`` iterations under the default config, so the flag — 

249 documented as "at least one source in the blend failed to 

250 converge" — must be `False`. Regression test for finding C-1 

251 of the ``audits/audit-2026-05-05.md`` audit (the flag was 

252 previously stored with inverted semantics, reporting a 

253 converged blend as failed). 

254 """ 

255 defaultMaxIter = ScarletDeblendTask.ConfigClass().maxIter 

256 bundle = self._deblend(SCENES["three_source_blend"]) 

257 parents = bundle.result.objectParents 

258 parents = parents[parents["parent"] == 0] 

259 self.assertEqual(len(parents), 1) 

260 parent = parents[0] 

261 # The blend stopped because it converged, not because it hit 

262 # the iteration cap. 

263 self.assertLess(parent.get("deblend_iterations"), defaultMaxIter) 

264 self.assertFalse(parent.get("deblend_blendConvergenceFailedFlag")) 

265 

266 def test_convergence_flag_true_when_not_converged(self): 

267 """``deblend_blendConvergenceFailedFlag`` is set for a blend 

268 that exhausts ``maxIter`` without converging. 

269 

270 Capping ``maxIter`` at 2 forces the ``three_source_blend`` 

271 parent to stop at the iteration limit before the relative-error 

272 criterion is met, so the flag must be `True`. Regression test 

273 for finding C-1 of the ``audits/audit-2026-05-05.md`` audit. 

274 """ 

275 config = ScarletDeblendTask.ConfigClass() 

276 config.maxIter = 2 

277 config.catchFailures = False 

278 bundle = self._deblend(SCENES["three_source_blend"], config=config) 

279 parents = bundle.result.objectParents 

280 parents = parents[parents["parent"] == 0] 

281 self.assertEqual(len(parents), 1) 

282 parent = parents[0] 

283 # The blend stopped at the iteration cap, i.e. it did not converge. 

284 self.assertEqual(parent.get("deblend_iterations"), config.maxIter) 

285 self.assertTrue(parent.get("deblend_blendConvergenceFailedFlag")) 

286 

287 def test_sub_blend_progress_log_distinguishes_inner_loop(self): 

288 """The periodic logger inside the sub-blend loop names the 

289 sub-blend, not the outer parent count. 

290 

291 Under the bug the inner-loop ``periodicLog.log`` reused the 

292 outer loop's format string, so long sub-blend runs of a 

293 single top-level parent kept emitting "Deblended N out of M 

294 parents" with stale ``N`` — making the task look stalled. 

295 Patches ``PeriodicLogger.LOGGING_INTERVAL`` so every ``log()`` 

296 call fires regardless of wall time, then asserts at least one 

297 captured INFO message identifies a sub-blend. Regression test 

298 for finding DB-2 of the ``audits/audit-2026-05-05.md`` audit. 

299 """ 

300 # Perturb a config field that no other test touches so the 

301 # ``pipeline.deblend`` memoization cache misses and the task 

302 # actually runs inside the ``assertLogs`` block. 

303 config = ScarletDeblendTask.ConfigClass() 

304 config.minIter = 1 

305 config.catchFailures = False 

306 

307 # Negative interval guarantees ``time.time() > next_log_time`` 

308 # on every call. 

309 with patch( 

310 "lsst.utils.logging.PeriodicLogger.LOGGING_INTERVAL", -1.0 

311 ): 

312 with self.assertLogs(level="INFO") as logs: 

313 self._deblend(SCENES["multi-blend"], config=config) 

314 

315 self.assertTrue( 

316 any("sub-blend" in record for record in logs.output), 

317 f"No sub-blend progress message in: {logs.output}", 

318 ) 

319 

320 def test_is_masked_combines_bands_with_and(self): 

321 """``_isMasked`` counts a pixel toward the ``maskLimits`` 

322 fraction only when *every* band has that mask bit set. 

323 

324 This mirrors ``buildObservation``'s per-band weight zeroing 

325 — a pixel only becomes truly unconstrained when masked in 

326 all bands; if even one band leaves it unmasked, that band's 

327 data still constrains it. Sets ``INTRP`` in band 0 only 

328 across a multi-peak parent's footprint and asserts the 

329 parent is not flagged masked. Under the bug (OR across 

330 bands) every footprint pixel had ``INTRP`` set in *some* 

331 band and the parent was skipped. A second pass sets 

332 ``INTRP`` in every band and asserts the parent is then 

333 flagged, confirming AND still triggers when the bit is 

334 universal. Regression test for finding DB-3 of the 

335 ``audits/audit-2026-05-05.md`` audit. 

336 """ 

337 bundle = self._deblend(SCENES["multi-blend"]) 

338 parents = bundle.result.objectParents 

339 parents = parents[parents["parent"] == 0] 

340 parent = next(p for p in parents if p["deblend_nPeaks"] > 1) 

341 footprint = parent.getFootprint() 

342 

343 # Deep-copy the per-band exposures so the cached bundle's 

344 # mask plane isn't corrupted for downstream tests. 

345 bands = bundle.image.mCoadd.bands 

346 mExposure = afwImage.MultibandExposure.fromExposures( 

347 bands, [exp.clone() for exp in bundle.image.mCoadd] 

348 ) 

349 

350 config = ScarletDeblendTask.ConfigClass() 

351 config.maskLimits = {"INTRP": 0.05} 

352 task = ScarletDeblendTask( 

353 schema=Schema(bundle.deconvolved.detection.schema), 

354 config=config, 

355 ) 

356 

357 intrp = mExposure.mask.getPlaneBitMask("INTRP") 

358 

359 # Sanity: a clean mask plane is not flagged. 

360 self.assertFalse(task._isMasked(footprint, mExposure)) 

361 

362 # INTRP in band 0 only — not flagged under AND, was under OR. 

363 mExposure.mask.array[0] |= intrp 

364 self.assertFalse(task._isMasked(footprint, mExposure)) 

365 

366 # INTRP in every band — flagged under AND as well. 

367 for b in range(len(bands)): 

368 mExposure.mask.array[b] |= intrp 

369 self.assertTrue(task._isMasked(footprint, mExposure)) 

370 

371 def test_init_warns_on_unknown_pseudo_columns(self): 

372 """``ScarletDeblendTask.__init__`` warns when a configured 

373 ``pseudoColumn`` is not present on the deblender schema. 

374 

375 ``isPseudoSource`` deliberately swallows ``KeyError`` so the 

376 same ``pseudoColumns`` list can be evaluated against 

377 ``PeakRecord`` and ``SourceRecord``. The schema-presence check 

378 runs once at task construction: any pseudoColumn missing from 

379 the input source schema (and from ``peakSchema`` if provided) 

380 is reported as a likely config typo or missing-upstream-task 

381 signal, since the science pipeline always adds the relevant 

382 columns before the deblender runs. Regression test for 

383 finding DB-17 of the ``audits/audit-2026-05-05.md`` audit. 

384 """ 

385 schema = afwTable.SourceTable.makeMinimalSchema() 

386 config = ScarletDeblendTask.ConfigClass() 

387 config.pseudoColumns = ["definitely_not_a_real_column"] 

388 

389 with self.assertLogs(level="WARNING") as logs: 

390 ScarletDeblendTask(schema=schema, config=config) 

391 

392 self.assertTrue( 

393 any( 

394 "definitely_not_a_real_column" in record 

395 for record in logs.output 

396 ), 

397 f"No unknown-pseudoColumn warning in: {logs.output}", 

398 ) 

399 

400 def test_max_iter_zero_skips_fit_cleanly(self): 

401 """With ``maxIter=0`` the optimizer never runs, so the parent 

402 record's ``deblend_iterations`` must report 0 (not the synthetic 

403 2 the old code produced by stuffing two equal entries into 

404 ``blend.loss`` to keep downstream consumers indexable). 

405 

406 Also pins ``deblend_blendConvergenceFailedFlag`` False: a blend 

407 that was never fit hasn't failed convergence, and 

408 ``_checkBlendConvergence`` must handle the empty-loss case 

409 without crashing. Regression test for finding DB-1 of the 

410 ``audits/audit-2026-05-05.md`` audit. 

411 """ 

412 config = ScarletDeblendTask.ConfigClass() 

413 config.maxIter = 0 

414 config.catchFailures = False 

415 

416 bundle = self._deblend(SCENES["three_source_blend"], config=config) 

417 parents = bundle.result.objectParents 

418 parents = parents[parents["parent"] == 0] 

419 self.assertEqual(len(parents), 1) 

420 parent = parents[0] 

421 

422 self.assertEqual(parent.get("deblend_iterations"), 0) 

423 self.assertFalse(parent.get("deblend_blendConvergenceFailedFlag")) 

424 

425 def test_detected_peak_skips_pseudo_peaks(self): 

426 """Each deblended source's ``detectedPeak`` is the real peak at 

427 its own center, even when a pseudo peak precedes it in the 

428 footprint's peak list. 

429 

430 ``deblend`` filters pseudo peaks (e.g. sky objects) out of the 

431 list it initializes sources from, so the back-pointer to the 

432 ``PeakRecord`` must be indexed in that *filtered* list. Indexing 

433 the unfiltered ``footprint.peaks`` shifts every source's 

434 ``detectedPeak`` by the number of preceding pseudo peaks. 

435 Regression test for finding C-4 of the 

436 ``audits/audit-2026-05-05.md`` audit. 

437 """ 

438 deconv = pipeline.deconvolve( 

439 pipeline.detect(pipeline.build_image(SCENES["three_source_blend"])) 

440 ) 

441 config = ScarletDeblendTask.ConfigClass() 

442 context = ScarletDeblendContext.build( 

443 deconv.image.mCoadd, deconv.mDeconvolved, deconv.detection.catalog, config 

444 ) 

445 

446 # Rebuild the (single, three-peak) parent footprint with a sky 

447 # pseudo peak prepended ahead of the three real peaks. 

448 origFootprint = deconv.detection.catalog[0].getFootprint() 

449 peakSchema = PeakTable.makeMinimalSchema() 

450 skyKey = peakSchema.addField( 

451 "merge_peak_sky", type="Flag", doc="sky pseudo peak" 

452 ) 

453 footprint = afwDet.Footprint(origFootprint.spans, peakSchema) 

454 pseudoPeak = footprint.addPeak(12, 14, 1.0) 

455 pseudoPeak.set(skyKey, True) 

456 pseudoId = pseudoPeak.getId() 

457 for peak in origFootprint.peaks: 

458 footprint.addPeak(peak.getIx(), peak.getIy(), 10.0) 

459 

460 blend = deblend(context, footprint, config, spectrumInit=False) 

461 

462 self.assertEqual(len(blend.sources), len(origFootprint.peaks)) 

463 for source in blend.sources: 

464 detectedPeak = source.detectedPeak 

465 # The real peaks sit exactly on the source centers, so a 

466 # correctly-mapped detectedPeak lands on its source's center 

467 # (center is ordered (y, x)). 

468 self.assertEqual(detectedPeak.getIy(), source.center[0]) 

469 self.assertEqual(detectedPeak.getIx(), source.center[1]) 

470 self.assertNotEqual(detectedPeak.getId(), pseudoId) 

471 

472 def test_parent_peak_mapper_propagates_extra_peak_fields(self): 

473 """``ScarletDeblendTask.parentPeakSchemaMapper`` carries extra 

474 ``merge_peak_*`` peak fields onto a parent record. 

475 

476 The task copies a parent footprint's first peak onto the 

477 parent source record via this mapper. Without the extra-field 

478 mappings any ``merge_peak_*`` flags (for example 

479 ``merge_peak_sky``) on that peak are silently dropped, 

480 breaking propagation of pseudo-source flags to deconvolved 

481 sub-blend parents. Regression test for finding C-8 of the 

482 ``audits/audit-2026-05-05.md`` audit. 

483 """ 

484 schema = afwTable.SourceTable.makeMinimalSchema() 

485 peakSchema = PeakTable.makeMinimalSchema() 

486 skyKey = peakSchema.addField( 

487 "merge_peak_sky", type="Flag", doc="sky pseudo peak" 

488 ) 

489 task = ScarletDeblendTask(schema=schema, peakSchema=peakSchema) 

490 

491 peakCat = afwDet.PeakCatalog(afwDet.PeakTable.make(peakSchema)) 

492 peak = peakCat.addNew() 

493 peak.set(skyKey, True) 

494 

495 parentCatalog = afwTable.SourceCatalog( 

496 afwTable.SourceTable.make(task.parentSchema) 

497 ) 

498 parent = parentCatalog.addNew() 

499 parent.assign(peak, task.parentPeakSchemaMapper) 

500 

501 self.assertTrue(parent.get("merge_peak_sky")) 

502 

503 def test_build_intersecting_footprints_rejects_out_of_bounds_peak(self): 

504 """A peak whose pixel coordinates fall outside the bbox of 

505 ``footprintImage`` raises ``RuntimeError``, in every direction. 

506 

507 ``_buildIntersectingFootprints`` indexed ``footprintImage.data`` 

508 with the peak's offset from the image bbox origin inside a 

509 ``try/except IndexError``. NumPy only raises ``IndexError`` for 

510 out-of-range *positive* indices; negative indices silently wrap 

511 from the opposite edge of the array, so a peak west or south 

512 of the bbox origin used to be silently mapped to an unrelated 

513 pixel on the opposite edge and processed as if it lay there. 

514 The bounds check must reject all four directions. Regression 

515 test for finding C-9 of the ``audits/audit-2026-05-05.md`` 

516 audit. 

517 """ 

518 schema = afwTable.SourceTable.makeMinimalSchema() 

519 task = ScarletDeblendTask(schema=schema) 

520 # 5x5 footprintImage filled with zeros at origin (10, 10), so 

521 # NumPy wraparound from a negative index lands on a 0 cell and 

522 # the inner loop would silently skip the peak under the bug. 

523 footprintImage = scl.Image( 

524 np.zeros((5, 5), dtype=np.int32), yx0=(10, 10) 

525 ) 

526 parentCatalog = afwTable.SourceCatalog( 

527 afwTable.SourceTable.make(task.parentSchema) 

528 ) 

529 

530 # One peak per out-of-bounds direction relative to the 

531 # ``[10..14, 10..14]`` bbox. ``south`` and ``west`` exercise 

532 # the negative-wraparound branches the bug missed; ``north`` 

533 # and ``east`` exercise the positive-overflow branches that 

534 # NumPy raised on. 

535 directions = { 

536 "west": (5, 12), 

537 "south": (12, 5), 

538 "east": (20, 12), 

539 "north": (12, 20), 

540 } 

541 for name, (x, y) in directions.items(): 

542 with self.subTest(direction=name): 

543 footprint = afwDet.Footprint( 

544 SpanSet(), PeakTable.makeMinimalSchema() 

545 ) 

546 footprint.addPeak(x, y, 1.0) 

547 with self.assertRaises(RuntimeError): 

548 task._buildIntersectingFootprints( 

549 parentId=0, 

550 afwFootprint=footprint, 

551 parentCatalog=parentCatalog, 

552 sclFootprints=[], 

553 footprintImage=footprintImage, 

554 ) 

555 

556 def test_addDeblendedSource_chi2_masks_to_source_footprint(self): 

557 """``_addDeblendedSource`` sums ``chi2`` over only the source's 

558 own positive-model pixels before dividing by the source area. 

559 

560 The ``chi2`` image passed in is the *blend's* chi2, masked at 

561 construction by ``blendModel.data > 0`` — the union of every 

562 source's positive-model footprint. Within one source's 

563 ``scarletSource.bbox`` that mask therefore leaks contributions 

564 from neighbors whose models extend into the same bbox. 

565 Summing the full bbox and dividing by just this source's 

566 positive-pixel count conflates those neighbor residuals into 

567 this source's reduced chi2. The correct calculation masks 

568 ``chi2`` by ``scarletSource.get_model().data > 0`` before 

569 summing so the numerator and denominator are over the same 

570 pixel set. 

571 

572 Builds a synthetic two-band setup with a 4x4 source bbox: the 

573 source's model is positive only in the top-left 2x2 quadrant. 

574 The chi2 image has a small contribution in that quadrant (the 

575 source's own residual) and a much larger contribution in the 

576 disjoint bottom-right 2x2 quadrant (a hypothetical neighbor's 

577 residual leaking into this bbox). The masked reduced chi2 is 

578 1.0; the unmasked formula yields 11.0. 

579 """ 

580 bands = ("g", "r") 

581 nBands = len(bands) 

582 

583 # Source model: positive only in the top-left 2x2 quadrant. 

584 modelData = np.zeros((nBands, 4, 4), dtype=np.float32) 

585 modelData[:, :2, :2] = 1.0 

586 sourceModel = scl.Image(modelData, yx0=(0, 0), bands=bands) 

587 component = scl.component.CubeComponent(model=sourceModel, peak=(0, 0)) 

588 scarletSource = scl.Source([component]) 

589 

590 # chi2: 1.0 per pixel in the source's own quadrant (sum = 8 

591 # across both bands), 10.0 per pixel in the disjoint 

592 # bottom-right "neighbor" quadrant (sum = 80 across both 

593 # bands). area = 2 bands * 4 positive pixels = 8. 

594 # Masked reduced chi2 = 8/8 = 1.0; unmasked = 88/8 = 11.0. 

595 chi2Data = np.zeros((nBands, 4, 4), dtype=np.float32) 

596 chi2Data[:, :2, :2] = 1.0 

597 chi2Data[:, 2:, 2:] = 10.0 

598 chi2 = scl.Image(chi2Data, yx0=(0, 0), bands=bands) 

599 

600 schema = afwTable.SourceTable.makeMinimalSchema() 

601 task = ScarletDeblendTask(schema=schema) 

602 

603 parentCatalog = afwTable.SourceCatalog( 

604 afwTable.SourceTable.make(task.parentSchema) 

605 ) 

606 blendRecord = parentCatalog.addNew() 

607 footprint = afwDet.Footprint( 

608 SpanSet(Box2I(Point2I(0, 0), Point2I(3, 3))), 

609 PeakTable.makeMinimalSchema(), 

610 ) 

611 peak = footprint.addPeak(0, 0, 1.0) 

612 blendRecord.setFootprint(footprint) 

613 

614 objectCatalog = afwTable.SourceCatalog( 

615 afwTable.SourceTable.make(task.objectSchema) 

616 ) 

617 

618 task._addDeblendedSource( 

619 parentId=0, 

620 blendRecord=blendRecord, 

621 peak=peak, 

622 objectCatalog=objectCatalog, 

623 scarletSource=scarletSource, 

624 chi2=chi2, 

625 ) 

626 

627 src = objectCatalog[0] 

628 self.assertAlmostEqual(src.get("deblend_chi2"), 1.0, places=6) 

629 

630 def test_catalog_total_count(self): 

631 """The deblended catalog has one row per input model.""" 

632 bundle = self._deblend(SCENES["multi-blend"]) 

633 catalog = bundle.result.deblendedCatalog 

634 self.assertEqual(len(catalog), len(SCENES["multi-blend"].models)) 

635 

636 def test_catalog_sorted_by_parent_id(self): 

637 """Catalog rows are sorted by parent id.""" 

638 bundle = self._deblend(SCENES["multi-blend"]) 

639 catalog = bundle.result.deblendedCatalog 

640 np.testing.assert_array_equal(sorted(catalog["parent"]), catalog["parent"]) 

641 

642 def test_child_ids_above_parent_ids(self): 

643 """Every child source's id exceeds the maximum parent id.""" 

644 bundle = self._deblend(SCENES["multi-blend"]) 

645 catalog = bundle.result.deblendedCatalog 

646 objectParents = bundle.result.objectParents 

647 maxId = np.max(objectParents["id"]) 

648 for src in catalog: 

649 if src["parent"] > 0: 

650 self.assertGreater(src["id"], maxId) 

651 

652 def test_every_source_has_one_peak(self): 

653 """After ``updateCatalogFootprints``, every catalog source's 

654 footprint has exactly one peak. 

655 

656 Isolated parents keep their (single-peak) detection footprint; 

657 children get HeavyFootprints attached by 

658 ``updateCatalogFootprints`` whose peaks come from the 

659 scarlet model. The deblend invariant is "one peak per 

660 source row" after the catalog has been hydrated. 

661 """ 

662 bundle = self._deblend(SCENES["multi-blend"]) 

663 catalog = bundle.result.deblendedCatalog 

664 modelData = bundle.result.scarletModelData 

665 mes.io.updateCatalogFootprints( 

666 modelData, 

667 catalog, 

668 band=bundle.image.bands[0], 

669 imageForRedistribution=None, 

670 removeScarletData=False, 

671 updateFluxColumns=True, 

672 ) 

673 

674 for src in catalog: 

675 self.assertEqual(len(src.getFootprint().peaks), 1) 

676 

677 def test_nChild_consistency(self): 

678 """``deblend_nChild`` over multi-peak parents accounts for every child. 

679 

680 Sum of ``deblend_nChild`` over the top-level multi-peak parents 

681 equals (total catalog size) − (number of top-level isolated 

682 parents) — i.e. every non-parent row in the catalog is a child 

683 of exactly one multi-peak parent. 

684 """ 

685 bundle = self._deblend(SCENES["multi-blend"]) 

686 catalog = bundle.result.deblendedCatalog 

687 objectParents = bundle.result.objectParents 

688 isolated = catalog[catalog["parent"] == 0] 

689 parents = objectParents[ 

690 (objectParents["parent"] == 0) & (objectParents["deblend_nPeaks"] > 1) 

691 ] 

692 self.assertEqual( 

693 np.sum(parents["deblend_nChild"]), len(catalog) - len(isolated) 

694 ) 

695 

696 def test_isolated_parents_marked(self): 

697 """``deblend_skipped_isolatedParent`` is set on every single-peak 

698 parent. 

699 """ 

700 bundle = self._deblend(SCENES["multi-blend"]) 

701 objectParents = bundle.result.objectParents 

702 isolatedParents = objectParents[ 

703 (objectParents["parent"] == 0) 

704 & (objectParents["deblend_nPeaks"] == 1) 

705 ] 

706 self.assertEqual( 

707 np.sum(objectParents["deblend_skipped_isolatedParent"]), 

708 len(isolatedParents), 

709 ) 

710 

711 def test_isolated_source_persisted(self): 

712 """Each isolated parent has a matching ``SourceData`` in 

713 ``modelData.isolated``: same span array, same origin. 

714 """ 

715 bundle = self._deblend(SCENES["multi-blend"]) 

716 catalog = bundle.result.deblendedCatalog 

717 modelData = bundle.result.scarletModelData 

718 isolated = catalog[catalog["parent"] == 0] 

719 self.assertEqual(len(isolated), len(modelData.isolated)) 

720 for sid, source in modelData.isolated.items(): 

721 catalog_footprint = catalog.find(sid).getFootprint() 

722 np.testing.assert_array_equal( 

723 source.span_array, catalog_footprint.spans.asArray() 

724 ) 

725 self.assertTupleEqual( 

726 source.origin[::-1], tuple(catalog_footprint.getBBox().getMin()) 

727 ) 

728 

729 def test_heavy_footprint_flux_at_peak(self): 

730 """Heavy-footprint flux at the source's peak equals 

731 ``deblend_peak_instFlux`` in every band, with and without 

732 flux redistribution. 

733 

734 This works in the test image because the detected peak is in 

735 the same location as the scarlet peak; if the peak were 

736 shifted the flux value would still be correct but 

737 ``deblend_peak_center`` would not be the right pixel. 

738 """ 

739 bundle = self._deblend(SCENES["multi-blend"]) 

740 for useFlux in [False, True]: 

741 for band in bundle.image.bands: 

742 with self.subTest(band=band, useFlux=useFlux): 

743 self._attach_band_footprints(bundle, band, useFlux) 

744 for _, child in self._iter_multipeak_children(bundle): 

745 fp = child.getFootprint() 

746 img = fp.extractImage(fill=0.0) 

747 px = child.get("deblend_peak_center_x") 

748 py = child.get("deblend_peak_center_y") 

749 flux = img[Point2I(px, py)] 

750 self.assertEqual( 

751 flux, child.get("deblend_peak_instFlux") 

752 ) 

753 

754 def test_heavy_footprint_band_columns_populated(self): 

755 """``updateCatalogFootprints`` populates every band-dependent 

756 ``deblend_*`` column it owns on each deblended child. 

757 

758 Pins both the four pre-existing writes 

759 (``deblend_zeroFlux``, ``deblend_dataCoverage``, 

760 ``deblend_scarletFlux``, ``deblend_peak_instFlux``) and the 

761 four blendedness/overlap metric writes 

762 (``deblend_maxOverlap``, ``deblend_fluxOverlap``, 

763 ``deblend_fluxOverlapFraction``, ``deblend_blendedness``). 

764 The metric writes were previously missing: their values were 

765 computed onto ``source.metrics`` by ``setDeblenderMetrics`` 

766 but never pushed onto the catalog record, so every child 

767 carried the schema's default ``NaN`` for the four 

768 ``np.float32`` metric fields. Runs with ``useFlux=True`` so 

769 the ``deblend_dataCoverage`` branch is also covered. 

770 Regression test for finding DB-7 of the 

771 ``audits/audit-2026-05-05.md`` audit. 

772 

773 The overlap columns are compared against neighbor flux 

774 recomputed from the scarlet sources, since a child in a 

775 multi-peak parent need not overlap any neighbor in a given band 

776 (its neighbors' models can be zero on its support). 

777 """ 

778 bundle = self._deblend(SCENES["multi-blend"]) 

779 band = bundle.image.bands[0] 

780 self._attach_band_footprints(bundle, band, useFlux=True) 

781 

782 nOverlapping = 0 

783 for parent, child in self._iter_multipeak_children(bundle): 

784 with self.subTest(childId=child.getId()): 

785 self.assertFalse(child.get("deblend_zeroFlux")) 

786 self.assertGreater(child.get("deblend_dataCoverage"), 0.0) 

787 self.assertGreater(child.get("deblend_scarletFlux"), 0.0) 

788 self.assertFalse( 

789 np.isnan(child.get("deblend_peak_instFlux")) 

790 ) 

791 blend, source, _ = self._scarlet_blend_for_child( 

792 bundle, parent, child, band 

793 ) 

794 maxOverlap, fluxOverlap = self._neighbor_overlap(blend, source) 

795 # Tolerances allow for the float32 catalog columns. 

796 np.testing.assert_allclose( 

797 child.get("deblend_maxOverlap"), maxOverlap[0], 

798 rtol=1e-5, atol=1e-7, 

799 ) 

800 np.testing.assert_allclose( 

801 child.get("deblend_fluxOverlap"), fluxOverlap[0], 

802 rtol=1e-5, atol=1e-6, 

803 ) 

804 self.assertFalse( 

805 np.isnan(child.get("deblend_fluxOverlapFraction")) 

806 ) 

807 self.assertFalse(np.isnan(child.get("deblend_blendedness"))) 

808 nOverlapping += fluxOverlap[0] > 0 

809 # Guard against a scene in which no child overlaps a neighbor, 

810 # which would make the comparisons above trivially zero. 

811 self.assertGreater(nOverlapping, 0) 

812 

813 def test_heavy_footprint_peak_position(self): 

814 """The HeavyFootprint's peak position and the scarlet model's 

815 source center both match ``deblend_peak_center_{x,y}``. 

816 

817 The assertion is geometric (independent of band / useFlux), 

818 so this runs once against the first band's HeavyFootprints. 

819 """ 

820 bundle = self._deblend(SCENES["multi-blend"]) 

821 band = bundle.image.bands[0] 

822 self._attach_band_footprints(bundle, band, useFlux=False) 

823 

824 for parent, child in self._iter_multipeak_children(bundle): 

825 fp = child.getFootprint() 

826 px = child.get("deblend_peak_center_x") 

827 py = child.get("deblend_peak_center_y") 

828 self.assertEqual(px, fp.getPeaks()[0].getIx())