Coverage for tests/test_subtractTask.py: 99%

1044 statements  

« prev     ^ index     » next       coverage.py v7.16.2, created at 2026-10-14 19:10 +0000

1# This file is part of ip_diffim. 

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 

22import unittest 

23import uuid 

24 

25from astropy import units as u 

26 

27import lsst.afw.detection as afwDetection 

28import lsst.afw.image as afwImage 

29import lsst.afw.cameraGeom.testUtils 

30import lsst.afw.math as afwMath 

31import lsst.afw.table as afwTable 

32import lsst.geom 

33import lsst.meas.algorithms as measAlg 

34from lsst.daf.butler import DataCoordinate, DimensionUniverse 

35from lsst.images import (Box, DifferenceImage, DifferenceImageTemplateInfo, Polygon, VisitImage, 

36 get_legacy_difference_image_mask_planes) 

37from lsst.images.convolution_kernels import ImageBasisConvolutionKernel 

38from lsst.images.fields import ChebyshevField 

39from lsst.images.tests import compare_masked_image_to_legacy 

40from lsst.ip.diffim import subtractImages, InsufficientKernelSourcesError 

41from lsst.pex.config import FieldValidationError 

42from lsst.pipe.base import InvalidQuantumError, NoWorkFound 

43import lsst.utils.tests 

44import numpy as np 

45from lsst.ip.diffim.utils import (computeRobustStatistics, computePSFNoiseEquivalentArea, 

46 evaluateMeanPsfFwhm, getPsfFwhm) 

47from lsst.pex.exceptions import InvalidParameterError 

48 

49from utils import (makeStats, makeTestExposureRecord, makeTestImage, makeTestVisitInfo, 

50 CustomCoaddPsf) 

51 

52 

53class AlardLuptonSubtractTestBase: 

54 goodPsfSize = 2.0 

55 midPsfSize = 2.4 

56 badPsfSize = 2.8 

57 

58 def _setup_subtraction(self, fluxField="truth_instFlux", errField="truth_instFluxErr", **kwargs): 

59 """Setup and configure the image subtraction PipelineTask. 

60 

61 Parameters 

62 ---------- 

63 fluxField : `str`, optional 

64 Name of the flux field in the source catalog. 

65 errField : `str`, optional 

66 Name of the flux error field in the source catalog. 

67 **kwargs 

68 Any additional config parameters to set. 

69 

70 Returns 

71 ------- 

72 `lsst.pipe.base.PipelineTask` 

73 The configured Task to use for detection and measurement. 

74 """ 

75 config = self.subtractTask.ConfigClass() 

76 config.restrictKernelEdgeSources = False 

77 config.sourceSelector.signalToNoise.fluxField = fluxField 

78 config.sourceSelector.signalToNoise.errField = errField 

79 config.sourceSelector.doUnresolved = True 

80 config.sourceSelector.doIsolated = True 

81 config.sourceSelector.doRequirePrimary = True 

82 config.sourceSelector.doFlags = True 

83 config.sourceSelector.doSignalToNoise = True 

84 config.sourceSelector.flags.bad = ["base_PsfFlux_flag", ] 

85 config.update(**kwargs) 

86 

87 return self.subtractTask(config=config) 

88 

89 

90class AlardLuptonSubtractTest(AlardLuptonSubtractTestBase, lsst.utils.tests.TestCase): 

91 subtractTask = subtractImages.AlardLuptonSubtractTask 

92 

93 def test_allowed_config_modes(self): 

94 """Verify the allowable modes for convolution. 

95 """ 

96 config = subtractImages.AlardLuptonSubtractTask.ConfigClass() 

97 config.mode = 'auto' 

98 config.mode = 'convolveScience' 

99 config.mode = 'convolveTemplate' 

100 

101 with self.assertRaises(FieldValidationError): 

102 config.mode = 'aotu' 

103 

104 def test_mismatched_template(self): 

105 """Test that an error is raised if the template 

106 does not fully contain the science image. 

107 """ 

108 xSize = 200 

109 ySize = 200 

110 science, sources = makeTestImage(psfSize=self.midPsfSize, xSize=xSize + 20, ySize=ySize + 20) 

111 template, _ = makeTestImage(psfSize=self.midPsfSize, xSize=xSize, ySize=ySize, 

112 doApplyCalibration=True) 

113 task = self._setup_subtraction() 

114 with self.assertRaises(AssertionError): 

115 task.run(template, science, sources) 

116 

117 def test_mismatched_filter(self): 

118 """Test that an error is raised if the science and template have 

119 different bands. 

120 """ 

121 xSize = 200 

122 ySize = 200 

123 science, sources = makeTestImage(psfSize=self.midPsfSize, xSize=xSize + 20, ySize=ySize + 20, 

124 band="g", physicalFilter="g noCamera") 

125 template, _ = makeTestImage(psfSize=self.midPsfSize, xSize=xSize, ySize=ySize, 

126 doApplyCalibration=True, band="not-g", physicalFilter="not-g noCamera") 

127 task = self._setup_subtraction() 

128 with self.assertRaises(AssertionError): 

129 task.run(template, science, sources) 

130 

131 def test_incomplete_template_coverage(self): 

132 noiseLevel = 1. 

133 border = 20 

134 xSize = 400 

135 ySize = 400 

136 nSources = 80 

137 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6, 

138 nSrc=nSources, xSize=xSize, ySize=ySize) 

139 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7, 

140 nSrc=nSources, templateBorderSize=border, doApplyCalibration=True, 

141 xSize=xSize, ySize=ySize) 

142 

143 science_height = science.getBBox().getDimensions().getY() 

144 

145 def _run_and_check_coverage(template_coverage, 

146 requiredTemplateFraction=0.1, 

147 minTemplateFractionForExpectedSuccess=0.2): 

148 template_cut = template.clone() 

149 template_height = int(science_height*template_coverage + border) 

150 template_cut.image.array[:, template_height:] = 0 

151 template_cut.mask.array[:, template_height:] = template_cut.mask.getPlaneBitMask('NO_DATA') 

152 task = self._setup_subtraction( 

153 requiredTemplateFraction=requiredTemplateFraction, 

154 minTemplateFractionForExpectedSuccess=minTemplateFractionForExpectedSuccess 

155 ) 

156 if template_coverage < requiredTemplateFraction: 

157 doRaise = True 

158 elif template_coverage < minTemplateFractionForExpectedSuccess: 

159 doRaise = True 

160 else: 

161 doRaise = False 

162 if doRaise: 

163 with self.assertRaises(NoWorkFound): 

164 task.run(template_cut, science.clone(), sources.copy(deep=True)) 

165 else: 

166 task.run(template_cut, science.clone(), sources.copy(deep=True)) 

167 _run_and_check_coverage(template_coverage=0.09) 

168 _run_and_check_coverage(template_coverage=0.15) 

169 _run_and_check_coverage(template_coverage=0.7) 

170 

171 def test_clear_template_mask(self): 

172 noiseLevel = 1. 

173 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6) 

174 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7, 

175 templateBorderSize=20, doApplyCalibration=True) 

176 diffimEmptyMaskPlanes = ["DETECTED", "DETECTED_NEGATIVE"] 

177 task = self._setup_subtraction(mode="convolveTemplate") 

178 # Ensure that each each mask plane is set for some pixels 

179 mask = template.mask 

180 x0 = 50 

181 x1 = 75 

182 y0 = 150 

183 y1 = 175 

184 scienceMaskCheck = {} 

185 for maskPlane in mask.getMaskPlaneDict().keys(): 

186 scienceMaskCheck[maskPlane] = np.sum(science.mask.array & mask.getPlaneBitMask(maskPlane) > 0) 

187 mask.array[x0: x1, y0: y1] |= mask.getPlaneBitMask(maskPlane) 

188 self.assertTrue(np.sum(mask.array & mask.getPlaneBitMask(maskPlane) > 0)) 

189 

190 output = task.run(template, science, sources) 

191 # Verify that the template mask has been modified in place 

192 for maskPlane in mask.getMaskPlaneDict().keys(): 

193 if maskPlane in diffimEmptyMaskPlanes: 

194 self.assertTrue(np.sum(mask.array & mask.getPlaneBitMask(maskPlane) == 0)) 

195 elif maskPlane in task.config.preserveTemplateMask: 

196 self.assertTrue(np.sum(mask.array & mask.getPlaneBitMask(maskPlane) > 0)) 

197 else: 

198 self.assertTrue(np.sum(mask.array & mask.getPlaneBitMask(maskPlane) == 0)) 

199 # Mask planes set in the science image should also be set in the difference 

200 # Except the "DETECTED" planes should have been cleared 

201 diffimMask = output.difference.mask 

202 for maskPlane, scienceSum in scienceMaskCheck.items(): 

203 diffimSum = np.sum(diffimMask.array & mask.getPlaneBitMask(maskPlane) > 0) 

204 if maskPlane in diffimEmptyMaskPlanes: 

205 self.assertEqual(diffimSum, 0) 

206 else: 

207 self.assertTrue(diffimSum >= scienceSum) 

208 

209 def test_renamed_template_mask_survives_registered_plane(self): 

210 """A renamed template plane keeps its pixels whether or not the plane 

211 it is renamed to is already registered. 

212 

213 afw's mask plane dictionary is shared by every mask in the process, so 

214 reading an image that defines SAT_TEMPLATE puts that plane in the 

215 template's dictionary as well, with no pixels set. 

216 """ 

217 for preRegistered in (False, True): 

218 with self.subTest(preRegistered=preRegistered): 

219 science, _ = makeTestImage(psfSize=self.midPsfSize, noiseLevel=1., noiseSeed=6) 

220 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=1., noiseSeed=7, 

221 templateBorderSize=20, doApplyCalibration=True) 

222 template.mask.addMaskPlane("SAT") 

223 satBitMask = template.mask.getPlaneBitMask("SAT") 

224 template.mask.array[10:20, 10:20] |= satBitMask 

225 expected = np.count_nonzero(template.mask.array & satBitMask) 

226 self.assertGreater(expected, 0) 

227 

228 if preRegistered: 

229 afwImage.Mask(1, 1).addMaskPlane("SAT_TEMPLATE") 

230 

231 task = self._setup_subtraction(mode="convolveTemplate") 

232 task.updateMasks(template, science) 

233 

234 satTemplateBitMask = template.mask.getPlaneBitMask("SAT_TEMPLATE") 

235 self.assertEqual(np.count_nonzero(template.mask.array & satTemplateBitMask), expected) 

236 # The plane the pixels were renamed from is still cleared. 

237 self.assertEqual(np.count_nonzero(template.mask.array & satBitMask), 0) 

238 

239 def test_equal_images(self): 

240 """Test that running with enough sources produces reasonable output, 

241 with the same size psf in the template and science. 

242 """ 

243 noiseLevel = 1. 

244 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6) 

245 template, _ = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=7, 

246 templateBorderSize=20, doApplyCalibration=True) 

247 task = self._setup_subtraction() 

248 output = task.run(template, science, sources) 

249 # There shoud be no NaN values in the difference image 

250 self.assertTrue(np.all(np.isfinite(output.difference.image.array))) 

251 # Mean of difference image should be close to zero. 

252 meanError = noiseLevel/np.sqrt(output.difference.image.array.size) 

253 # Make sure to include pixels with the DETECTED mask bit set. 

254 statsCtrl = makeStats(badMaskPlanes=("EDGE", "BAD", "NO_DATA", "DETECTED", "DETECTED_NEGATIVE")) 

255 differenceMean = computeRobustStatistics(output.difference.image, output.difference.mask, statsCtrl) 

256 self.assertFloatsAlmostEqual(differenceMean, 0, atol=5*meanError) 

257 # stddev of difference image should be close to expected value. 

258 differenceStd = computeRobustStatistics(output.difference.image, output.difference.mask, 

259 makeStats(), statistic=afwMath.STDEV) 

260 self.assertFloatsAlmostEqual(differenceStd, np.sqrt(2)*noiseLevel, rtol=0.1) 

261 

262 def test_equal_images_missing_mask_planes(self): 

263 """Test that running with enough sources produces reasonable output, 

264 with the same size psf in the template and science and with missing 

265 mask planes. 

266 """ 

267 noiseLevel = 1. 

268 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6, 

269 addMaskPlanes=[]) 

270 template, _ = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=7, 

271 templateBorderSize=20, doApplyCalibration=True, addMaskPlanes=[]) 

272 task = self._setup_subtraction() 

273 output = task.run(template, science, sources) 

274 # There shoud be no NaN values in the difference image 

275 self.assertTrue(np.all(np.isfinite(output.difference.image.array))) 

276 # Mean of difference image should be close to zero. 

277 meanError = noiseLevel/np.sqrt(output.difference.image.array.size) 

278 # Make sure to include pixels with the DETECTED mask bit set. 

