Coverage for tests/test_deblend_task.py: 13%

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) 

49 

50import pipeline 

51from scenes import SCENES 

52 

53 

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

55 """Tests for ``ScarletDeblendTask`` in 

56 ``lsst.meas.extensions.scarlet.scarletDeblendTask``. 

57 

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

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

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

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

62 on targeted single-blend scenes. 

63 """ 

64 

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

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

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

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

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

70 image = pipeline.build_image(scene) 

71 detection = pipeline.detect(image) 

72 deconv = pipeline.deconvolve(detection) 

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

74 

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

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

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

78 image = bundle.image 

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

80 mes.io.updateCatalogFootprints( 

81 bundle.result.scarletModelData, 

82 bundle.result.deblendedCatalog, 

83 band=band, 

84 imageForRedistribution=imageForRedistribution, 

85 removeScarletData=False, 

86 updateFluxColumns=True, 

87 ) 

88 

89 def _iter_multipeak_children(self, bundle): 

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

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

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

93 catalog = bundle.result.deblendedCatalog 

94 objectParents = bundle.result.objectParents 

95 parents = objectParents[ 

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

97 ] 

98 for parent in parents: 

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

100 yield parent, child 

101 

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

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

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

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

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

107 modelData = bundle.result.scarletModelData 

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

109 parentFootprint = parent.getFootprint() 

110 

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

112 full_blend = blendData.minimal_data_to_blend( 

113 model_psf=modelData.metadata["model_psf"][None, :, :], 

114 psf=modelData.metadata["psf"], 

115 bands=modelData.metadata["bands"], 

116 ) 

117 blend = full_blend[band] 

118 source = next( 

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

120 ) 

121 return blend, source, parentFootprint 

122 

123 def test_skip_too_big(self): 

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

125 with the ``deblend_skipped_parentTooBig`` flag set. 

126 

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

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

129 deconvolved footprint exceeds that limit so the parent is 

130 skipped. 

131 """ 

132 config = ScarletDeblendTask.ConfigClass() 

133 config.maxFootprintArea = 2000 

134 config.catchFailures = False 

135 

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

137 

138 catalog = bundle.result.objectParents 

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

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

141 parent = parents[0] 

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

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

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

145 

146 def test_skip_too_many_peaks(self): 

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

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

149 

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

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

152 

153 Also pins that a skipped blend leaves 

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

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

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

157 """ 

158 config = ScarletDeblendTask.ConfigClass() 

159 config.maxNumberOfPeaks = 2 

160 config.catchFailures = False 

161 

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

163 

164 catalog = bundle.result.objectParents 

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

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

167 parent = parents[0] 

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

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

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

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

172 

173 def test_all_subblends_failed_updates_parent_record(self): 

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

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

176 record, not on the last sub-blend. 

177 

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

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

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

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

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

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

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

185 Regression test for finding C-5 of the 

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

187 """ 

188 config = ScarletDeblendTask.ConfigClass() 

189 config.maxNumberOfPeaks = 2 

190 config.catchFailures = False 

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

192 catalog = bundle.result.objectParents 

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

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

195 parent = parents[0] 

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

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

198 

199 def test_skip_doesnt_affect_other_parents(self): 

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

201 other parents. 

202 

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

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

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

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

207 — the deblend invariant. 

208 """ 

209 config = ScarletDeblendTask.ConfigClass() 

210 config.maxNumberOfPeaks = 2 

211 config.catchFailures = False 

212 

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

214 

215 catalog = bundle.result.objectParents 

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

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

218 

219 skipped = parents[parents["deblend_skipped"]] 

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

221 

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

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

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

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

226 for p in blend_parents: 

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

228 

229 def test_convergence_flag_false_when_converged(self): 

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

231 that reaches convergence. 

232 

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

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

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

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

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

238 previously stored with inverted semantics, reporting a 

239 converged blend as failed). 

240 """ 

241 defaultMaxIter = ScarletDeblendTask.ConfigClass().maxIter 

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

243 parents = bundle.result.objectParents 

244 parents = parents[parents["parent"] == 0] 

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

246 parent = parents[0] 

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

248 # the iteration cap. 

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

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

251 

252 def test_convergence_flag_true_when_not_converged(self): 

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

254 that exhausts ``maxIter`` without converging. 

255 

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

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

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

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

260 """ 

261 config = ScarletDeblendTask.ConfigClass() 

262 config.maxIter = 2 

