Coverage for tests/test_getTemplate.py: 89%

245 statements  

« prev     ^ index     » next       coverage.py v7.16.2, created at 2026-09-29 10:36 +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 collections 

23import itertools 

24import unittest 

25 

26import numpy as np 

27 

28import lsst.afw.geom 

29import lsst.afw.image 

30import lsst.afw.math 

31import lsst.afw.table 

32from lsst.daf.butler import DataCoordinate, DimensionUniverse 

33import lsst.geom 

34import lsst.ip.diffim 

35import lsst.meas.algorithms 

36import lsst.meas.base.tests 

37import lsst.pipe.base as pipeBase 

38import lsst.skymap 

39import lsst.utils.tests 

40 

41from utils import generate_data_id 

42 

43# Change this to True, `setup display_ds9`, and open ds9 (or use another afw 

44# display backend) to show the tract/patch layouts on the image. 

45debug = False 

46if debug: 46 ↛ 47line 46 didn't jump to line 47 because the condition on line 46 was never true

47 import lsst.afw.display 

48 display = lsst.afw.display.Display() 

49 display.frame = 1 

50 

51 

52def _showTemplate(box, template): 

53 """Show the corners of the template we made in this test.""" 

54 for point in box.getCorners(): 

55 display.dot("+", point.x, point.y, ctype="orange", size=40) 

56 display.frame = 2 

57 display.image(template, "warped template") 

58 display.frame = 3 

59 display.image(template.variance, "warped variance") 

60 

61 

62class GetTemplateTaskTestCase(lsst.utils.tests.TestCase): 

63 """Test that GetTemplateTask works on both one tract and multiple tract 

64 input coadd exposures. 

65 

66 Makes a synthetic exposure large enough to fit four small tracts with 2x2 

67 (300x300 pixel) patches each, extracts pixels for those patches by warping, 

68 and tests GetTemplateTask's output against boxes that overlap various 

69 combinations of one or multiple tracts. 

70 """ 

71 def setUp(self): 

72 self.scale = 0.2 # arcsec/pixel 

73 self.skymap = self._makeSkymap() 

74 self.patches = collections.defaultdict(list) 

75 self.dataIds = collections.defaultdict(list) 

76 self.exposure = self._makeExposure() 

77 self.varianceBox = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)) 

78 

79 if debug: 79 ↛ 80line 79 didn't jump to line 80 because the condition on line 79 was never true

80 display.image(self.exposure, "base exposure") 

81 

82 for tract_id in range(4): 

83 tract = self.skymap.generateTract(tract_id) 

84 self._makePatches(tract) 

85 

86 def _makeSkymap(self): 

87 """Make a Skymap with 4 tracts with 4 patches each. 

88 """ 

89 tractScale = 0.02 # degrees 

90 # On-sky coordinates of the tract centers. 

91 coords = [(0, 0), 

92 (0, tractScale), 

93 (tractScale, 0), 

94 (tractScale, tractScale), 

95 ] 

96 config = lsst.skymap.DiscreteSkyMap.ConfigClass() 

97 config.raList = [c[0] for c in coords] 

98 config.decList = [c[1] for c in coords] 

99 # Half the tract center step size, to keep the tract overlap small. 

100 config.radiusList = [tractScale/2 for c in coords] 

101 config.projection = "TAN" 

102 config.pixelScale = self.scale 

103 config.tractOverlap = 0.0005 

104 config.tractBuilder = "legacy" 

105 config.tractBuilder["legacy"].patchInnerDimensions = (300, 300) 

106 config.tractBuilder["legacy"].patchBorder = 10 

107 return lsst.skymap.DiscreteSkyMap(config=config) 

108 

109 def _makeExposure(self): 

110 """Create a large image to break up into tracts and patches. 

111 

112 The image will have a source every 100 pixels in x and y, and a WCS 

113 that results in the tracts all fitting in the image, with tract=0 

114 in the lower left, tract=1 to the right, tract=2 above, and tract=3 

115 to the upper right. 

116 """ 