279 statsCtrl = makeStats(badMaskPlanes=("EDGE", "BAD", "NO_DATA", "DETECTED", "DETECTED_NEGATIVE")) 

280 differenceMean = computeRobustStatistics(output.difference.image, output.difference.mask, statsCtrl) 

281 self.assertFloatsAlmostEqual(differenceMean, 0, atol=5*meanError) 

282 # stddev of difference image should be close to expected value. 

283 differenceStd = computeRobustStatistics(output.difference.image, output.difference.mask, 

284 makeStats(), statistic=afwMath.STDEV) 

285 self.assertFloatsAlmostEqual(differenceStd, np.sqrt(2)*noiseLevel, rtol=0.1) 

286 

287 def test_psf_size(self): 

288 """Test that the image subtract task runs without failing, if 

289 fwhmExposureBuffer and fwhmExposureGrid parameters are set. 

290 """ 

291 noiseLevel = 1. 

292 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6) 

293 template, _ = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=7, 

294 templateBorderSize=20, doApplyCalibration=True) 

295 

296 schema = afwTable.ExposureTable.makeMinimalSchema() 

297 weightKey = schema.addField("weight", type="D", doc="Coadd weight") 

298 exposureCatalog = afwTable.ExposureCatalog(schema) 

299 kernel = measAlg.DoubleGaussianPsf(7, 7, 2.0).getKernel() 

300 psf = measAlg.KernelPsf(kernel, template.getBBox().getCenter()) 

301 

302 record = exposureCatalog.addNew() 

303 record.setPsf(psf) 

304 record.setWcs(template.wcs) 

305 record.setD(weightKey, 1.0) 

306 record.setBBox(template.getBBox()) 

307 

308 customPsf = CustomCoaddPsf(exposureCatalog, template.wcs) 

309 template.setPsf(customPsf) 

310 

311 # Test that we get an exception if we simply get the FWHM at center. 

312 with self.assertRaises(InvalidParameterError): 

313 getPsfFwhm(template.psf, True) 

314 

315 with self.assertRaises(InvalidParameterError): 

316 getPsfFwhm(template.psf, False) 

317 

318 # Test that evaluateMeanPsfFwhm runs successfully on the template. 

319 evaluateMeanPsfFwhm(template, fwhmExposureBuffer=0.05, fwhmExposureGrid=10) 

320 

321 # Since the PSF is spatially invariant, the FWHM should be the same at 

322 # all points in the science image. 

323 fwhm1 = getPsfFwhm(science.psf, False) 

324 fwhm2 = evaluateMeanPsfFwhm(science, fwhmExposureBuffer=0.05, fwhmExposureGrid=10) 

325 self.assertAlmostEqual(fwhm1[0], fwhm2, places=13) 

326 self.assertAlmostEqual(fwhm1[1], fwhm2, places=13) 

327 

328 self.assertAlmostEqual(evaluateMeanPsfFwhm(science, fwhmExposureBuffer=0.05, 

329 fwhmExposureGrid=10), 

330 getPsfFwhm(science.psf, True), places=7 

331 ) 

332 

333 # Test that the image subtraction task runs successfully. 

334 task = self._setup_subtraction() 

335 

336 # Test that the task runs if we take the mean FWHM on a grid. 

337 with self.assertLogs(level="INFO") as cm: 

338 task.run(template, science, sources) 

339 

340 # Check that evaluateMeanPsfFwhm was called. 

341 # This tests that getPsfFwhm failed raising InvalidParameterError, 

342 # that is caught and handled appropriately. 

343 logMessage = ("INFO:lsst.alardLuptonSubtract:Unable to evaluate PSF at the average position. " 

344 "Evaluting PSF on a grid of points." 

345 ) 

346 self.assertIn(logMessage, cm.output) 

347 

348 def test_auto_convolveTemplate(self): 

349 """Test that auto mode gives the same result as convolveTemplate when 

350 the template psf is the smaller. 

351 """ 

352 noiseLevel = 1. 

353 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6) 

354 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7, 

355 templateBorderSize=20, doApplyCalibration=True) 

356 task = self._setup_subtraction(mode="convolveTemplate") 

357 output = task.run(template.clone(), science.clone(), sources) 

358 

359 task = self._setup_subtraction(mode="auto") 

360 outputAuto = task.run(template, science, sources) 

361 self.assertMaskedImagesEqual(output.difference.maskedImage, outputAuto.difference.maskedImage) 

362 

363 def test_auto_convolveScience(self): 

364 """Test that auto mode gives the same result as convolveScience when 

365 the science psf is the smaller. 

366 """ 

367 noiseLevel = 1. 

368 science, sources = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=6) 

369 template, _ = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=7, 

370 templateBorderSize=20, doApplyCalibration=True) 

371 task = self._setup_subtraction(mode="convolveScience") 

372 output = task.run(template.clone(), science.clone(), sources) 

373 

374 task = self._setup_subtraction(mode="auto") 

375 outputAuto = task.run(template, science, sources) 

376 self.assertMaskedImagesEqual(output.difference.maskedImage, outputAuto.difference.maskedImage) 

377 

378 def test_science_better(self): 

379 """Test that running with enough sources produces reasonable output, 

380 with the science psf being smaller than the template. 

381 """ 

382 statsCtrl = makeStats() 

383 statsCtrlDetect = makeStats(badMaskPlanes=("EDGE", "BAD", "NO_DATA")) 

384 

385 def _run_and_check_images(statsCtrl, statsCtrlDetect, scienceNoiseLevel, templateNoiseLevel): 

386 science, sources = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=scienceNoiseLevel, 

387 noiseSeed=6) 

388 template, _ = makeTestImage(psfSize=self.midPsfSize, noiseLevel=templateNoiseLevel, noiseSeed=7, 

389 templateBorderSize=20, doApplyCalibration=True) 

390 task = self._setup_subtraction(mode="convolveScience") 

391 output = task.run(template, science, sources) 

392 self.assertFloatsAlmostEqual(task.metadata["scaleScienceVarianceFactor"], 1., atol=.05) 

393 # Mean of difference image should be close to zero. 

394 nGoodPix = np.sum(np.isfinite(output.difference.image.array)) 

395 meanError = (scienceNoiseLevel + templateNoiseLevel)/np.sqrt(nGoodPix) 

396 diffimMean = computeRobustStatistics(output.difference.image, output.difference.mask, 

397 statsCtrlDetect) 

398 

399 self.assertFloatsAlmostEqual(diffimMean, 0, atol=5*meanError) 

400 # stddev of difference image should be close to expected value. 

401 noiseLevel = np.sqrt(scienceNoiseLevel**2 + templateNoiseLevel**2) 

402 varianceMean = computeRobustStatistics(output.difference.variance, output.difference.mask, 

403 statsCtrl) 

404 diffimStd = computeRobustStatistics(output.difference.image, output.difference.mask, 

405 statsCtrl, statistic=afwMath.STDEV) 

406 self.assertFloatsAlmostEqual(varianceMean, noiseLevel**2, rtol=0.1) 

407 self.assertFloatsAlmostEqual(diffimStd, noiseLevel, rtol=0.1) 

408 

409 _run_and_check_images(statsCtrl, statsCtrlDetect, scienceNoiseLevel=1., templateNoiseLevel=1.) 

410 _run_and_check_images(statsCtrl, statsCtrlDetect, scienceNoiseLevel=1., templateNoiseLevel=.1) 

411 _run_and_check_images(statsCtrl, statsCtrlDetect, scienceNoiseLevel=.1, templateNoiseLevel=.1) 

412 

413 def test_template_better(self): 

414 """Test that running with enough sources produces reasonable output, 

415 with the template psf being smaller than the science. 

416 """ 

417 statsCtrl = makeStats() 

418 statsCtrlDetect = makeStats(badMaskPlanes=("EDGE", "BAD", "NO_DATA")) 

419 

420 def _run_and_check_images(statsCtrl, statsCtrlDetect, scienceNoiseLevel, templateNoiseLevel): 

421 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=scienceNoiseLevel, 

422 noiseSeed=6) 

423 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=templateNoiseLevel, noiseSeed=7, 

424 templateBorderSize=20, doApplyCalibration=True) 

425 task = self._setup_subtraction() 

426 output = task.run(template, science, sources) 

427 self.assertFloatsAlmostEqual(task.metadata["scaleScienceVarianceFactor"], 1., atol=.05) 

428 # There should be no NaNs in the image if we convolve the template with a buffer 

429 self.assertTrue(np.all(np.isfinite(output.difference.image.array))) 

430 # Mean of difference image should be close to zero. 

431 meanError = (scienceNoiseLevel + templateNoiseLevel)/np.sqrt(output.difference.image.array.size) 

432 

433 diffimMean = computeRobustStatistics(output.difference.image, output.difference.mask, 

434 statsCtrlDetect) 

435 self.assertFloatsAlmostEqual(diffimMean, 0, atol=5*meanError) 

436 # stddev of difference image should be close to expected value. 

437 noiseLevel = np.sqrt(scienceNoiseLevel**2 + templateNoiseLevel**2) 

438 varianceMean = computeRobustStatistics(output.difference.variance, output.difference.mask, 

439 statsCtrl) 

440 diffimStd = computeRobustStatistics(output.difference.image, output.difference.mask, 

441 statsCtrl, statistic=afwMath.STDEV) 

442 self.assertFloatsAlmostEqual(varianceMean, noiseLevel**2, rtol=0.1) 

443 self.assertFloatsAlmostEqual(diffimStd, noiseLevel, rtol=0.1) 

444 

445 _run_and_check_images(statsCtrl, statsCtrlDetect, scienceNoiseLevel=1., templateNoiseLevel=1.) 

446 _run_and_check_images(statsCtrl, statsCtrlDetect, scienceNoiseLevel=1., templateNoiseLevel=.1) 

447 _run_and_check_images(statsCtrl, statsCtrlDetect, scienceNoiseLevel=.1, templateNoiseLevel=.1) 

448 

449 def test_symmetry(self): 

450 """Test that convolving the science and convolving the template are 

451 symmetric: if the psfs are switched between them, the difference image 

452 should be nearly the same. 

453 """ 

454 noiseLevel = 1. 

455 # Don't include a border for the template, in order to make the results 

456 # comparable when we swap which image is treated as the "science" image. 

457 science, sources = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, 

458 noiseSeed=6, templateBorderSize=0, doApplyCalibration=True) 

459 template, _ = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, 

460 noiseSeed=7, templateBorderSize=0, doApplyCalibration=True) 

461 task = self._setup_subtraction(mode='auto') 

462 

463 # The science image will be modified in place, so use a copy for the second run. 

464 science_better = task.run(template.clone(), science.clone(), sources) 

465 template_better = task.run(science, template, sources) 

466 

467 delta = template_better.difference.clone() 

468 delta.image -= science_better.difference.image 

469 delta.variance -= science_better.difference.variance 

470 delta.mask.array -= science_better.difference.mask.array 

471 

472 statsCtrl = makeStats() 

473 # Mean of delta should be very close to zero. 

474 nGoodPix = np.sum(np.isfinite(delta.image.array)) 

475 meanError = 2*noiseLevel/np.sqrt(nGoodPix) 

476 deltaMean = computeRobustStatistics(delta.image, delta.mask, statsCtrl) 

477 deltaStd = computeRobustStatistics(delta.image, delta.mask, statsCtrl, statistic=afwMath.STDEV) 

478 self.assertFloatsAlmostEqual(deltaMean, 0, atol=5*meanError) 

479 # stddev of difference image should be close to expected value 

480 self.assertFloatsAlmostEqual(deltaStd, 2*np.sqrt(2)*noiseLevel, rtol=.1) 

481 

482 def test_few_sources(self): 

483 """Test with only 1 source, to check that we get a useful error. 

484 """ 

485 xSize = 256 

486 ySize = 256 

487 science, sources = makeTestImage(psfSize=self.midPsfSize, nSrc=10, xSize=xSize, ySize=ySize) 

488 template, _ = makeTestImage(psfSize=self.goodPsfSize, nSrc=10, xSize=xSize, ySize=ySize, 

489 doApplyCalibration=True) 

490 task = self._setup_subtraction() 

491 sources = sources[0:1] 

492 with self.assertRaises(InsufficientKernelSourcesError): 

493 task.run(template, science, sources) 

494 

495 def test_kernel_source_selector(self): 

496 """Check that kernel source selection behaves as expected. 

497 """ 

498 xSize = 256 

499 ySize = 256 

500 nSourcesSimulated = 20 

501 sciencePsfSize = self.midPsfSize 

502 templatePsfSize = self.goodPsfSize 

503 science, sources = makeTestImage(psfSize=sciencePsfSize, nSrc=nSourcesSimulated, 

504 xSize=xSize, ySize=ySize) 