263 config.catchFailures = False 

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

265 parents = bundle.result.objectParents 

266 parents = parents[parents["parent"] == 0] 

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

268 parent = parents[0] 

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

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

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

272 

273 def test_sub_blend_progress_log_distinguishes_inner_loop(self): 

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

275 sub-blend, not the outer parent count. 

276 

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

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

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

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

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

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

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

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

285 """ 

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

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

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

289 config = ScarletDeblendTask.ConfigClass() 

290 config.minIter = 1 

291 config.catchFailures = False 

292 

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

294 # on every call. 

295 with patch( 

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

297 ): 

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

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

300 

301 self.assertTrue( 

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

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

304 ) 

305 

306 def test_is_masked_combines_bands_with_and(self): 

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

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

309 

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

311 — a pixel only becomes truly unconstrained when masked in 

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

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

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

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

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

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

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

319 flagged, confirming AND still triggers when the bit is 

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

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

322 """ 

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

324 parents = bundle.result.objectParents 

325 parents = parents[parents["parent"] == 0] 

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

327 footprint = parent.getFootprint() 

328 

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

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

331 bands = bundle.image.mCoadd.bands 

332 mExposure = afwImage.MultibandExposure.fromExposures( 

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

334 ) 

335 

336 config = ScarletDeblendTask.ConfigClass() 

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

338 task = ScarletDeblendTask( 

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

340 config=config, 

341 ) 

342 

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

344 

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

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

347 

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

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

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

351 

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

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

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

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

356 

357 def test_init_warns_on_unknown_pseudo_columns(self): 

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

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

360 

361 ``isPseudoSource`` deliberately swallows ``KeyError`` so the 

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

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

364 runs once at task construction: any pseudoColumn missing from 

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

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

367 signal, since the science pipeline always adds the relevant 

368 columns before the deblender runs. Regression test for 

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

370 """ 

371 schema = afwTable.SourceTable.makeMinimalSchema() 

372 config = ScarletDeblendTask.ConfigClass() 

373 config.pseudoColumns = ["definitely_not_a_real_column"] 

374 

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

376 ScarletDeblendTask(schema=schema, config=config) 

377 

378 self.assertTrue( 

379 any( 

380 "definitely_not_a_real_column" in record 

381 for record in logs.output 

382 ), 

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

384 ) 

385 

386 def test_max_iter_zero_skips_fit_cleanly(self): 

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

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

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

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

391 

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

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

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

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

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

397 """ 

398 config = ScarletDeblendTask.ConfigClass() 

399 config.maxIter = 0 

400 config.catchFailures = False 

401 

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

403 parents = bundle.result.objectParents 

404 parents = parents[parents["parent"] == 0] 

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

406 parent = parents[0] 

407 

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

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

410 

411 def test_detected_peak_skips_pseudo_peaks(self): 

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

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

414 footprint's peak list. 

415 

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

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

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

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

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

421 Regression test for finding C-4 of the 

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

423 """ 

424 deconv = pipeline.deconvolve( 

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

426 ) 

427 config = ScarletDeblendTask.ConfigClass() 

428 context = ScarletDeblendContext.build( 

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

430 ) 

431 

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

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

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

435 peakSchema = PeakTable.makeMinimalSchema() 

436 skyKey = peakSchema.addField( 

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

438 ) 

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

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

441 pseudoPeak.set(skyKey, True) 

442 pseudoId = pseudoPeak.getId() 

443 for peak in origFootprint.peaks: 

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

445 

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

447 

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

449 for source in blend.sources: 

450 detectedPeak = source.detectedPeak 

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

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

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

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

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

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

457 

458 def test_parent_peak_mapper_propagates_extra_peak_fields(self): 

459 """``ScarletDeblendTask.parentPeakSchemaMapper`` carries extra 

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

461 

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

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

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

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

466 breaking propagation of pseudo-source flags to deconvolved 

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

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

469 """ 

470 schema = afwTable.SourceTable.makeMinimalSchema() 

471 peakSchema = PeakTable.makeMinimalSchema() 

472 skyKey = peakSchema.addField( 

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

474 ) 

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

476 

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

478 peak = peakCat.addNew() 

479 peak.set(skyKey, True) 

480 

481 parentCatalog = afwTable.SourceCatalog( 

482 afwTable.SourceTable.make(task.parentSchema) 

483 ) 

484 parent = parentCatalog.addNew() 

485 parent.assign(peak, task.parentPeakSchemaMapper) 

486 

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

488 

489 def test_build_intersecting_footprints_rejects_out_of_bounds_peak(self): 

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

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

492 

493 ``_buildIntersectingFootprints`` indexed ``footprintImage.data`` 

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

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

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

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

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

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

500 The bounds check must reject all four directions. Regression 

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

502 audit. 

503 """ 