117 box = lsst.geom.Box2I(lsst.geom.Point2I(-200, -200), lsst.geom.Point2I(800, 800)) 

118 # This WCS was constructed so that tract 0 mostly fills the lower left 

119 # quadrant of the image, and the other tracts fill the rest; slight 

120 # extra rotation as a check on the final warp layout, scaled by 5% 

121 # from the patch pixel scale. 

122 cd_matrix = lsst.afw.geom.makeCdMatrix(1.05*self.scale*lsst.geom.arcseconds, 93*lsst.geom.degrees) 

123 wcs = lsst.afw.geom.makeSkyWcs(lsst.geom.Point2D(120, 150), 

124 lsst.geom.SpherePoint(0, 0, lsst.geom.radians), 

125 cd_matrix) 

126 dataset = lsst.meas.base.tests.TestDataset(box, wcs=wcs) 

127 for x, y in itertools.product(np.arange(0, 500, 100), np.arange(0, 500, 100)): 

128 dataset.addSource(1e5, lsst.geom.Point2D(x, y)) 

129 exposure, _ = dataset.realize(2, dataset.makeMinimalSchema()) 

130 exposure.setFilter(lsst.afw.image.FilterLabel("a", "a_test")) 

131 return exposure 

132 

133 def _makePatches(self, tract): 

134 """Populate the patches and dataId dicts, keyed on tract id, with the 

135 warps of the main exposure and minimal dataIds, respectively. 

136 """ 

137 if debug: 137 ↛ 138line 137 didn't jump to line 138 because the condition on line 137 was never true

138 color = ['red', 'green', 'cyan', 'yellow'][tract.tract_id] 

139 point = self.exposure.wcs.skyToPixel(tract.ctr_coord) 

140 # Show the tract center, colored by tract id. 

141 display.dot("x", point.x, point.y, ctype=color, size=30) 

142 

143 # Use 5th order to minimize artifacts on the templates. 

144 config = lsst.afw.math.Warper.ConfigClass() 

145 config.warpingKernelName = "lanczos5" 

146 warper = lsst.afw.math.Warper.fromConfig(config) 

147 for patchId in range(tract.num_patches.x*tract.num_patches.y): 

148 patch = tract.getPatchInfo(patchId) 

149 box = patch.getOuterBBox() 

150 

151 if debug: 151 ↛ 153line 151 didn't jump to line 153 because the condition on line 151 was never true

152 # Show the patch corners as patch ids, colored by tract id. 

153 points = self.exposure.wcs.skyToPixel(patch.wcs.pixelToSky([lsst.geom.Point2D(x) 

154 for x in box.getCorners()])) 

155 for p in points: 

156 display.dot(patchId, p.x, p.y, ctype=color) 

157 

158 # This is mostly taken from drp_tasks makePsfMatchedWarp, but 

159 # ip_diffim cannot depend on drp_tasks. 

160 xyTransform = lsst.afw.geom.makeWcsPairTransform(self.exposure.wcs, patch.wcs) 

161 warpedPsf = lsst.meas.algorithms.WarpedPsf(self.exposure.psf, xyTransform) 

162 warped = warper.warpExposure(patch.wcs, self.exposure, destBBox=box) 

163 warped.setPsf(warpedPsf) 

164 

165 warped.getInfo().setCoaddInputs( 

166 self._makeCoaddInputs([(self.exposure.wcs, self.exposure.getBBox())])) 

167 dataRef = pipeBase.InMemoryDatasetHandle( 

168 warped, 

169 storageClass="ExposureF", 

170 copy=True, 

171 dataId=generate_data_id( 

172 tract=tract, 

173 patch=patch, 

174 ) 

175 ) 

176 self.patches[tract.tract_id].append(dataRef) 

177 dataCoordinate = DataCoordinate.standardize({"tract": tract.tract_id, 

178 "patch": patchId, 

179 "band": "a", 

180 "skymap": "skymap"}, 

181 universe=DimensionUniverse()) 