505 template, _ = makeTestImage(psfSize=templatePsfSize, nSrc=nSourcesSimulated, 

506 xSize=xSize, ySize=ySize, doApplyCalibration=True) 

507 

508 def _run_and_check_sources(sourcesIn, maxKernelSources=1000, minKernelSources=3): 

509 sources = sourcesIn.copy(deep=True) 

510 

511 task = self._setup_subtraction(maxKernelSources=maxKernelSources, 

512 minKernelSources=minKernelSources, 

513 ) 

514 task.templatePsfSize = templatePsfSize 

515 task.sciencePsfSize = sciencePsfSize 

516 task.matchedPsfSize = sciencePsfSize 

517 # Verify that source flags are not set in the input catalog 

518 # Note that this will use the last flag in the list for the rest of 

519 # the test. 

520 for badSourceFlag in task.sourceSelector.config.flags.bad: 

521 self.assertEqual(np.sum(sources[badSourceFlag]), 0) 

522 nSources = len(sources) 

523 # Flag a third of the sources 

524 sources[0:: 3][badSourceFlag] = True 

525 rejectRadius = 2*task.config.makeKernel.kernel.active.kernelSize 

526 bbox = science.getBBox() 

527 bbox.grow(-rejectRadius) 

528 edgeSources = ~bbox.contains(sources.getX(), sources.getY()) 

529 sources[edgeSources][badSourceFlag] = True 

530 nBadSources = np.sum(sources[badSourceFlag]) 

531 if maxKernelSources > 0: 

532 nGoodSources = np.minimum(nSources - nBadSources, maxKernelSources) 

533 else: 

534 nGoodSources = nSources - nBadSources 

535 

536 signalToNoise = sources.getPsfInstFlux()/sources.getPsfInstFluxErr() 

537 signalToNoise = signalToNoise[~sources[badSourceFlag]] 

538 signalToNoise.sort() 

539 selectSources = task._sourceSelector(template, science, sources) 

540 self.assertEqual(nGoodSources, len(selectSources)) 

541 signalToNoiseOut = selectSources.getPsfInstFlux()/selectSources.getPsfInstFluxErr() 

542 signalToNoiseOut.sort() 

543 self.assertFloatsAlmostEqual(signalToNoise[-nGoodSources:], signalToNoiseOut) 

544 

545 _run_and_check_sources(sources) 

546 _run_and_check_sources(sources, maxKernelSources=len(sources)//3) 

547 _run_and_check_sources(sources, maxKernelSources=-1) 

548 with self.assertRaises(RuntimeError): 

549 _run_and_check_sources(sources, minKernelSources=1000) 

550 

551 def test_order_equal_images(self): 

552 """Verify that the result is the same regardless of convolution mode 

553 if the images are equivalent. 

554 """ 

555 noiseLevel = .1 

556 seed1 = 6 

557 seed2 = 7 

558 for psfSize in [self.midPsfSize, self.goodPsfSize, self.badPsfSize]: 

559 science1, sources1 = makeTestImage(psfSize=psfSize, noiseLevel=noiseLevel, noiseSeed=seed1, 

560 clearEdgeMask=True) 

561 template1, _ = makeTestImage(psfSize=psfSize, noiseLevel=noiseLevel, noiseSeed=seed2, 

562 templateBorderSize=0, doApplyCalibration=True, 

563 clearEdgeMask=True) 

564 task1 = self._setup_subtraction(mode="convolveTemplate") 

565 results_convolveTemplate = task1.run(template1, science1, sources1) 

566 

567 science2, sources2 = makeTestImage(psfSize=psfSize, noiseLevel=noiseLevel, noiseSeed=seed1, 

568 clearEdgeMask=True) 

569 template2, _ = makeTestImage(psfSize=psfSize, noiseLevel=noiseLevel, noiseSeed=seed2, 

570 templateBorderSize=0, doApplyCalibration=True, 

571 clearEdgeMask=True) 

572 task2 = self._setup_subtraction(mode="convolveScience") 

573 results_convolveScience = task2.run(template2, science2, sources2) 

574 bbox = results_convolveTemplate.difference.getBBox().clippedTo( 

575 results_convolveScience.difference.getBBox()) 

576 diff1 = science1.maskedImage.clone()[bbox] 

577 diff1 -= template1.maskedImage[bbox] 

578 diff2 = science2.maskedImage.clone()[bbox] 

579 diff2 -= template2.maskedImage[bbox] 

580 self.assertFloatsAlmostEqual(results_convolveTemplate.difference[bbox].image.array, 

581 diff1.image.array, 

582 atol=noiseLevel*5.) 

583 self.assertFloatsAlmostEqual(results_convolveScience.difference[bbox].image.array, 

584 diff2.image.array, 

585 atol=noiseLevel*5.) 

586 diffErr = noiseLevel*2 

587 self.assertMaskedImagesAlmostEqual(results_convolveTemplate.difference[bbox].maskedImage, 

588 results_convolveScience.difference[bbox].maskedImage, 

589 atol=diffErr*5.) 

590 

591 def test_background_subtraction(self): 

592 """Check that we can recover the background, 

593 and that it is subtracted correctly in the difference image. 

594 

595 NOTE: Background subtraction is now turned off by default in 

596 subtractImages. It is now run in detectAndMeasure instead, but since the 

597 code to run background subtraction is not being removed this test should 

598 stay to make sure it continues functioning as intended. 

599 """ 

600 noiseLevel = 1. 

601 xSize = 512 

602 ySize = 512 

603 x0 = 123 

604 y0 = 456 

605 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7, 

606 templateBorderSize=20, 

607 xSize=xSize, ySize=ySize, x0=x0, y0=y0, 

608 doApplyCalibration=True) 

609 params = [2.2, 2.1, 2.0, 1.2, 1.1, 1.0] 

610 

611 bbox2D = lsst.geom.Box2D(lsst.geom.Point2D(x0, y0), lsst.geom.Extent2D(xSize, ySize)) 

612 background_model = afwMath.Chebyshev1Function2D(params, bbox2D) 

613 science, sources = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=6, 

614 background=background_model, 

615 xSize=xSize, ySize=ySize, x0=x0, y0=y0) 

616 # Don't use ``self._setup_subtraction()`` here. 

617 # Modifying the config of a subtask is messy. 

618 config = subtractImages.AlardLuptonSubtractTask.ConfigClass() 

619 

620 config.sourceSelector.signalToNoise.fluxField = "truth_instFlux" 

621 config.sourceSelector.signalToNoise.errField = "truth_instFluxErr" 

622 config.doSubtractBackground = True 

623 

624 config.makeKernel.kernel.name = "AL" 

625 config.makeKernel.kernel.active.fitForBackground = True 

626 config.makeKernel.kernel.active.spatialKernelOrder = 1 

627 config.makeKernel.kernel.active.spatialBgOrder = 2 

628 statsCtrl = makeStats() 

629 

630 def _run_and_check_images(config, statsCtrl, mode): 

631 """Check that the fit background matches the input model. 

632 """ 

633 config.mode = mode 

634 task = subtractImages.AlardLuptonSubtractTask(config=config) 

635 output = task.run(template.clone(), science.clone(), sources) 

636 

637 # We should be fitting the same number of parameters as were in the input model 

638 self.assertEqual(output.backgroundModel.getNParameters(), background_model.getNParameters()) 

639 

640 # The parameters of the background fit should be close to the input model 

641 self.assertFloatsAlmostEqual(np.array(output.backgroundModel.getParameters()), 

642 np.array(params), rtol=0.3) 

643 

644 # stddev of difference image should be close to expected value. 

645 # This will fail if we have mis-subtracted the background. 

646 stdVal = computeRobustStatistics(output.difference.image, output.difference.mask, 

647 statsCtrl, statistic=afwMath.STDEV) 

648 self.assertFloatsAlmostEqual(stdVal, np.sqrt(2)*noiseLevel, rtol=0.1) 

649 

650 _run_and_check_images(config, statsCtrl, "convolveTemplate") 

651 _run_and_check_images(config, statsCtrl, "convolveScience") 

652 

653 def test_scale_variance_convolve_template(self): 

654 """Check variance scaling of the image difference. 

655 """ 

656 scienceNoiseLevel = 4. 

657 templateNoiseLevel = 2. 

658 scaleFactor = 1.345 

659 # Make sure to include pixels with the DETECTED mask bit set. 

660 statsCtrl = makeStats(badMaskPlanes=("EDGE", "BAD", "NO_DATA")) 

661 

662 def _run_and_check_images(science, template, sources, statsCtrl, 

663 doDecorrelation, doScaleVariance, scaleFactor=1.): 

664 """Check that the variance plane matches the expected value for 

665 different configurations of ``doDecorrelation`` and ``doScaleVariance``. 

666 """ 

667 

668 task = self._setup_subtraction(doDecorrelation=doDecorrelation, 

669 doScaleVariance=doScaleVariance, 

670 ) 

671 output = task.run(template.clone(), science.clone(), sources) 

672 if doScaleVariance: 

673 self.assertFloatsAlmostEqual(task.metadata["scaleScienceVarianceFactor"], 

674 scaleFactor, atol=0.05) 

675 

676 scienceNoise = computeRobustStatistics(science.variance, science.mask, statsCtrl) 

677 if doDecorrelation: 

678 templateNoise = computeRobustStatistics(template.variance, template.mask, statsCtrl) 

679 else: 

680 templateNoise = computeRobustStatistics(output.matchedTemplate.variance, 

681 output.matchedTemplate.mask, 

682 statsCtrl) 

683 

684 if doScaleVariance: 

685 # Only the science variance is scaled here. The template 

686 # variance is scaled independently in ``GetTemplateTask``. 

687 scienceNoise *= scaleFactor 

688 varMean = computeRobustStatistics(output.difference.variance, output.difference.mask, statsCtrl) 

689 self.assertFloatsAlmostEqual(varMean, scienceNoise + templateNoise, rtol=0.1) 

690 

691 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=scienceNoiseLevel, noiseSeed=6) 

692 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=templateNoiseLevel, noiseSeed=7, 

693 templateBorderSize=20, doApplyCalibration=True) 

694 # Verify that the variance plane of the difference image is correct 

695 # when the template and science variance planes are correct 

696 _run_and_check_images(science, template, sources, statsCtrl, 

697 doDecorrelation=True, doScaleVariance=True) 

698 _run_and_check_images(science, template, sources, statsCtrl, 

699 doDecorrelation=True, doScaleVariance=False) 

700 _run_and_check_images(science, template, sources, statsCtrl, 

701 doDecorrelation=False, doScaleVariance=True) 

702 _run_and_check_images(science, template, sources, statsCtrl, 

703 doDecorrelation=False, doScaleVariance=False) 

704 

705 # Verify that the variance plane of the difference image is correct 

706 # when the input science variance plane is incorrect 

707 science.variance.array /= scaleFactor 

708 _run_and_check_images(science, template, sources, statsCtrl, 

709 doDecorrelation=True, doScaleVariance=True, scaleFactor=scaleFactor) 

710 _run_and_check_images(science, template, sources, statsCtrl, 

711 doDecorrelation=True, doScaleVariance=False, scaleFactor=scaleFactor) 

712 _run_and_check_images(science, template, sources, statsCtrl, 

713 doDecorrelation=False, doScaleVariance=True, scaleFactor=scaleFactor) 

714 _run_and_check_images(science, template, sources, statsCtrl, 

715 doDecorrelation=False, doScaleVariance=False, scaleFactor=scaleFactor) 

716 

717 def test_scale_variance_convolve_science(self): 

718 """Check variance scaling of the image difference. 

719 """ 

720 scienceNoiseLevel = 4. 

721 templateNoiseLevel = 2. 

722 scaleFactor = 1.345 

723 # Make sure to include pixels with the DETECTED mask bit set. 

724 statsCtrl = makeStats(badMaskPlanes=("EDGE", "BAD", "NO_DATA")) 

725 

726 def _run_and_check_images(science, template, sources, statsCtrl, 

727 doDecorrelation, doScaleVariance, scaleFactor=1.): 

728 """Check that the variance plane matches the expected value for 

729 different configurations of ``doDecorrelation`` and ``doScaleVariance``. 

730 """ 

731 

732 task = self._setup_subtraction(mode="convolveScience", 

733 doDecorrelation=doDecorrelation, 

734 doScaleVariance=doScaleVariance, 

735 restrictKernelEdgeSources=False, 

736 ) 

737 output = task.run(template.clone(), science.clone(), sources) 

738 if doScaleVariance: 

739 self.assertFloatsAlmostEqual(task.metadata["scaleScienceVarianceFactor"], 

740 scaleFactor, atol=0.05) 

741 

742 templateNoise = computeRobustStatistics(template.variance, template.mask, statsCtrl) 

743 if doDecorrelation: 