504 schema = afwTable.SourceTable.makeMinimalSchema() 

505 task = ScarletDeblendTask(schema=schema) 

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

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

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

509 footprintImage = scl.Image( 

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

511 ) 

512 parentCatalog = afwTable.SourceCatalog( 

513 afwTable.SourceTable.make(task.parentSchema) 

514 ) 

515 

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

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

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

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

520 # NumPy raised on. 

521 directions = { 

522 "west": (5, 12), 

523 "south": (12, 5), 

524 "east": (20, 12), 

525 "north": (12, 20), 

526 } 

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

528 with self.subTest(direction=name): 

529 footprint = afwDet.Footprint( 

530 SpanSet(), PeakTable.makeMinimalSchema() 

531 ) 

532 footprint.addPeak(x, y, 1.0) 

533 with self.assertRaises(RuntimeError): 

534 task._buildIntersectingFootprints( 

535 parentId=0, 

536 afwFootprint=footprint, 

537 parentCatalog=parentCatalog, 

538 sclFootprints=[], 

539 footprintImage=footprintImage, 

540 ) 

541 

542 def test_addDeblendedSource_chi2_masks_to_source_footprint(self): 

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

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

545 

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

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

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

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

550 from neighbors whose models extend into the same bbox. 

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

552 positive-pixel count conflates those neighbor residuals into 

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

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

555 summing so the numerator and denominator are over the same 

556 pixel set. 

557 

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

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

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

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

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

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

564 1.0; the unmasked formula yields 11.0. 

565 """ 

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

567 nBands = len(bands) 

568 

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

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

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

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

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

574 scarletSource = scl.Source([component]) 

575 

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

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

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

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

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

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

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

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

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

585 

586 schema = afwTable.SourceTable.makeMinimalSchema() 

587 task = ScarletDeblendTask(schema=schema) 

588 

589 parentCatalog = afwTable.SourceCatalog( 

590 afwTable.SourceTable.make(task.parentSchema) 

591 ) 

592 blendRecord = parentCatalog.addNew() 

593 footprint = afwDet.Footprint( 

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

595 PeakTable.makeMinimalSchema(), 

596 ) 

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

598 blendRecord.setFootprint(footprint) 

599 

600 objectCatalog = afwTable.SourceCatalog( 

601 afwTable.SourceTable.make(task.objectSchema) 

602 ) 

603 

604 task._addDeblendedSource( 

605 parentId=0, 

606 blendRecord=blendRecord, 

607 peak=peak, 

608 objectCatalog=objectCatalog, 

609 scarletSource=scarletSource, 

610 chi2=chi2, 

611 ) 

612 

613 src = objectCatalog[0] 

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

615 

616 def test_catalog_total_count(self): 

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

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

619 catalog = bundle.result.deblendedCatalog 

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

621 

622 def test_catalog_sorted_by_parent_id(self): 

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

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

625 catalog = bundle.result.deblendedCatalog 

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

627 

628 def test_child_ids_above_parent_ids(self): 

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

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

631 catalog = bundle.result.deblendedCatalog 

632 objectParents = bundle.result.objectParents 

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

634 for src in catalog: 

635 if src["parent"] > 0: 

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

637 

638 def test_every_source_has_one_peak(self): 

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

640 footprint has exactly one peak. 

641 

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

643 children get HeavyFootprints attached by 

644 ``updateCatalogFootprints`` whose peaks come from the 

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

646 source row" after the catalog has been hydrated. 

647 """ 

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

649 catalog = bundle.result.deblendedCatalog 

650 modelData = bundle.result.scarletModelData 

651 mes.io.updateCatalogFootprints( 

652 modelData, 

653 catalog, 

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

655 imageForRedistribution=None, 

656 removeScarletData=False, 

657 updateFluxColumns=True, 

658 ) 

659 

660 for src in catalog: 

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

662 

663 def test_nChild_consistency(self): 

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

665 

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

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

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

669 of exactly one multi-peak parent. 