182 self.dataIds[tract.tract_id].append(dataCoordinate) 

183 

184 def _checkMetadata(self, template, config, box, wcs, nPsfs): 

185 """Check that the various metadata components were set correctly. 

186 """ 

187 expectedBox = lsst.geom.Box2I(box) 

188 expectedBox.grow(config.templateBorderSize) 

189 self.assertEqual(template.getBBox(), expectedBox) 

190 # WCS should match our exposure, not any of the coadd tracts. 

191 for tract in self.patches: 

192 self.assertNotEqual(template.wcs, self.patches[tract][0].get().wcs) 

193 self.assertEqual(template.wcs, self.exposure.wcs) 

194 self.assertEqual(template.photoCalib, self.exposure.photoCalib) 

195 self.assertEqual(template.getXY0(), expectedBox.getMin()) 

196 self.assertEqual(template.filter.bandLabel, "a") 

197 self.assertEqual(template.filter.physicalLabel, "a_test") 

198 self.assertEqual(template.psf.getComponentCount(), nPsfs) 

199 self.assertTrue(template.getInfo().hasCoaddInputs()) 

200 self.assertEqual(len(template.getInfo().getCoaddInputs().ccds), nPsfs) 

201 

202 def _checkPixels(self, template, config, box): 

203 """Check that the pixel values in the template are close to the 

204 original image. 

205 """ 

206 # All pixels should have real values! 

207 expectedBox = lsst.geom.Box2I(box) 

208 expectedBox.grow(config.templateBorderSize) 

209 

210 if debug: 210 ↛ 211line 210 didn't jump to line 211 because the condition on line 210 was never true

211 _showTemplate(expectedBox, template) 

212 

213 # Check that we fully filled the template from the patches. 

214 self.assertTrue(np.all(np.isfinite(template.image.array))) 

215 # Because of the scale changes, there will be some ringing in the 

216 # difference between the template and the original image; pick 

217 # tolerances large enough to account for that. 

218 self.assertImagesAlmostEqual(template.image, self.exposure[expectedBox].image, 

219 rtol=.1, atol=4) 

220 # Variance plane ==4 in the original image (realize() takes a noise 

221 # sigma). Warping sets the level from the pixel areas and the two 

222 # warping kernels, which `_correctVariance` corrects up to the 

223 # difference between the configured coaddWarpKernel and the lanczos5 

224 # `_makePatches` really used. A per-pixel ripple that no scalar 

225 # correction can remove remains on top of that, so check the level 

226 # tightly and allow for the ripple around it. 

227 variance = template.variance.array[np.isfinite(template.variance.array)] 

228 median = np.median(variance) 

229 self.assertFloatsAlmostEqual(median, 

230 np.median(self.exposure[expectedBox].variance.array), 

231 rtol=0.35, msg="variance level differs") 

232 self.assertLess(np.percentile(variance, 99)/median, 1.6, msg="variance ripple too large") 

233 self.assertGreater(np.percentile(variance, 1)/median, 0.6, msg="variance ripple too large") 

234 # Not checking the mask, as warping changes the sizes of the masks. 

235 

236 def testRunOneTractInput(self): 

237 """Test a bounding box that fully fits inside one tract, with only 

238 that tract passed as input. This checks that the code handles a single 

239 tract input correctly. 

240 """ 

241 box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)) 

242 task = lsst.ip.diffim.GetTemplateTask() 

243 # Restrict to tract 0, since the box fits in just that tract. 

244 # Task modifies the input bbox, so pass a copy. 

245 result = task.run(coaddExposureHandles={0: self.patches[0]}, 

246 bbox=lsst.geom.Box2I(box), 

247 wcs=self.exposure.wcs, 

248 dataIds={0: self.dataIds[0]}, 

249 physical_filter="a_test") 

250 

251 # All 4 patches from tract 0 are included in this template. 

252 self._checkMetadata(result.template, task.config, box, self.exposure.wcs, 4) 