744 scienceNoise = computeRobustStatistics(science.variance, science.mask, statsCtrl) 

745 else: 

746 scienceNoise = computeRobustStatistics(output.matchedScience.variance, 

747 output.matchedScience.mask, 

748 statsCtrl) 

749 

750 if doScaleVariance: 

751 # Only the science variance is scaled here. The template 

752 # variance is scaled independently in ``GetTemplateTask``. 

753 scienceNoise *= scaleFactor 

754 

755 varMean = computeRobustStatistics(output.difference.variance, output.difference.mask, statsCtrl) 

756 self.assertFloatsAlmostEqual(varMean, scienceNoise + templateNoise, rtol=0.1) 

757 

758 science, sources = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=scienceNoiseLevel, noiseSeed=6) 

759 template, _ = makeTestImage(psfSize=self.midPsfSize, noiseLevel=templateNoiseLevel, noiseSeed=7, 

760 templateBorderSize=20, doApplyCalibration=True) 

761 # Verify that the variance plane of the difference image is correct 

762 # when the template and science variance planes are correct 

763 _run_and_check_images(science, template, sources, statsCtrl, 

764 doDecorrelation=True, doScaleVariance=True) 

765 _run_and_check_images(science, template, sources, statsCtrl, 

766 doDecorrelation=True, doScaleVariance=False) 

767 _run_and_check_images(science, template, sources, statsCtrl, 

768 doDecorrelation=False, doScaleVariance=True) 

769 _run_and_check_images(science, template, sources, statsCtrl, 

770 doDecorrelation=False, doScaleVariance=False) 

771 

772 # Verify that the variance plane of the difference image is correct 

773 # when the input science variance plane is incorrect 

774 science.variance.array /= scaleFactor 

775 _run_and_check_images(science, template, sources, statsCtrl, 

776 doDecorrelation=True, doScaleVariance=True, scaleFactor=scaleFactor) 

777 _run_and_check_images(science, template, sources, statsCtrl, 

778 doDecorrelation=True, doScaleVariance=False, scaleFactor=scaleFactor) 

779 _run_and_check_images(science, template, sources, statsCtrl, 

780 doDecorrelation=False, doScaleVariance=True, scaleFactor=scaleFactor) 

781 _run_and_check_images(science, template, sources, statsCtrl, 

782 doDecorrelation=False, doScaleVariance=False, scaleFactor=scaleFactor) 

783 

784 def test_exposure_properties_convolve_template(self): 

785 """Check that all necessary exposure metadata is included 

786 when the template is convolved. 

787 """ 

788 noiseLevel = 1. 

789 seed = 37 

790 rng = np.random.RandomState(seed) 

791 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6) 

792 psf = science.psf 

793 psfAvgPos = psf.getAveragePosition() 

794 psfSize = getPsfFwhm(science.psf) 

795 psfImg = psf.computeKernelImage(psfAvgPos) 

796 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7, 

797 templateBorderSize=20, doApplyCalibration=True) 

798 

799 # Generate a random aperture correction map 

800 apCorrMap = lsst.afw.image.ApCorrMap() 

801 for name in ("a", "b", "c"): 

802 apCorrMap.set(name, lsst.afw.math.ChebyshevBoundedField(science.getBBox(), rng.randn(3, 3))) 

803 science.info.setApCorrMap(apCorrMap) 

804 

805 def _run_and_check_images(doDecorrelation): 

806 """Check that the metadata is correct with or without decorrelation. 

807 """ 

808 task = self._setup_subtraction(mode="convolveTemplate", 

809 doDecorrelation=doDecorrelation, 

810 ) 

811 output = task.run(template.clone(), science.clone(), sources) 

812 psfOut = output.difference.psf 

813 psfAvgPos = psfOut.getAveragePosition() 

814 if doDecorrelation: 

815 # Decorrelation requires recalculating the PSF, 

816 # so it will not be the same as the input 

817 psfOutSize = getPsfFwhm(science.psf) 

818 self.assertFloatsAlmostEqual(psfSize, psfOutSize) 

819 else: 

820 psfOutImg = psfOut.computeKernelImage(psfAvgPos) 

821 self.assertImagesAlmostEqual(psfImg, psfOutImg) 

822 

823 # check PSF, WCS, bbox, filterLabel, photoCalib, aperture correction 

824 self._compare_apCorrMaps(apCorrMap, output.difference.info.getApCorrMap()) 

825 self.assertWcsAlmostEqualOverBBox(science.wcs, output.difference.wcs, science.getBBox()) 

826 self.assertEqual(science.filter, output.difference.filter) 

827 self.assertEqual(science.photoCalib, output.difference.photoCalib) 

828 _run_and_check_images(doDecorrelation=True) 

829 _run_and_check_images(doDecorrelation=False) 

830 

831 def test_exposure_properties_convolve_science(self): 

832 """Check that all necessary exposure metadata is included 

833 when the science image is convolved. 

834 """ 

835 noiseLevel = 1. 

836 seed = 37 

837 rng = np.random.RandomState(seed) 

838 science, sources = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=6) 

839 template, _ = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=7, 

840 templateBorderSize=20, doApplyCalibration=True) 

841 psf = template.psf 

842 psfAvgPos = psf.getAveragePosition() 

843 psfSize = getPsfFwhm(template.psf) 

844 psfImg = psf.computeKernelImage(psfAvgPos) 

845 

846 # Generate a random aperture correction map 

847 apCorrMap = lsst.afw.image.ApCorrMap() 

848 for name in ("a", "b", "c"): 

849 apCorrMap.set(name, lsst.afw.math.ChebyshevBoundedField(science.getBBox(), rng.randn(3, 3))) 

850 science.info.setApCorrMap(apCorrMap) 

851 

852 def _run_and_check_images(doDecorrelation): 

853 """Check that the metadata is correct with or without decorrelation. 

854 """ 

855 task = self._setup_subtraction(mode="convolveScience", 

856 doDecorrelation=doDecorrelation, 

857 ) 

858 output = task.run(template.clone(), science.clone(), sources) 

859 if doDecorrelation: 

860 # Decorrelation requires recalculating the PSF, 

861 # so it will not be the same as the input 

862 psfOutSize = getPsfFwhm(template.psf) 

863 self.assertFloatsAlmostEqual(psfSize, psfOutSize) 

864 else: 

865 psfOut = output.difference.psf 

866 psfAvgPos = psfOut.getAveragePosition() 

867 psfOutImg = psfOut.computeKernelImage(psfAvgPos) 

868 self.assertImagesAlmostEqual(psfImg, psfOutImg) 

869 

870 # check PSF, WCS, bbox, filterLabel, photoCalib, aperture correction 

871 self._compare_apCorrMaps(apCorrMap, output.difference.info.getApCorrMap()) 

872 self.assertWcsAlmostEqualOverBBox(science.wcs, output.difference.wcs, science.getBBox()) 

873 self.assertEqual(science.filter, output.difference.filter) 

874 self.assertEqual(science.photoCalib, output.difference.photoCalib) 

875 

876 _run_and_check_images(doDecorrelation=True) 

877 _run_and_check_images(doDecorrelation=False) 

878 

879 def _compare_apCorrMaps(self, a, b): 

880 """Compare two ApCorrMaps for equality, without assuming that their BoundedFields have the 

881 same addresses (i.e. so we can compare after serialization). 

882 

883 This function is taken from ``ApCorrMapTestCase`` in afw/tests/. 

884 

885 Parameters 

886 ---------- 

887 a, b : `lsst.afw.image.ApCorrMap` 

888 The two aperture correction maps to compare. 

889 """ 

890 self.assertEqual(len(a), len(b)) 

891 for name, value in list(a.items()): 

892 value2 = b.get(name) 

893 self.assertIsNotNone(value2) 

894 self.assertEqual(value.getBBox(), value2.getBBox()) 

895 self.assertFloatsAlmostEqual( 

896 value.getCoefficients(), value2.getCoefficients(), rtol=0.0) 

897 

898 def test_fake_mask_plane_propagation(self): 

899 """Test that we have the mask planes related to fakes in diffim images. 

900 This is testing method called updateMasks 

901 """ 

902 xSize = 200 

903 ySize = 200 

904 science, sources = makeTestImage(psfSize=self.midPsfSize, xSize=xSize, ySize=ySize) 

905 science_fake_img, science_fake_sources = makeTestImage( 

906 psfSize=self.midPsfSize, xSize=xSize, ySize=ySize, seed=7, nSrc=2, noiseLevel=0.25, fluxRange=1 

907 ) 

908 template, _ = makeTestImage(psfSize=self.midPsfSize, xSize=xSize, ySize=ySize, 

909 doApplyCalibration=True) 

910 tmplt_fake_img, tmplt_fake_sources = makeTestImage( 

911 psfSize=self.midPsfSize, xSize=xSize, ySize=ySize, seed=9, nSrc=2, noiseLevel=0.25, fluxRange=1 

912 ) 

913 # created fakes and added them to the images 

914 science.image.array += science_fake_img.image.array 

915 template.image.array += tmplt_fake_img.image.array 

916 

917 # TODO: DM-40796 update to INJECTED names when source injection gets refactored 

918 # adding mask planes to both science and template images 

919 science_mask_planes = science.mask.addMaskPlane("FAKE") 

920 template_mask_planes = template.mask.addMaskPlane("FAKE") 

921 

922 for a_science_source in science_fake_sources: 

923 # 3 x 3 masking of the source locations is fine 

924 bbox = lsst.geom.Box2I( 

925 lsst.geom.Point2I(a_science_source.getX(), a_science_source.getY()), lsst.geom.Extent2I(3, 3) 

926 ) 

927 science[bbox].mask.array |= science_mask_planes 

928 

929 for a_template_source in tmplt_fake_sources: 

930 # 3 x 3 masking of the source locations is fine 

931 bbox = lsst.geom.Box2I( 

932 lsst.geom.Point2I(a_template_source.getX(), a_template_source.getY()), 

933 lsst.geom.Extent2I(3, 3) 

934 ) 

935 template[bbox].mask.array |= template_mask_planes 

936 

937 science_fake_masked = (science.mask.array & science.mask.getPlaneBitMask("FAKE")) > 0 

938 template_fake_masked = (template.mask.array & template.mask.getPlaneBitMask("FAKE")) > 0 

939 

940 task = self._setup_subtraction() 

941 subtraction = task.run(template, science, sources) 

942 

943 # check subtraction mask plane is set where we set the previous masks 

944 diff_mask = subtraction.difference.mask 

945 

946 # science mask should be now in INJECTED 

947 inj_masked = (diff_mask.array & diff_mask.getPlaneBitMask("INJECTED")) > 0 

948 

949 # template mask should be now in INJECTED_TEMPLATE 

950 injTmplt_masked = (diff_mask.array & diff_mask.getPlaneBitMask("INJECTED_TEMPLATE")) > 0 

951 

952 self.assertEqual(np.sum(inj_masked.astype(int)-science_fake_masked.astype(int)), 0) 

953 self.assertEqual(np.sum(injTmplt_masked.astype(int)-template_fake_masked.astype(int)), 0) 

954 

955 def test_calculateMagLim(self): 

956 """Test that the limiting magnitude matches the analytic value for a 

957 Gaussian PSF and flat noise, and does not depend on the pixel units 

958 or on masked pixels. 

959 """ 

960 sigma = 1.5 

961 psfSigma = 2.0 

962 calibration = 3.0 

963 bbox = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Extent2I(200, 200)) 

964 exposure = afwImage.ExposureF(bbox) 

965 exposure.variance.array[:, :] = sigma**2 

966 exposure.setPsf(afwDetection.GaussianPsf(41, 41, psfSigma)) 

967 exposure.setPhotoCalib(afwImage.PhotoCalib(calibration)) 

968 task = self._setup_subtraction() 

969 

970 maglim = task._calculateMagLim(exposure) 

971 fluxLim = 5*sigma*np.sqrt(4*np.pi*psfSigma**2) 

972 self.assertFloatsAlmostEqual(maglim, afwImage.PhotoCalib(calibration).instFluxToMagnitude(fluxLim), 

973 atol=1e-3) 

974 

975 # Pixels in other units, with a matching calibration. 

976 scaled = exposure.clone() 

977 scaled.maskedImage *= 10.0 

978 scaled.setPhotoCalib(afwImage.PhotoCalib(calibration/10.0)) 

979 self.assertFloatsAlmostEqual(task._calculateMagLim(scaled), maglim, atol=1e-6) 

980 

981 # Masked pixels do not count toward the noise. 

982 masked = exposure.clone() 

983 masked.variance.array[:, :100] = 1e6 

984 masked.mask.array[:, :100] |= masked.mask.getPlaneBitMask("NO_DATA") 

985 self.assertFloatsAlmostEqual(task._calculateMagLim(masked), maglim, atol=1e-6) 