670 """ 

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

672 catalog = bundle.result.deblendedCatalog 

673 objectParents = bundle.result.objectParents 

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

675 parents = objectParents[ 

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

677 ] 

678 self.assertEqual( 

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

680 ) 

681 

682 def test_isolated_parents_marked(self): 

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

684 parent. 

685 """ 

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

687 objectParents = bundle.result.objectParents 

688 isolatedParents = objectParents[ 

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

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

691 ] 

692 self.assertEqual( 

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

694 len(isolatedParents), 

695 ) 

696 

697 def test_isolated_source_persisted(self): 

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

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

700 """ 

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

702 catalog = bundle.result.deblendedCatalog 

703 modelData = bundle.result.scarletModelData 

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

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

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

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

708 np.testing.assert_array_equal( 

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

710 ) 

711 self.assertTupleEqual( 

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

713 ) 

714 

715 def test_heavy_footprint_flux_at_peak(self): 

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

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

718 flux redistribution. 

719 

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

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

722 shifted the flux value would still be correct but 

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

724 """ 

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

726 for useFlux in [False, True]: 

727 for band in bundle.image.bands: 

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

729 self._attach_band_footprints(bundle, band, useFlux) 

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

731 fp = child.getFootprint() 

732 img = fp.extractImage(fill=0.0) 

733 px = child.get("deblend_peak_center_x") 

734 py = child.get("deblend_peak_center_y") 

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

736 self.assertEqual( 

737 flux, child.get("deblend_peak_instFlux") 

738 ) 

739 

740 def test_heavy_footprint_band_columns_populated(self): 

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

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

743 

744 Pins both the four pre-existing writes 

745 (``deblend_zeroFlux``, ``deblend_dataCoverage``, 

746 ``deblend_scarletFlux``, ``deblend_peak_instFlux``) and the 

747 four blendedness/overlap metric writes 

748 (``deblend_maxOverlap``, ``deblend_fluxOverlap``, 

749 ``deblend_fluxOverlapFraction``, ``deblend_blendedness``). 

750 The metric writes were previously missing: their values were 

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

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

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

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

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

756 Regression test for finding DB-7 of the 

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

758 """ 

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

760 band = bundle.image.bands[0] 

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

762 

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

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

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

766 self.assertGreater(child.get("deblend_dataCoverage"), 0) 

767 self.assertGreater(child.get("deblend_scarletFlux"), 0) 

768 self.assertFalse( 

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

770 ) 

771 self.assertGreater(child.get("deblend_maxOverlap"), 0) 

772 self.assertGreater(child.get("deblend_fluxOverlap"), 0) 

773 self.assertGreater( 

774 child.get("deblend_fluxOverlapFraction"), 0 

775 ) 

776 self.assertGreater(child.get("deblend_blendedness"), 0) 

777 

778 def test_heavy_footprint_peak_position(self): 

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

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

781 

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

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

784 """ 

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

786 band = bundle.image.bands[0] 

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

788 

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

790 fp = child.getFootprint() 

791 px = child.get("deblend_peak_center_x") 

792 py = child.get("deblend_peak_center_y") 

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

794 self.assertEqual(py, fp.getPeaks()[0].getIy()) 

795 

796 _, source, _ = self._scarlet_blend_for_child( 

797 bundle, parent, child, band 

798 ) 

799 self.assertEqual(source.center[1], px) 

800 self.assertEqual(source.center[0], py) 

801 

802 def test_heavy_footprint_matches_model(self): 

803 """The HeavyFootprint pixel data matches the scarlet model. 

804 

805 Two variants per band: 

806 

807 - ``useFlux=False``: convolve the source's bare scarlet model 

808 with the band PSF and compare to the HeavyFootprint image. 

809 - ``useFlux=True``: rebind the blend's observation to the 

810 observed coadd, call ``conserve_flux`` to redistribute, and 

811 compare the flux-weighted model to the HeavyFootprint 

812 (projected onto the redistributed image's frame). 

813 """ 

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

815 image = bundle.image 

816 for useFlux in [False, True]: 

817 for band in image.bands: 

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

819 self._attach_band_footprints(bundle, band, useFlux) 

820 imageForRedistribution = ( 

821 image.mCoadd[band] if useFlux else None 

822 ) 

823 

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

825 fp = child.getFootprint() 

826 img = fp.extractImage(fill=0.0) 

827 blend, source, parentFootprint = ( 

828 self._scarlet_blend_for_child( 

829 bundle, parent, child, band