253 self._checkPixels(result.template, task.config, box) 

254 

255 def testRunOneTractMultipleInputs(self): 

256 """Test a bounding box that fully fits inside one tract but where 

257 multiple tracts were passed in. This checks that patches that are 

258 mostly NaN after warping are merged correctly in the output. 

259 """ 

260 box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)) 

261 task = lsst.ip.diffim.GetTemplateTask() 

262 # Task modifies the input bbox, so pass a copy. 

263 result = task.run(coaddExposureHandles=self.patches, 

264 bbox=lsst.geom.Box2I(box), 

265 wcs=self.exposure.wcs, 

266 dataIds=self.dataIds, 

267 physical_filter="a_test") 

268 

269 # All 4 patches from two tracts are included in this template. 

270 self._checkMetadata(result.template, task.config, box, self.exposure.wcs, 6) 

271 self._checkPixels(result.template, task.config, box) 

272 

273 def testRunTwoTracts(self): 

274 """Test a bounding box that crosses tract boundaries. 

275 """ 

276 box = lsst.geom.Box2I(lsst.geom.Point2I(200, 200), lsst.geom.Point2I(600, 600)) 

277 task = lsst.ip.diffim.GetTemplateTask() 

278 # Task modifies the input bbox, so pass a copy. 

279 result = task.run(coaddExposureHandles=self.patches, 

280 bbox=lsst.geom.Box2I(box), 

281 wcs=self.exposure.wcs, 

282 dataIds=self.dataIds, 

283 physical_filter="a_test") 

284 

285 # All 4 patches from all 4 tracts are included in this template 

286 self._checkMetadata(result.template, task.config, box, self.exposure.wcs, 9) 

287 self._checkPixels(result.template, task.config, box) 

288 

289 def testRunNoTemplate(self): 

290 """A bounding box that doesn't overlap the patches will raise. 

291 """ 

292 box = lsst.geom.Box2I(lsst.geom.Point2I(1200, 1200), lsst.geom.Point2I(1600, 1600)) 

293 task = lsst.ip.diffim.GetTemplateTask() 

294 with self.assertRaisesRegex(lsst.pipe.base.NoWorkFound, "No patches found"): 

295 task.run(coaddExposureHandles=self.patches, 

296 bbox=lsst.geom.Box2I(box), 

297 wcs=self.exposure.wcs, 

298 dataIds=self.dataIds, 

299 physical_filter="a_test") 

300 

301 def testMissingPatches(self): 

302 """Test that a missing patch results in an appropriate mask. 

303 

304 This fixes the bug reported on DM-44997 (image and variance were NaN 

305 but the mask was not set to NO_DATA for those pixels). 

306 """ 

307 # tract=0, patch=1 is the lower-left corner, as displayed in DS9. 

308 self.patches[0].pop(1) 

309 box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)) 

310 task = lsst.ip.diffim.GetTemplateTask() 

311 # Task modifies the input bbox, so pass a copy. 

312 result = task.run(coaddExposureHandles=self.patches, 

313 bbox=lsst.geom.Box2I(box), 

314 wcs=self.exposure.wcs, 

315 dataIds=self.dataIds, 

316 physical_filter="a_test") 

317 no_data = (result.template.mask.array & result.template.mask.getPlaneBitMask("NO_DATA")) != 0 

318 self.assertTrue(np.isfinite(result.template.image.array).all()) 

319 self.assertTrue(np.isfinite(result.template.variance.array).all()) 

320 self.assertEqual(no_data.sum(), 20990) 

321 

322 @lsst.utils.tests.methodParameters( 

323 box=[ 

324 lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)), 

325 lsst.geom.Box2I(lsst.geom.Point2I(200, 200), lsst.geom.Point2I(600, 600)), 

326 ], 

327 nInput=[8, 16], 

328 ) 

329 def testNanInputs(self, box=None, nInput=None): 