986 

987 masked.mask.array[:, :] |= masked.mask.getPlaneBitMask("NO_DATA") 

988 self.assertTrue(np.isnan(task._calculateMagLim(masked))) 

989 

990 def test_metadata_metrics(self): 

991 """Verify fields are added to metadata when subtraction is run, and 

992 that the difference image limiting magnitude is calculated correctly, 

993 both with a "good" and "bad" seeing template. 

994 """ 

995 science, sources = makeTestImage(psfSize=self.badPsfSize, noiseLevel=1) 

996 template_good, _ = makeTestImage(psfSize=self.midPsfSize, doApplyCalibration=True, noiseLevel=0.25, 

997 templateBorderSize=20) 

998 template_bad, _ = makeTestImage(psfSize=9.5, doApplyCalibration=True, noiseLevel=0.25, 

999 templateBorderSize=20) 

1000 

1001 # Add a few sky objects; sky footprints are needed for some metrics. 

1002 config = measAlg.SkyObjectsTask.ConfigClass() 

1003 config.nSources = 3 

1004 skyTask = measAlg.SkyObjectsTask(config=config, name="skySources") 

1005 skyTask.skySourceKey = sources.schema["sky_source"].asKey() 

1006 skyTask.run(science.mask, 10, catalog=sources) 

1007 sources = sources.copy(deep=True) 

1008 # Add centroids, since these sources were added post-measurement. 

1009 for record in sources[sources["sky_source"]]: 

1010 record["truth_x"] = record.getFootprint().getPeaks()[0].getFx() 

1011 record["truth_y"] = record.getFootprint().getPeaks()[0].getFy() 

1012 

1013 # The metadata fields are attached to the subtractTask, so we do 

1014 # need to run that; run it for both "good" and "bad" seeing templates 

1015 

1016 subtractTask_good = self._setup_subtraction() 

1017 _ = subtractTask_good.run(template_good.clone(), science.clone(), sources) 

1018 subtractTask_bad = self._setup_subtraction() 

1019 _ = subtractTask_bad.run(template_bad.clone(), science.clone(), sources) 

1020 

1021 # Test that the diffim limiting magnitude combines the science and 

1022 # template limits in flux. The limits depend on the variance and mask 

1023 # planes, which `run` modifies, so take them from the metadata. 

1024 def combineMagLims(metadata): 

1025 fluxlim_science = (metadata['scienceLimitingMagnitude']*u.ABmag).to_value(u.nJy) 

1026 fluxlim_template = (metadata['templateLimitingMagnitude']*u.ABmag).to_value(u.nJy) 

1027 return (np.hypot(fluxlim_science, fluxlim_template)*u.nJy).to(u.ABmag).value 

1028 

1029 for task in (subtractTask_good, subtractTask_bad): 

1030 self.assertFloatsAlmostEqual(task.metadata['diffimLimitingMagnitude'], 

1031 combineMagLims(task.metadata), atol=1e-6) 

1032 # The template with the wider PSF is shallower. 

1033 self.assertLess(subtractTask_bad.metadata['templateLimitingMagnitude'], 

1034 subtractTask_good.metadata['templateLimitingMagnitude']) 

1035 

1036 # Create a template with a PSF that is not defined at the image center. 

1037 # First, make an exposure catalog so we can force the template to have 

1038 # a bad (off-image) PSF. It must have a record with a weight field 

1039 # and a BBox in order to let us set the PSF manually. 

1040 template_offimage, _ = makeTestImage() 

1041 schema = afwTable.ExposureTable.makeMinimalSchema() 

1042 weightKey = schema.addField("weight", type="D", doc="Coadd weight") 

1043 exposureCatalog = afwTable.ExposureCatalog(schema) 

1044 record = exposureCatalog.addNew() 

1045 record.setD(weightKey, 1.0) 

1046 record.setBBox(template_offimage.getBBox()) 

1047 kernel = measAlg.DoubleGaussianPsf(7, 7, 2.0).getKernel() 

1048 psf = measAlg.KernelPsf(kernel, template_offimage.getBBox().getCenter()) 

1049 record.setPsf(psf) 

1050 record.setWcs(template_offimage.wcs) 

1051 custom_offimage_psf = CustomCoaddPsf(exposureCatalog, template_offimage.wcs) 

1052 template_offimage.setPsf(custom_offimage_psf) 

1053 

1054 # Test that template PSF size edge cases are handled correctly. 

1055 subtractTask_offimage = self._setup_subtraction() 

1056 _ = subtractTask_offimage.run(template_offimage.clone(), science.clone(), sources) 

1057 # Test that providing no fallbackPsfSize results in a nan template 

1058 # limiting magnitude. 

1059 maglim_template_offimage = subtractTask_offimage._calculateMagLim(template_offimage) 

1060 self.assertTrue(np.isnan(maglim_template_offimage)) 

1061 # Test that the task falls back to the provided fallbackPsfSize, so 

1062 # the template and diffim limiting magnitudes are finite. 

1063 self.assertTrue(np.isfinite(subtractTask_offimage.metadata['templateLimitingMagnitude'])) 

1064 self.assertFloatsAlmostEqual(subtractTask_offimage.metadata['diffimLimitingMagnitude'], 

1065 combineMagLims(subtractTask_offimage.metadata), atol=1e-6) 

1066 

1067 # Test that several other expected metadata metrics exist 

1068 self.assertIn('scienceLimitingMagnitude', subtractTask_good.metadata) 

1069 self.assertIn('templateLimitingMagnitude', subtractTask_good.metadata) 

1070 

1071 # The mean ratio metric should be much worse on the "bad" subtraction. 

1072 self.assertLess(subtractTask_good.metadata['differenceFootprintRatioMean'], 0.02) 

1073 self.assertGreater(subtractTask_bad.metadata['differenceFootprintRatioMean'], 0.12) 

1074 

1075 

1076class AlardLuptonPreconvolveSubtractTest(AlardLuptonSubtractTestBase, lsst.utils.tests.TestCase): 

1077 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask 

1078 

1079 def test_mismatched_template(self): 

1080 """Test that an error is raised if the template 

1081 does not fully contain the science image. 

1082 """ 

1083 xSize = 200 

1084 ySize = 200 

1085 science, sources = makeTestImage(psfSize=self.midPsfSize, xSize=xSize + 20, ySize=ySize + 20) 

1086 template, _ = makeTestImage(psfSize=self.midPsfSize, xSize=xSize, ySize=ySize, 

1087 doApplyCalibration=True) 

1088 task = self._setup_subtraction() 

1089 with self.assertRaises(AssertionError): 

1090 task.run(template, science, sources) 

1091 

1092 def test_equal_images(self): 

1093 """Test that running with enough sources produces reasonable output, 

1094 with the same size psf in the template and science. 

1095 """ 

1096 noiseLevel = 1. 

1097 xSize = 400 

1098 ySize = 400 

1099 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6, 

1100 xSize=xSize, ySize=ySize) 

1101 template, _ = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=7, 

1102 templateBorderSize=20, doApplyCalibration=True, 

1103 xSize=xSize, ySize=ySize) 

1104 task = self._setup_subtraction() 

1105 output = task.run(template, science, sources) 

1106 # There shoud be no NaN values in the Score image 

1107 self.assertTrue(np.all(np.isfinite(output.scoreExposure.image.array))) 

1108 # Mean of Score image should be close to zero. 

1109 meanError = noiseLevel/np.sqrt(output.scoreExposure.image.array.size) 

1110 # Make sure to include pixels with the DETECTED mask bit set. 

1111 statsCtrl = makeStats(badMaskPlanes=("EDGE", "BAD", "NO_DATA")) 

1112 scoreMean = computeRobustStatistics(output.scoreExposure.image, 

1113 output.scoreExposure.mask, 

1114 statsCtrl) 

1115 self.assertFloatsAlmostEqual(scoreMean, 0, atol=5*meanError) 

1116 nea = computePSFNoiseEquivalentArea(science.psf) 

1117 # stddev of Score image should be close to expected value. 

1118 scoreStd = computeRobustStatistics(output.scoreExposure.image, output.scoreExposure.mask, 

1119 statsCtrl=statsCtrl, statistic=afwMath.STDEV) 

1120 self.assertFloatsAlmostEqual(scoreStd, np.sqrt(2)*noiseLevel/np.sqrt(nea), rtol=0.1) 

1121 

1122 def test_incomplete_template_coverage(self): 

1123 noiseLevel = 1. 

1124 border = 20 

1125 xSize = 400 

1126 ySize = 400 

1127 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6, 

1128 xSize=xSize, ySize=ySize) 

1129 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7, 

1130 templateBorderSize=border, doApplyCalibration=True, 

1131 xSize=xSize, ySize=ySize) 

1132 

1133 science_height = science.getBBox().getDimensions().getY() 

1134 

1135 def _run_and_check_coverage(template_coverage, 

1136 requiredTemplateFraction=0.1, 

1137 minTemplateFractionForExpectedSuccess=0.2): 

1138 template_cut = template.clone() 

1139 template_height = int(science_height*template_coverage + border) 

1140 template_cut.image.array[:, template_height:] = 0 

1141 template_cut.mask.array[:, template_height:] = template_cut.mask.getPlaneBitMask('NO_DATA') 

1142 task = self._setup_subtraction( 

1143 requiredTemplateFraction=requiredTemplateFraction, 

1144 minTemplateFractionForExpectedSuccess=minTemplateFractionForExpectedSuccess 

1145 ) 

1146 if template_coverage < task.config.requiredTemplateFraction: 

1147 doRaise = True 

1148 elif template_coverage < task.config.minTemplateFractionForExpectedSuccess: 

1149 doRaise = True 

1150 else: 

1151 doRaise = False 

1152 if doRaise: 

1153 with self.assertRaises(NoWorkFound): 

1154 task.run(template_cut, science.clone(), sources.copy(deep=True)) 

1155 else: 

1156 task.run(template_cut, science.clone(), sources.copy(deep=True)) 

1157 _run_and_check_coverage(template_coverage=0.09) 

1158 _run_and_check_coverage(template_coverage=0.15) 

1159 _run_and_check_coverage(template_coverage=.7) 

1160 

1161 def test_clear_template_mask(self): 

1162 noiseLevel = 1. 

1163 xSize = 400 

1164 ySize = 400 

1165 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6, 

1166 xSize=xSize, ySize=ySize) 

1167 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7, 

1168 templateBorderSize=20, doApplyCalibration=True, 

1169 xSize=xSize, ySize=ySize) 

1170 diffimEmptyMaskPlanes = ["DETECTED", "DETECTED_NEGATIVE"] 

1171 task = self._setup_subtraction() 

1172 # Ensure that each each mask plane is set for some pixels 

1173 mask = template.mask 

1174 x0 = 50 

1175 x1 = 75 

1176 y0 = 150 

1177 y1 = 175 

1178 scienceMaskCheck = {} 

1179 for maskPlane in mask.getMaskPlaneDict().keys(): 

1180 scienceMaskCheck[maskPlane] = np.sum(science.mask.array & mask.getPlaneBitMask(maskPlane) > 0) 

1181 mask.array[x0: x1, y0: y1] |= mask.getPlaneBitMask(maskPlane) 

1182 self.assertTrue(np.sum(mask.array & mask.getPlaneBitMask(maskPlane) > 0)) 

1183 

1184 output = task.run(template, science, sources) 

1185 # Verify that the template mask has been modified in place 

1186 for maskPlane in mask.getMaskPlaneDict().keys(): 

1187 if maskPlane in diffimEmptyMaskPlanes: 

1188 self.assertTrue(np.sum(mask.array & mask.getPlaneBitMask(maskPlane) == 0)) 

1189 elif maskPlane in task.config.preserveTemplateMask: 

1190 self.assertTrue(np.sum(mask.array & mask.getPlaneBitMask(maskPlane) > 0)) 

1191 else: 

1192 self.assertTrue(np.sum(mask.array & mask.getPlaneBitMask(maskPlane) == 0)) 

1193 # Mask planes set in the science image should also be set in the difference 

1194 # Except the "DETECTED" planes should have been cleared 

1195 diffimMask = output.scoreExposure.mask 

1196 for maskPlane, scienceSum in scienceMaskCheck.items(): 

1197 diffimSum = np.sum(diffimMask.array & mask.getPlaneBitMask(maskPlane) > 0) 

1198 if maskPlane in diffimEmptyMaskPlanes: 

1199 self.assertEqual(diffimSum, 0) 

1200 else: 

1201 self.assertTrue(diffimSum >= scienceSum) 

1202 

1203 def test_agnostic_template_psf(self): 

1204 """Test that the Score image is the same whether the template PSF is 

1205 larger or smaller than the science image PSF. 

1206 """ 

1207 noiseLevel = .3 