330 """Test that the template has finite values when some of the input 

331 pixels have NaN as variance. 

332 """ 

333 for tract, patchRefs in self.patches.items(): 

334 for patchRef in patchRefs: 

335 patchCoadd = patchRef.get() 

336 bbox = lsst.geom.Box2I() 

337 bbox.include(lsst.geom.Point2I(patchCoadd.getBBox().getCenter())) 

338 bbox.grow(3) 

339 patchCoadd.variance[bbox].array *= np.nan 

340 

341 box = lsst.geom.Box2I(lsst.geom.Point2I(200, 200), lsst.geom.Point2I(600, 600)) 

342 task = lsst.ip.diffim.GetTemplateTask() 

343 result = task.run(coaddExposureHandles=self.patches, 

344 bbox=lsst.geom.Box2I(box), 

345 wcs=self.exposure.wcs, 

346 dataIds=self.dataIds, 

347 physical_filter="a_test") 

348 if debug: 348 ↛ 349line 348 didn't jump to line 349 because the condition on line 348 was never true

349 _showTemplate(box, result.template) 

350 self._checkMetadata(result.template, task.config, box, self.exposure.wcs, 9) 

351 # We just check that the pixel values are all finite. We cannot check that pixel values 

352 # in the template are closer to the original anymore. 

353 self.assertTrue(np.isfinite(result.template.image.array).all()) 

354 

355 def _runCorrection(self, doCorrectVariancePlateScale=False, doScaleVariance=False, 

356 handles=None): 

357 """Run the task on tract 0 with the named variance corrections. 

358 

359 Everything defaults to off, so each test turns on exactly what it is 

360 exercising. 

361 """ 

362 config = lsst.ip.diffim.GetTemplateTask.ConfigClass() 

363 config.doCorrectVariancePlateScale = doCorrectVariancePlateScale 

364 config.doScaleVariance = doScaleVariance 

365 task = lsst.ip.diffim.GetTemplateTask(config=config) 

366 result = task.run(coaddExposureHandles={0: handles or self.patches[0]}, 

367 bbox=lsst.geom.Box2I(self.varianceBox), 

368 wcs=self.exposure.wcs, 

369 dataIds={0: self.dataIds[0]}, 

370 physical_filter="a_test") 

371 return task, result.template 

372 

373 def _makeScaledWcs(self, factor): 

374 """Make a WCS like the base exposure's, but with its pixel scale 

375 multiplied by ``factor``. 

376 """ 

377 cdMatrix = lsst.afw.geom.makeCdMatrix(factor*1.05*self.scale*lsst.geom.arcseconds, 

378 93*lsst.geom.degrees) 

379 return lsst.afw.geom.makeSkyWcs(lsst.geom.Point2D(120, 150), 

380 lsst.geom.SpherePoint(0, 0, lsst.geom.radians), 

381 cdMatrix) 

382 

383 @staticmethod 

384 def _makeCoaddInputs(records): 

385 """Make a CoaddInputs holding the given input records. 

386 

387 Parameters 

388 ---------- 

389 records : `list` [`tuple` [`lsst.afw.geom.SkyWcs` or `None`, \ 

390 `lsst.geom.Box2I`]] 

391 The WCS and bbox to record for each input. A `None` WCS makes a 

392 record that `_plateScaleFactor` has to skip. May be empty, to 

393 simulate a coadd whose plate scale cannot be reconstructed. 

394 """ 

395 ccdSchema = lsst.afw.table.ExposureTable.makeMinimalSchema() 

396 weightKey = ccdSchema.addField("weight", type=float, doc="Coadd weight") 

397 coaddInputs = lsst.afw.image.CoaddInputs( 

398 lsst.afw.table.ExposureTable.makeMinimalSchema(), ccdSchema) 

399 for wcs, bbox in records: 

400 record = coaddInputs.ccds.addNew() 

401 record.setWcs(wcs) 

402 record.setBBox(bbox) 

403 # Included because real coadds have it, though a single-record 

404 # correction does not use it. 