1208 xSize = 400 

1209 ySize = 400 

1210 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, 

1211 noiseSeed=6, templateBorderSize=0, 

1212 xSize=xSize, ySize=ySize) 

1213 template1, _ = makeTestImage(psfSize=self.badPsfSize, noiseLevel=noiseLevel, 

1214 noiseSeed=7, doApplyCalibration=True, 

1215 xSize=xSize, ySize=ySize) 

1216 template2, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, 

1217 noiseSeed=8, doApplyCalibration=True, 

1218 xSize=xSize, ySize=ySize) 

1219 task = self._setup_subtraction() 

1220 

1221 science_better = task.run(template1, science.clone(), sources) 

1222 template_better = task.run(template2, science, sources) 

1223 bbox = science_better.scoreExposure.getBBox().clippedTo(template_better.scoreExposure.getBBox()) 

1224 

1225 delta = template_better.scoreExposure[bbox].clone() 

1226 delta.image -= science_better.scoreExposure[bbox].image 

1227 delta.variance -= science_better.scoreExposure[bbox].variance 

1228 delta.mask.array &= science_better.scoreExposure[bbox].mask.array 

1229 

1230 statsCtrl = makeStats(badMaskPlanes=("EDGE", "BAD", "NO_DATA")) 

1231 # Mean of delta should be very close to zero. 

1232 nGoodPix = np.sum(np.isfinite(delta.image.array)) 

1233 meanError = 2*noiseLevel/np.sqrt(nGoodPix) 

1234 deltaMean = computeRobustStatistics(delta.image, delta.mask, statsCtrl) 

1235 deltaStd = computeRobustStatistics(delta.image, delta.mask, statsCtrl, 

1236 statistic=afwMath.STDEV) 

1237 self.assertFloatsAlmostEqual(deltaMean, 0, atol=5*meanError) 

1238 nea = computePSFNoiseEquivalentArea(science.psf) 

1239 # stddev of Score image should be close to expected value 

1240 self.assertFloatsAlmostEqual(deltaStd, np.sqrt(2)*noiseLevel/np.sqrt(nea), rtol=.1) 

1241 

1242 def test_few_sources(self): 

1243 """Test with only 1 source, to check that we get a useful error. 

1244 """ 

1245 xSize = 256 

1246 ySize = 256 

1247 science, sources = makeTestImage(psfSize=self.midPsfSize, nSrc=10, xSize=xSize, ySize=ySize) 

1248 template, _ = makeTestImage(psfSize=self.goodPsfSize, nSrc=10, xSize=xSize, ySize=ySize, 

1249 doApplyCalibration=True) 

1250 task = self._setup_subtraction() 

1251 sources = sources[0:1] 

1252 with self.assertRaises(InsufficientKernelSourcesError): 

1253 task.run(template, science, sources) 

1254 

1255 def test_background_subtraction(self): 

1256 """Check that we can recover the background, 

1257 and that it is subtracted correctly in the Score image. 

1258 """ 

1259 noiseLevel = 1. 

1260 xSize = 512 

1261 ySize = 512 

1262 x0 = 123 

1263 y0 = 456 

1264 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7, 

1265 templateBorderSize=20, 

1266 xSize=xSize, ySize=ySize, x0=x0, y0=y0, 

1267 doApplyCalibration=True) 

1268 params = [2.2, 2.1, 2.0, 1.2, 1.1, 1.0] 

1269 

1270 bbox2D = lsst.geom.Box2D(lsst.geom.Point2D(x0, y0), lsst.geom.Extent2D(xSize, ySize)) 

1271 background_model = afwMath.Chebyshev1Function2D(params, bbox2D) 

1272 science, sources = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=6, 

1273 background=background_model, 

1274 xSize=xSize, ySize=ySize, x0=x0, y0=y0) 

1275 # Don't use ``self._setup_subtraction()`` here. 

1276 # Modifying the config of a subtask is messy. 

1277 config = subtractImages.AlardLuptonPreconvolveSubtractTask.ConfigClass() 

1278 

1279 config.sourceSelector.signalToNoise.fluxField = "truth_instFlux" 

1280 config.sourceSelector.signalToNoise.errField = "truth_instFluxErr" 

1281 config.doSubtractBackground = True 

1282 

1283 config.makeKernel.kernel.name = "AL" 

1284 config.makeKernel.kernel.active.fitForBackground = True 

1285 config.makeKernel.kernel.active.spatialKernelOrder = 1 

1286 config.makeKernel.kernel.active.spatialBgOrder = 2 

1287 statsCtrl = makeStats(badMaskPlanes=("EDGE", "BAD", "NO_DATA")) 

1288 

1289 task = subtractImages.AlardLuptonPreconvolveSubtractTask(config=config) 

1290 output = task.run(template.clone(), science.clone(), sources) 

1291 

1292 # We should be fitting the same number of parameters as were in the input model 

1293 self.assertEqual(output.backgroundModel.getNParameters(), background_model.getNParameters()) 

1294 

1295 # The parameters of the background fit should be close to the input model 

1296 self.assertFloatsAlmostEqual(np.array(output.backgroundModel.getParameters()), 

1297 np.array(params), rtol=0.2) 

1298 

1299 # stddev of Score image should be close to expected value. 

1300 # This will fail if we have mis-subtracted the background. 

1301 stdVal = computeRobustStatistics(output.scoreExposure.image, output.scoreExposure.mask, 

1302 statsCtrl, statistic=afwMath.STDEV) 

1303 # get the img psf Noise Equivalent Area value 

1304 nea = computePSFNoiseEquivalentArea(science.psf) 

1305 self.assertFloatsAlmostEqual(stdVal, np.sqrt(2)*noiseLevel/np.sqrt(nea), rtol=0.12) 

1306 

1307 def test_scale_variance(self): 

1308 """Check variance scaling of the Score image. 

1309 """ 

1310 scienceNoiseLevel = 4. 

1311 templateNoiseLevel = 2. 

1312 scaleFactor = 1.345 

1313 xSize = 400 

1314 ySize = 400 

1315 # Make sure to include pixels with the DETECTED mask bit set. 

1316 statsCtrl = makeStats(badMaskPlanes=("EDGE", "BAD", "NO_DATA")) 

1317 

1318 def _run_and_check_images(science, template, sources, statsCtrl, 

1319 doDecorrelation, doScaleVariance, scaleFactor=1.): 

1320 """Check that the variance plane matches the expected value for 

1321 different configurations of ``doDecorrelation`` and ``doScaleVariance``. 

1322 """ 

1323 

1324 task = self._setup_subtraction(doDecorrelation=doDecorrelation, 

1325 doScaleVariance=doScaleVariance, 

1326 ) 

1327 output = task.run(template.clone(), science.clone(), sources) 

1328 if doScaleVariance: 

1329 self.assertFloatsAlmostEqual(task.metadata["scaleScienceVarianceFactor"], 

1330 scaleFactor, atol=0.05) 

1331 

1332 scienceNoise = computeRobustStatistics(science.variance, science.mask, statsCtrl) 

1333 # get the img psf Noise Equivalent Area value 

1334 nea = computePSFNoiseEquivalentArea(science.psf) 

1335 scienceNoise /= nea 

1336 if doDecorrelation: 

1337 templateNoise = computeRobustStatistics(template.variance, template.mask, statsCtrl) 

1338 templateNoise /= nea 

1339 else: 

1340 # Don't divide by NEA in this case, since the template is convolved 

1341 # and in the same units as the Score exposure. 

1342 templateNoise = computeRobustStatistics(output.matchedTemplate.variance, 

1343 output.matchedTemplate.mask, 

1344 statsCtrl) 

1345 if doScaleVariance: 

1346 # Only the science variance is scaled here. The template 

1347 # variance is scaled independently in ``GetTemplateTask``. 

1348 scienceNoise *= scaleFactor 

1349 varMean = computeRobustStatistics(output.scoreExposure.variance, 

1350 output.scoreExposure.mask, 

1351 statsCtrl) 

1352 self.assertFloatsAlmostEqual(varMean, scienceNoise + templateNoise, rtol=0.1) 

1353 

1354 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=scienceNoiseLevel, noiseSeed=6, 

1355 xSize=xSize, ySize=ySize) 

1356 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=templateNoiseLevel, noiseSeed=7, 

1357 templateBorderSize=20, doApplyCalibration=True, 

1358 xSize=xSize, ySize=ySize) 

1359 # Verify that the variance plane of the Score image is correct 

1360 # when the template and science variance planes are correct 

1361 _run_and_check_images(science, template, sources, statsCtrl, 

1362 doDecorrelation=True, doScaleVariance=True) 

1363 _run_and_check_images(science, template, sources, statsCtrl, 

1364 doDecorrelation=True, doScaleVariance=False) 

1365 _run_and_check_images(science, template, sources, statsCtrl, 

1366 doDecorrelation=False, doScaleVariance=True) 

1367 _run_and_check_images(science, template, sources, statsCtrl, 

1368 doDecorrelation=False, doScaleVariance=False) 

1369 

1370 # Verify that the variance plane of the Score image is correct 

1371 # when the input science variance plane is incorrect 

1372 science.variance.array /= scaleFactor 

1373 _run_and_check_images(science, template, sources, statsCtrl, 

1374 doDecorrelation=True, doScaleVariance=True, scaleFactor=scaleFactor) 

1375 _run_and_check_images(science, template, sources, statsCtrl, 

1376 doDecorrelation=True, doScaleVariance=False, scaleFactor=scaleFactor) 

1377 _run_and_check_images(science, template, sources, statsCtrl, 

1378 doDecorrelation=False, doScaleVariance=True, scaleFactor=scaleFactor) 

1379 _run_and_check_images(science, template, sources, statsCtrl, 

1380 doDecorrelation=False, doScaleVariance=False, scaleFactor=scaleFactor) 

1381 

1382 def test_exposure_properties(self): 

1383 """Check that all necessary exposure metadata is included 

1384 with the Score image. 

1385 """ 

1386 noiseLevel = 1. 

1387 xSize = 400 

1388 ySize = 400 

1389 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6, 

1390 xSize=xSize, ySize=ySize) 

1391 psf = science.psf 

1392 psfAvgPos = psf.getAveragePosition() 

1393 psfSize = getPsfFwhm(science.psf) 

1394 psfImg = psf.computeKernelImage(psfAvgPos) 

1395 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7, 

1396 templateBorderSize=20, doApplyCalibration=True, 

1397 xSize=xSize, ySize=ySize) 

1398 

1399 def _run_and_check_images(doDecorrelation): 

1400 """Check that the metadata is correct with or without decorrelation. 

1401 """ 

1402 task = self._setup_subtraction(doDecorrelation=doDecorrelation) 

1403 output = task.run(template.clone(), science.clone(), sources) 

1404 psfOut = output.scoreExposure.psf 

1405 psfAvgPos = psfOut.getAveragePosition() 

1406 if doDecorrelation: 

1407 # Decorrelation requires recalculating the PSF, 

1408 # so it will not be the same as the input 

1409 psfOutSize = getPsfFwhm(science.psf) 

1410 self.assertFloatsAlmostEqual(psfSize, psfOutSize) 

1411 else: 

1412 psfOutImg = psfOut.computeKernelImage(psfAvgPos) 

1413 self.assertImagesAlmostEqual(psfImg, psfOutImg) 

1414 

1415 # check PSF, WCS, bbox, filterLabel, photoCalib 

1416 self.assertWcsAlmostEqualOverBBox(science.wcs, output.scoreExposure.wcs, science.getBBox()) 

1417 self.assertEqual(science.filter, output.scoreExposure.filter) 

1418 self.assertEqual(science.photoCalib, output.scoreExposure.photoCalib) 

1419 _run_and_check_images(doDecorrelation=True) 

1420 _run_and_check_images(doDecorrelation=False) 

1421 

1422 

1423class SimplifiedSubtractTest(AlardLuptonSubtractTestBase, lsst.utils.tests.TestCase): 

1424 subtractTask = subtractImages.SimplifiedSubtractTask 

1425 

1426 def test_runSimplifiedTaskWithExistingKernel(self): 

1427 """Test that the simplified task produces the same output as 

1428 `AlardLuptonSubtractTask` if it uses the AL kernel. 

1429 """ 

1430 noiseLevel = 1. 

1431 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6) 

1432 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7, 

1433 templateBorderSize=20, doApplyCalibration=True) 

1434 alTask = AlardLuptonSubtractTest._setup_subtraction(AlardLuptonSubtractTest()) 

1435 task = self._setup_subtraction(useExistingKernel=True) 

1436 

1437 alResults = alTask.run(template.clone(), science.clone(), sources) 

1438 results = task.run(template.clone(), science.clone(), 

1439 inputPsfMatchingKernel=alResults.psfMatchingKernel) 

1440 