405 record.set(weightKey, 1.0) 

406 return coaddInputs 

407 

408 def _patchHandles(self, tract, records): 

409 """Return handles for a tract's patches, with their CoaddInputs 

410 replaced by ``records``. 

411 """ 

412 handles = [] 

413 for ref in self.patches[tract]: 

414 coadd = ref.get() 

415 coadd.getInfo().setCoaddInputs(self._makeCoaddInputs(records)) 

416 handles.append(pipeBase.InMemoryDatasetHandle( 

417 coadd, storageClass="ExposureF", copy=True, dataId=ref.dataId)) 

418 return handles 

419 

420 def testCorrectVariancePlateScale(self): 

421 """The plate scale correction is the total pixel area change from the 

422 images the coadds were built from to the science image. 

423 """ 

424 _, off = self._runCorrection() 

425 

426 # The fixture's records are the science image itself, so there is no 

427 # net change in pixel area: the same-instrument case. 

428 task, on = self._runCorrection(doCorrectVariancePlateScale=True) 

429 self.assertFloatsAlmostEqual(task.metadata["variancePlateScaleFactor"], 1.0, rtol=1e-6) 

430 

431 # Coarser original pixels than science pixels, the DECam-template 

432 # case: the correction is the ratio of their areas. 

433 for scaleFactor in (1.315, 0.5): 

434 with self.subTest(scaleFactor=scaleFactor): 

435 handles = self._patchHandles( 

436 0, [(self._makeScaledWcs(scaleFactor), self.exposure.getBBox())]) 

437 task, on = self._runCorrection(doCorrectVariancePlateScale=True, 

438 handles=handles) 

439 factor = task.metadata["variancePlateScaleFactor"] 

440 self.assertFloatsAlmostEqual(factor, scaleFactor**2, rtol=1e-6) 

441 self.assertFloatsAlmostEqual(on.variance.array, off.variance.array*factor, 

442 rtol=1e-5, ignoreNaNs=True) 

443 

444 def testCorrectVariancePlateScaleUsesOneRecord(self): 

445 """Only the first usable coadd input is read. 

446 

447 The spread between records is just the local pixel scale, a few 

448 tenths of a percent across a real focal plane, so one stands for all 

449 of them. 

450 """ 

451 handles = self._patchHandles( 

452 0, [(None, self.exposure.getBBox()), 

453 (self._makeScaledWcs(2.0), self.exposure.getBBox()), 

454 (self.exposure.wcs, self.exposure.getBBox())]) 

455 task, _ = self._runCorrection(doCorrectVariancePlateScale=True, handles=handles) 

456 

457 # The first record with a WCS: the one without is skipped, and the 

458 # last (which would give 1.0) is never reached. 

459 self.assertFloatsAlmostEqual(task.metadata["variancePlateScaleFactor"], 4.0, rtol=1e-6) 

460 

461 def testCorrectVariancePlateScaleNeedsCoaddInputs(self): 

462 """Without usable coadd inputs the plate scale cannot be 

463 reconstructed, and the task must say so rather than silently applying 

464 only part of the correction. 

465 """ 

466 handles = self._patchHandles(0, []) 

467 with self.assertRaisesRegex(RuntimeError, "doCorrectVariancePlateScale"): 

468 self._runCorrection(doCorrectVariancePlateScale=True, handles=handles) 

469 

470 def _scaleInputVariance(self, tract, factor): 

471 """Return fresh handles for one tract's patches, with their variance 

472 planes multiplied by ``factor``. 

473 

474 Parameters 

475 ---------- 

476 tract : `int` 

477 Id of the tract whose patches should be copied. 

478 factor : `float` 

479 Factor to multiply the input variance planes by. 

480 

481 Returns 

482 ------- 

483 handles : `list` [`lsst.pipe.base.InMemoryDatasetHandle`] 

484 Handles to the modified patches. 

485 """ 

486 handles = [] 

487 for handle in self.patches[tract]: 

488 # ``copy=True`` on the original handles means this is a copy, so 

489 # the patches shared with the other tests are left untouched. 

490 patch = handle.get() 

491 patch.variance.array *= factor 

492 handles.append(pipeBase.InMemoryDatasetHandle(patch, 

493 storageClass="ExposureF", 

494 copy=True, 

495 dataId=handle.dataId)) 

496 return handles 

497 

498 def testScaleVariance(self): 

499 """Test that the template variance plane is rescaled to match the 

500 empirical pixel noise, and that the factor used is recorded in the 

501 task metadata. 

502 """ 

503 scaleFactor = 1.345 

504 box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)) 

505 

506 def _configureAndRunTask(doScaleVariance, varianceScale=1.): 

507 """Build a template from tract 0, optionally rescaling the input 

508 variance planes by ``varianceScale`` first. 

509 """ 

510 config = lsst.ip.diffim.GetTemplateTask.ConfigClass() 

511 config.doScaleVariance = doScaleVariance 

512 task = lsst.ip.diffim.GetTemplateTask(config=config) 

513 # Task modifies the input bbox, so pass a copy. 

514 result = task.run(coaddExposureHandles={0: self._scaleInputVariance(0, varianceScale)}, 

515 bbox=lsst.geom.Box2I(box), 

516 wcs=self.exposure.wcs, 

517 dataIds={0: self.dataIds[0]}, 

518 physical_filter="a_test") 

519 return task, result.template 

520 

521 # With scaling disabled the subtask is never constructed, and nothing 

522 # is recorded in the metadata. 

523 taskOff, templateOff = _configureAndRunTask(False) 

524 self.assertFalse(hasattr(taskOff, "scaleVariance")) 

525 self.assertNotIn("scaleTemplateVarianceFactor", taskOff.metadata) 

526 

527 # Both warps -- lanczos5 in ``_makePatches`` and lanczos3 in the 

528 # task -- correlate the noise. The variance plane tracks only the 

529 # per-pixel diagonal, which the second warp leaves too low, so 

530 # ``scaleVariance`` measures a factor well above 1 even though the 

531 # input variance planes are correct. 

532 # 

533 taskOn, templateOn = _configureAndRunTask(True) 

534 factor = taskOn.metadata["scaleTemplateVarianceFactor"] 

535 # TODO DM-55879: this value is pinned on purpose. The lanczos warping 

536 # kernels introduce small correlations that artificially suppress the 

537 # image pixel stddev and inflate the variance scaling factor. This 

538 # should be changed to 1.0 after DM-55879 is merged. 

539 self.assertFloatsAlmostEqual(factor, 1.1465, atol=0.01, 

540 msg="Measured template variance scaling changed; see the" 

541 " comment above if the correlation correction landed.") 

542 # The only difference from the unscaled template is the constant 

543 # factor applied to the variance plane. 

544 self.assertFloatsAlmostEqual(templateOn.variance.array, 

545 templateOff.variance.array*factor, rtol=1e-5) 

546 # Tolerance here is float32 round-off: repeated runs of the task are 

547 # not bitwise identical. 

548 self.assertImagesAlmostEqual(templateOn.image, templateOff.image, rtol=1e-5, atol=1e-5) 

549 

550 # If the input variance planes under-estimate the noise by a known 

551 # factor, the measured factor grows by that amount and the same 

552 # output variance plane is recovered. 

553 taskLow, templateLow = _configureAndRunTask(True, varianceScale=1/scaleFactor) 

554 self.assertFloatsAlmostEqual(taskLow.metadata["scaleTemplateVarianceFactor"], 

555 factor*scaleFactor, rtol=1e-5) 

556 self.assertImagesAlmostEqual(templateLow.variance, templateOn.variance, rtol=1e-5) 

557 

558 

559def setup_module(module): 

560 lsst.utils.tests.init() 

561 

562 

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

564 pass 

565 

566 

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

568 lsst.utils.tests.init() 

569 unittest.main()