1441 self.assertMaskedImagesEqual(alResults.difference, results.difference) 

1442 

1443 def test_runSimplifiedTaskWithSourceDetection(self): 

1444 """Test that the simplified task with source detection produces 

1445 reasonable output. 

1446 """ 

1447 noiseLevel = 1. 

1448 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6) 

1449 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7, 

1450 templateBorderSize=20, doApplyCalibration=True) 

1451 task = self._setup_subtraction(useExistingKernel=False, 

1452 fluxField="base_PsfFlux_instFlux", 

1453 errField="base_PsfFlux_instFluxErr", 

1454 ) 

1455 

1456 output = task.run(template, science) 

1457 

1458 # There shoud be no NaN values in the difference image 

1459 self.assertTrue(np.all(np.isfinite(output.difference.image.array))) 

1460 # Mean of difference image should be close to zero. 

1461 meanError = noiseLevel/np.sqrt(output.difference.image.array.size) 

1462 # Make sure to include pixels with the DETECTED mask bit set. 

1463 statsCtrl = makeStats(badMaskPlanes=("EDGE", "BAD", "NO_DATA", "DETECTED", "DETECTED_NEGATIVE")) 

1464 differenceMean = computeRobustStatistics(output.difference.image, output.difference.mask, statsCtrl) 

1465 self.assertFloatsAlmostEqual(differenceMean, 0, atol=5*meanError) 

1466 # stddev of difference image should be close to expected value. 

1467 differenceStd = computeRobustStatistics(output.difference.image, output.difference.mask, 

1468 makeStats(), statistic=afwMath.STDEV) 

1469 self.assertFloatsAlmostEqual(differenceStd, np.sqrt(2)*noiseLevel, rtol=0.1) 

1470 

1471 

1472class OutputImageTypeTest(lsst.utils.tests.TestCase): 

1473 """Tests of the optional conversion of the image outputs to 

1474 `lsst.images` types. 

1475 """ 

1476 def setUp(self): 

1477 # `makeTestImage` does not attach a detector or a VisitInfo, so 

1478 # supply a detector with real amplifier geometry, a data id, and the 

1479 # exposure record the conversion takes its metadata from. 

1480 self.detector = list(lsst.afw.cameraGeom.testUtils.CameraWrapper().camera)[0] 

1481 self.dataId = DataCoordinate.standardize( 

1482 instrument="testCam", 

1483 visit=12345, 

1484 detector=self.detector.getId(), 

1485 universe=DimensionUniverse(), 

1486 ) 

1487 self.exposureRecord = makeTestExposureRecord(self.dataId.universe, 

1488 instrument=self.dataId["instrument"], 

1489 visit=self.dataId["visit"]) 

1490 

1491 def _setup_subtraction(self, taskClass, **kwargs): 

1492 """Configure one of the subtraction tasks with the shared test 

1493 settings. 

1494 

1495 Parameters 

1496 ---------- 

1497 taskClass : `lsst.pipe.base.PipelineTask` 

1498 The subtraction task class to configure. 

1499 **kwargs 

1500 Any additional config parameters to set. 

1501 

1502 Returns 

1503 ------- 

1504 task : `lsst.pipe.base.PipelineTask` 

1505 The configured task. 

1506 """ 

1507 helper = AlardLuptonSubtractTestBase() 

1508 helper.subtractTask = taskClass 

1509 return helper._setup_subtraction(**kwargs) 

1510 

1511 def _make_images(self): 

1512 """Make a science and template image pair to subtract. 

1513 

1514 Returns 

1515 ------- 

1516 science : `lsst.afw.image.ExposureF` 

1517 The science image, with a detector attached. 

1518 template : `lsst.afw.image.ExposureF` 

1519 The template image, with no detector attached. 

1520 sources : `lsst.afw.table.SourceCatalog` 

1521 Sources detected on the science image. 

1522 """ 

1523 noiseLevel = 1. 

1524 science, sources = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel, noiseSeed=6) 

1525 template, _ = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel, noiseSeed=7, 

1526 templateBorderSize=20, doApplyCalibration=True) 

1527 science.setDetector(self.detector) 

1528 for exposure in (science, template): 

1529 exposure.info.setVisitInfo(makeTestVisitInfo(self.dataId["visit"])) 

1530 return science, template, sources 

1531 

1532 def _check_converted_kernel(self, kernel, legacyKernel, bbox): 

1533 """Check that a converted kernel evaluates the same as the legacy 

1534 kernel it came from. 

1535 

1536 Parameters 

1537 ---------- 

1538 kernel : `lsst.images.convolution_kernels.ImageBasisConvolutionKernel` 

1539 The converted kernel to check. 

1540 legacyKernel : `lsst.afw.math.LinearCombinationKernel` 

1541 The kernel before conversion. 

1542 bbox : `lsst.geom.Box2I` 

1543 Region of the image the kernel was fit on. 

1544 """ 

1545 self.assertIsInstance(kernel, ImageBasisConvolutionKernel) 

1546 roundTrip = kernel.to_legacy() 

1547 for point in (bbox.getMin(), bbox.getCenter(), bbox.getMax()): 

1548 expected = afwImage.ImageD(legacyKernel.getDimensions()) 

1549 legacyKernel.computeImage(expected, False, point.x, point.y) 

1550 actual = afwImage.ImageD(roundTrip.getDimensions()) 

1551 roundTrip.computeImage(actual, False, point.x, point.y) 

1552 np.testing.assert_allclose(actual.array, expected.array, rtol=0, atol=1e-12) 

1553 

1554 def _check_converted(self, image, kernelExpected=True, photometricScaling=None): 

1555 """Check that one converted image is a fully-formed DifferenceImage. 

1556 

1557 Parameters 

1558 ---------- 

1559 image : `lsst.images.DifferenceImage` 

1560 The converted output to check. 

1561 kernelExpected : `bool`, optional 

1562 Should the PSF matching kernel be attached? 

1563 photometricScaling : `lsst.images.fields.BaseField`, optional 

1564 Photometric scaling the image should carry, or `None` if it should 

1565 have none. 

1566 """ 

1567 self.assertIsInstance(image, DifferenceImage) 

1568 self.assertEqual(image.unit, u.nJy) 

1569 # The instrument and visit come from the exposure record passed to 

1570 # the conversion. 

1571 frame = image.sky_projection.pixel_frame 

1572 self.assertEqual(frame.instrument, self.dataId["instrument"]) 

1573 self.assertEqual(frame.visit, self.dataId["visit"]) 

1574 self.assertEqual(image.detector.id, self.detector.getId()) 

1575 if kernelExpected: 1575 ↛ 1579line 1575 didn't jump to line 1579 because the condition on line 1575 was always true

1576 self.assertIsNotNone(image.kernel) 

1577 else: 

1578 # An unset kernel raises, rather than returning None. 

1579 with self.assertRaises(AttributeError): 

1580 image.kernel 

1581 if photometricScaling is None: 

1582 self.assertIsNone(image.photometric_scaling) 

1583 else: 

1584 self.assertIs(image.photometric_scaling, photometricScaling) 

1585 

1586 def test_connections_legacy(self): 

1587 """Test that the default configuration leaves the output connections 

1588 with the legacy storage classes. 

1589 """ 

1590 for configClass, connectionsClass in ( 

1591 (subtractImages.AlardLuptonSubtractConfig, 

1592 subtractImages.AlardLuptonSubtractConnections), 

1593 (subtractImages.AlardLuptonPreconvolveSubtractConfig, 

1594 subtractImages.AlardLuptonPreconvolveSubtractConnections), 

1595 (subtractImages.SimplifiedSubtractConfig, 

1596 subtractImages.SimplifiedSubtractConnections), 

1597 ): 

1598 with self.subTest(configClass=configClass.__name__): 

1599 config = configClass() 

1600 self.assertEqual(config.image_type, "legacy") 

1601 connections = connectionsClass(config=config) 

1602 for name in connections.outputs: 

1603 if name in ("difference", "matchedTemplate", "scoreExposure"): 

1604 self.assertEqual(getattr(connections, name).storageClass, "ExposureF") 

1605 self.assertEqual(connections.science.storageClass, "ExposureF") 

1606 for name in ("psfMatchingKernel", "inputPsfMatchingKernel"): 

1607 if name in connections.outputs or name in connections.inputs: 

1608 self.assertEqual(getattr(connections, name).storageClass, "MatchingKernel") 

1609 

1610 def test_connections_future(self): 

1611 """Test that ``image_type="future"`` swaps the storage class of 

1612 every image and kernel output. 

1613 """ 

1614 for configClass, connectionsClass, expected in ( 

1615 (subtractImages.AlardLuptonSubtractConfig, 

1616 subtractImages.AlardLuptonSubtractConnections, 

1617 ("difference", "matchedTemplate")), 

1618 (subtractImages.AlardLuptonPreconvolveSubtractConfig, 

1619 subtractImages.AlardLuptonPreconvolveSubtractConnections, 

1620 ("scoreExposure",)), 

1621 (subtractImages.SimplifiedSubtractConfig, 

1622 subtractImages.SimplifiedSubtractConnections, 

1623 ("difference", "matchedTemplate")), 

1624 ): 

1625 with self.subTest(configClass=configClass.__name__): 

1626 config = configClass() 

1627 config.image_type = "future" 

1628 connections = connectionsClass(config=config) 

1629 for name in expected: 

1630 self.assertEqual(getattr(connections, name).storageClass, "DifferenceImage") 

1631 # The science image is read in the new type to get its 

1632 # photometric scaling. 

1633 self.assertEqual(connections.science.storageClass, "VisitImage") 

1634 if "psfMatchingKernel" in connections.outputs: 

1635 self.assertEqual(connections.psfMatchingKernel.storageClass, "ConvolutionKernel") 

1636 

1637 def test_connections_future_simplified_existing_kernel(self): 

1638 """Test the storage class swap for the connection variant that 

1639 consumes an existing kernel. 

1640 """ 

1641 config = subtractImages.SimplifiedSubtractConfig() 

1642 config.image_type = "future" 

1643 config.useExistingKernel = True 

1644 connections = subtractImages.SimplifiedSubtractConnections(config=config) 

1645 self.assertEqual(connections.difference.storageClass, "DifferenceImage") 

1646 self.assertEqual(connections.matchedTemplate.storageClass, "DifferenceImage") 

1647 self.assertEqual(connections.inputPsfMatchingKernel.storageClass, "ConvolutionKernel") 

1648 self.assertNotIn("psfMatchingKernel", connections.outputs) 

1649 

1650 def test_convert_outputs_to_future(self): 

1651 """Test converting the outputs of `AlardLuptonSubtractTask`. 

1652 """ 

1653 science, template, sources = self._make_images() 

1654 task = self._setup_subtraction(subtractImages.AlardLuptonSubtractTask, 

1655 image_type="future") 

1656 results = task.run(template, science, sources) 

1657 # The template has no detector, so neither does the matched template. 

1658 self.assertIsNone(results.matchedTemplate.getDetector()) 

1659 legacy = {name: getattr(results, name).clone() for name in ("difference", "matchedTemplate")} 

1660 legacyKernel = results.psfMatchingKernel 

1661 

1662 task.convert_outputs_to_future(results, self.exposureRecord) 

1663 

1664 for name, exposure in legacy.items(): 

1665 with self.subTest(name=name): 

1666 self._check_converted(getattr(results, name)) 

1667 compare_masked_image_to_legacy(getattr(results, name), exposure.maskedImage, 

1668 plane_map=get_legacy_difference_image_mask_planes()) 

1669 self.assertEqual(results.difference.bbox.to_legacy(), science.getBBox()) 

1670 self._check_converted_kernel(results.psfMatchingKernel, legacyKernel, science.getBBox()) 

1671 self.assertIsInstance(results.kernelSources, afwTable.SourceCatalog) 

1672 

1673 def test_convert_outputs_to_future_score(self): 

1674 """Test converting the outputs of 

1675 `AlardLuptonPreconvolveSubtractTask`, which writes a score image. 

1676 """ 

1677 science, template, sources = self._make_images() 

1678 task = self._setup_subtraction(subtractImages.AlardLuptonPreconvolveSubtractTask, 

1679 image_type="future") 

1680 results = task.run(template, science, sources) 

1681 legacyKernel = results.psfMatchingKernel 

1682 

1683 task.convert_outputs_to_future(results, self.exposureRecord) 

1684 

1685 self._check_converted(results.scoreExposure) 

1686 self._check_converted_kernel(results.psfMatchingKernel, legacyKernel, science.getBBox()) 

1687 

1688 def test_convert_outputs_to_future_simplified(self): 

1689 """Test converting the outputs of `SimplifiedSubtractTask` running 

1690 with a pre-existing kernel. 

1691 """ 

1692 science, template, sources = self._make_images() 

1693 alTask = self._setup_subtraction(subtractImages.AlardLuptonSubtractTask) 

1694 alResults = alTask.run(template.clone(), science.clone(), sources) 

1695 task = self._setup_subtraction(subtractImages.SimplifiedSubtractTask, 

1696 useExistingKernel=True, 

1697 image_type="future") 

1698 results = task.run(template.clone(), science.clone(), 

1699 inputPsfMatchingKernel=alResults.psfMatchingKernel) 

1700 

1701 task.convert_outputs_to_future(results, self.exposureRecord) 

1702 

1703 self._check_converted(results.difference) 

1704 self._check_converted(results.matchedTemplate) 

1705 

1706 def test_convert_outputs_to_future_photometric_scaling(self): 

1707 """Test that every converted output carries the photometric scaling 

1708 of the science image. 

1709 """ 

1710 science, template, sources = self._make_images() 

1711 photometricScaling = self._makePhotometricScaling(science) 

1712 for taskClass, names in ((subtractImages.AlardLuptonSubtractTask, ("difference", "matchedTemplate")), 

1713 (subtractImages.AlardLuptonPreconvolveSubtractTask, ("scoreExposure",))): 

1714 with self.subTest(taskClass=taskClass.__name__): 

1715 task = self._setup_subtraction(taskClass, image_type="future") 

1716 results = task.run(template.clone(), science.clone(), sources) 

1717 

1718 task.convert_outputs_to_future(results, self.exposureRecord, 

1719 photometricScaling=photometricScaling) 

1720 

1721 for name in names: 

1722 self._check_converted(getattr(results, name), photometricScaling=photometricScaling) 

1723 

1724 def test_convert_outputs_to_future_missing_output(self): 

1725 """Test that conversion skips outputs that are absent from the result 

1726 struct, instead of failing. 

1727 """ 

1728 science, template, sources = self._make_images() 

1729 task = self._setup_subtraction(subtractImages.AlardLuptonSubtractTask, 

1730 image_type="future") 

1731 results = task.run(template, science, sources) 

1732 del results.matchedTemplate 

1733 

1734 task.convert_outputs_to_future(results, self.exposureRecord) 

1735 

1736 self._check_converted(results.difference) 

1737 self.assertFalse(hasattr(results, "matchedTemplate")) 

1738 

1739 @staticmethod 

1740 def _makeTemplateInfo(tract=9813, patch=42): 

1741 """Return one record of a coadd that went into a template. 

1742 

1743 Parameters 

1744 ---------- 

1745 tract, patch : `int`, optional 

1746 Identifiers of the coadd the record describes. 

1747 

1748 Returns 

1749 ------- 

1750 info : `lsst.images.DifferenceImageTemplateInfo` 

1751 The record. 

1752 """ 

1753 return DifferenceImageTemplateInfo( 

1754 skymap="skymap", tract=tract, patch=patch, dataset_id=uuid.uuid4(), 

1755 dataset_run="a/template/run", bounds=Polygon.from_box(Box.factory[0:50, 0:50]), 

1756 psf_shape_xx=4.0, psf_shape_yy=4.0, psf_shape_xy=0.0, psf_shape_flag=False) 

1757 

1758 @staticmethod 

1759 def _makePhotometricScaling(science): 

1760 """Return a photometric scaling that is easy to tell apart from the 

1761 one converted from the science image's own calibration. 

1762 

1763 Parameters 

1764 ---------- 

1765 science : `lsst.afw.image.ExposureF` 

1766 Science image from `_make_images`. 

1767 

1768 Returns 

1769 ------- 

1770 photometricScaling : `lsst.images.fields.ChebyshevField` 

1771 A constant scaling over the science image. 

1772 """ 

1773 return ChebyshevField(Box.from_legacy(science.getBBox()), np.array([[0.25]]), 

1774 unit=u.nJy/u.electron) 

1775 

1776 def _make_future_science(self, science, photometricScaling=None): 

1777 """Convert a legacy science image to the type the ``future`` mode 

1778 reads. 

1779 

1780 Parameters 

1781 ---------- 

1782 science : `lsst.afw.image.ExposureF` 

1783 Science image from `_make_images`. 

1784 photometricScaling : `lsst.images.fields.BaseField`, optional 

1785 Photometric scaling to attach in place of the converted one. 

1786 

1787 Returns 

1788 ------- 

1789 science : `lsst.images.VisitImage` 

1790 The converted science image. 

1791 """ 

1792 image = VisitImage.from_legacy(science.clone(), unit=u.nJy, exposure_record=self.exposureRecord) 

1793 if photometricScaling is not None: 

1794 image.photometric_scaling = photometricScaling 

1795 return image 

1796 

1797 def _make_future_template(self, template, templateInfo=None): 

1798 """Convert a legacy template to the type the ``future`` mode reads. 

1799 

1800 Parameters 

1801 ---------- 

1802 template : `lsst.afw.image.ExposureF` 

1803 Template from `_make_images`, which has no detector. 

1804 templateInfo : `list` [`lsst.images.DifferenceImageTemplateInfo`], 

1805 optional 

1806 Record of the coadds that went into it, left off if `None`. 

1807 

1808 Returns 

1809 ------- 

1810 template : `lsst.images.DifferenceImage` 

1811 The converted template. 

1812 """ 

1813 template = template.clone() 

1814 template.setDetector(self.detector) 

1815 image = DifferenceImage.from_legacy(template, unit=u.nJy, exposure_record=self.exposureRecord) 

1816 if templateInfo is not None: 

1817 image.templates = templateInfo 

1818 return image 

1819 

1820 def test_convert_outputs_to_future_templates(self): 

1821 """The record of the coadds that went into the template is copied to 

1822 every converted output. 

1823 """ 

1824 science, template, sources = self._make_images() 

1825 task = self._setup_subtraction(subtractImages.AlardLuptonSubtractTask, 

1826 image_type="future") 

1827 results = task.run(template, science, sources) 

1828 templateInfo = [self._makeTemplateInfo(patch=41), self._makeTemplateInfo(patch=42)] 

1829 

1830 task.convert_outputs_to_future(results, self.exposureRecord, templateInfo=templateInfo) 

1831 

1832 for image in (results.difference, results.matchedTemplate): 

1833 self.assertEqual([(t.tract, t.patch) for t in image.templates], 

1834 [(9813, 41), (9813, 42)]) 

1835 

1836 def test_convert_outputs_to_future_no_templates(self): 

1837 """A template with no record of its coadds leaves the outputs without 

1838 one, rather than with an empty list. 

1839 """ 

1840 science, template, sources = self._make_images() 

1841 task = self._setup_subtraction(subtractImages.AlardLuptonSubtractTask, 

1842 image_type="future") 

1843 results = task.run(template, science, sources) 

1844 

1845 task.convert_outputs_to_future(results, self.exposureRecord) 

1846 

1847 with self.assertRaises(AttributeError): 

1848 results.difference.templates 

1849 

1850 def test_run_quantum_puts_future_types(self): 

1851 """``runQuantum`` converts before putting in future mode, and does not 

1852 convert in legacy mode. 

1853 """ 

1854 science, template, sources = self._make_images() 

1855 templateInfo = [self._makeTemplateInfo()] 

1856 photometricScaling = self._makePhotometricScaling(science) 

1857 for imageType, expectedType in (("legacy", afwImage.ExposureF), 

1858 ("future", DifferenceImage)): 

1859 with self.subTest(imageType=imageType): 

1860 task = self._setup_subtraction(subtractImages.AlardLuptonSubtractTask, 

1861 image_type=imageType) 

1862 if imageType == "future": 

1863 inputTemplate = self._make_future_template(template, templateInfo) 

1864 inputScience = self._make_future_science(science, photometricScaling) 

1865 else: 

1866 inputTemplate = template.clone() 

1867 inputScience = science.clone() 

1868 butlerQC = _RecordingQuantumContext(self.dataId) 

1869 task.runQuantum(butlerQC, 

1870 _FakeRefs(template=inputTemplate, science=inputScience, 

1871 sources=sources), 

1872 _FakeRefs()) 

1873 self.assertIsInstance(butlerQC.put_values.difference, expectedType) 

1874 self.assertIsInstance(butlerQC.put_values.matchedTemplate, expectedType) 

1875 # The kernel is a separate output in both modes. 

1876 self.assertIsInstance(butlerQC.put_values.psfMatchingKernel, 

1877 ImageBasisConvolutionKernel if imageType == "future" 

1878 else afwMath.LinearCombinationKernel) 

1879 if imageType == "future": 

1880 for name in ("difference", "matchedTemplate"): 

1881 image = getattr(butlerQC.put_values, name) 

1882 self.assertEqual(image.templates, templateInfo) 

1883 self.assertIs(image.photometric_scaling, photometricScaling) 

1884 

1885 def test_run_quantum_template_without_provenance(self): 

1886 """A template with no record of its coadds raises an error that 

1887 names what is missing, and writes nothing. 

1888 """ 

1889 science, template, sources = self._make_images() 

1890 task = self._setup_subtraction(subtractImages.AlardLuptonSubtractTask, image_type="future") 

1891 butlerQC = _RecordingQuantumContext(self.dataId) 

1892 refs = _FakeRefs(template=self._make_future_template(template), 

1893 science=self._make_future_science(science), sources=sources) 

1894 

1895 with self.assertRaisesRegex(InvalidQuantumError, "'template'.*template coadd records"): 

1896 task.runQuantum(butlerQC, refs, _FakeRefs()) 

1897 self.assertIsNone(butlerQC.put_values) 

1898 

1899 def test_run_quantum_simplified_puts_future_types(self): 

1900 """``runQuantum`` converts the outputs of `SimplifiedSubtractTask`, 

1901 which loads an existing kernel rather than fitting one. 

1902 """ 

1903 science, template, sources = self._make_images() 

1904 templateInfo = [self._makeTemplateInfo()] 

1905 photometricScaling = self._makePhotometricScaling(science) 

1906 alTask = self._setup_subtraction(subtractImages.AlardLuptonSubtractTask) 

1907 alResults = alTask.run(template.clone(), science.clone(), sources) 

1908 for useExistingKernel in (True, False): 

1909 with self.subTest(useExistingKernel=useExistingKernel): 

1910 task = self._setup_subtraction(subtractImages.SimplifiedSubtractTask, 

1911 useExistingKernel=useExistingKernel, 

1912 image_type="future") 

1913 refs = {"template": self._make_future_template(template, templateInfo), 

1914 "science": self._make_future_science(science, photometricScaling)} 

1915 if useExistingKernel: 

1916 refs["inputPsfMatchingKernel"] = ImageBasisConvolutionKernel.from_legacy( 

1917 alResults.psfMatchingKernel 

1918 ) 

1919 butlerQC = _RecordingQuantumContext(self.dataId) 

1920 

1921 task.runQuantum(butlerQC, _FakeRefs(**refs), _FakeRefs()) 

1922 

1923 self._check_converted(butlerQC.put_values.difference, 

1924 photometricScaling=photometricScaling) 

1925 self._check_converted(butlerQC.put_values.matchedTemplate, 

1926 photometricScaling=photometricScaling) 

1927 self.assertEqual(butlerQC.put_values.difference.templates, templateInfo) 

1928 self.assertIsInstance(butlerQC.put_values.psfMatchingKernel, ImageBasisConvolutionKernel) 

1929 if useExistingKernel: 

1930 # The input kernel is the one the subtraction used. 

1931 self._check_converted_kernel(butlerQC.put_values.difference.kernel, 

1932 alResults.psfMatchingKernel, science.getBBox()) 

1933 

1934 

1935class _FakeRefs: 

1936 """Stand-in for an input or output ``QuantizedConnection``.""" 

1937 def __init__(self, **kwargs): 

1938 self.__dict__.update(kwargs) 

1939 

1940 

1941class _RecordingQuantumContext: 

1942 """Minimal `~lsst.pipe.base.QuantumContext` that records what was put. 

1943 

1944 Parameters 

1945 ---------- 

1946 dataId : `dict` 

1947 Data ID to report as the quantum data ID. 

1948 """ 

1949 def __init__(self, dataId): 

1950 self.quantum = _FakeRefs(dataId=dataId) 

1951 self.put_values = None 

1952 

1953 def get(self, refs): 

1954 return dict(refs.__dict__) 

1955 

1956 def put(self, values, refs): 

1957 self.put_values = values 

1958 

1959 

1960def setup_module(module): 

1961 lsst.utils.tests.init() 

1962 

1963 

1964class MemoryTestCase(lsst.utils.tests.MemoryTestCase): 

1965 pass 

1966 

1967 

1968if __name__ == "__main__": 1968 ↛ 1969line 1968 didn't jump to line 1969 because the condition on line 1968 was never true

1969 lsst.utils.tests.init() 

1970 unittest.main()