Coverage for tests/test_fields.py: 97%

1# This file is part of lsst-images. 

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# Use of this source code is governed by a 3-clause BSD-style 

10# license that can be found in the LICENSE file. 

11 

12from __future__ import annotations 

13 

14import os 

15import unittest 

16 

17import astropy.units as u 

18import numpy as np 

19 

20from lsst.images import Box, Image 

21from lsst.images.fields import ( 

22 BaseField, 

23 ChebyshevField, 

24 ProductField, 

25 SplineField, 

26 SumField, 

27 field_from_legacy, 

28 field_from_legacy_background, 

29) 

30from lsst.images.tests import ( 

31 RoundtripFits, 

32 assert_close, 

33 assert_images_equal, 

34 compare_field_to_legacy, 

35) 

36 

37try: 

38 from lsst.afw.math import BackgroundList as LegacyBackgroundList 

39 from lsst.afw.math import Chebyshev1Function2D as LegacyChebyshev1Function2D 

40 from lsst.afw.math import ChebyshevBoundedField as LegacyChebyshevBoundedField 

41 from lsst.afw.math import ProductBoundedField as LegacyProductBoundedField 

42 from lsst.geom import Box2D as LegacyBox2D 

43except ImportError: 

44 HAVE_LEGACY = False 

45else: 

46 HAVE_LEGACY = True 

47 

48DATA_DIR = os.environ.get("TESTDATA_IMAGES_DIR", None) 

49 

50 

51class FieldTestCase(unittest.TestCase): 

52 """Tests for the Field classes that do not require legacy code.""" 

53 

54 def setUp(self) -> None: 

55 self.rng = np.random.default_rng(10) 

56 self.box = Box.factory[6:32, -7:26] 

57 self.cheby = ChebyshevField(self.box, np.array([[0.5, -0.25], [0.40, 0.0]])) 

58 self.spline = SplineField( 

59 self.box, 

60 self.rng.standard_normal(size=(6, 7)), 

61 y=self.box.y.linspace(6), 

62 x=self.box.x.linspace(7), 

63 ) 

64 self.product = ProductField([self.spline, self.cheby]) 

65 self.sum = SumField([self.spline, self.cheby]) 

66 

67 def test_chebyshev_call_limits(self) -> None: 

68 """Test that ChebyshevField.__call__ evaluates correctly at low order 

69 at the corners of its box. 

70 

71 Testing at first order is enough to catch (y, x) swap errors, and is 

72 easy because ``T_0(x) = 1`` and ``T_1(x) = x``. Higher-order 

73 evaluation is covered in the legacy comparison tests. 

74 """ 

75 result = self.cheby(x=np.array([-7.5, 25.5, 25.5, -7.5]), y=np.array([5.5, 5.5, 31.5, 31.5])) 

76 self.assertEqual(result[0], 0.5 + 0.25 - 0.4) # [x=-1, y=-1] after remap 

77 self.assertEqual(result[1], 0.5 - 0.25 - 0.4) # [x=1, y=-1] after remap 

78 self.assertEqual(result[2], 0.5 - 0.25 + 0.4) # [x=1, y=1] after remap 

79 self.assertEqual(result[3], 0.5 + 0.25 + 0.4) # [x=-1, y=1] after remap 

80 

81 def test_chebyshev_attributes(self) -> None: 

82 """Test the basic properties of a Chebyshev field.""" 

83 self.assertEqual(self.cheby.bounds, self.box) 

84 self.assertIsNone(self.cheby.unit) 

85 self.assertEqual(self.cheby.x_order, 1) 

86 self.assertEqual(self.cheby.y_order, 1) 

87 self.assertEqual(self.cheby.order, 1) 

88 np.testing.assert_array_equal(self.cheby.coefficients, np.array([[0.5, -0.25], [0.40, 0.0]])) 

89 

90 def test_chebyshev_fit(self) -> None: 

91 """Test that we can fit a ChebyshevField to gridded data that should 

92 have zero residuals. 

93 """ 

94 data_image = self.cheby.render() 

95 cheby2 = ChebyshevField.fit( 

96 self.box, data_image.array, order=1, y=self.box.y.arange, x=self.box.x.arange 

97 ) 

98 self.assertEqual(cheby2.bounds, self.box) 

99 np.testing.assert_array_almost_equal(cheby2.coefficients, self.cheby.coefficients) 

100 # Repeat the experiment in a few different ways that don't affect the 

101 # result (at least in this scenario where we have an exact fit). 

102 # Fit to order 2 in x (will get us extra zero-valued coefficients): 

103 cheby3 = ChebyshevField.fit( 

104 self.box, data_image.array, x_order=2, y_order=1, y=self.box.y.arange, x=self.box.x.arange 

105 ) 

106 self.assertEqual(cheby3.bounds, self.box) 

107 np.testing.assert_array_almost_equal( 

108 cheby3.coefficients, 

109 np.array([[0.5, -0.25, 0.0], [0.40, 0.0, 0.0]], dtype=np.float64), 

110 ) 

111 # Fit with triangular=False (would allow the (1, 1) term to be nonzero, 

112 # but it will still be zero here): 

113 cheby4 = ChebyshevField.fit( 

114 self.box, data_image.array, order=1, y=self.box.y.arange, x=self.box.x.arange, triangular=False 

115 ) 

116 self.assertEqual(cheby4.bounds, self.box) 

117 np.testing.assert_array_almost_equal(cheby4.coefficients, self.cheby.coefficients) 

118 # Fit with weights. 

119 cheby5 = ChebyshevField.fit( 

120 self.box, 

121 data_image.array, 

122 order=1, 

123 y=self.box.y.arange, 

124 x=self.box.x.arange, 

125 weight=self.rng.uniform(0.8, 1.2, size=data_image.array.shape), 

126 ) 

127 self.assertEqual(cheby5.bounds, self.box) 

128 np.testing.assert_array_almost_equal(cheby5.coefficients, self.cheby.coefficients) 

129 # Fit to a Quantity. 

130 cheby6 = ChebyshevField.fit( 

131 self.box, data_image.array * u.nJy, order=1, y=self.box.y.arange, x=self.box.x.arange 

132 ) 

133 self.assertEqual(cheby6.bounds, self.box) 

134 self.assertEqual(cheby6.unit, u.nJy) 

135 np.testing.assert_array_almost_equal(cheby6.coefficients, self.cheby.coefficients) 

136 # Fit with units provided separately. 

137 cheby7 = ChebyshevField.fit( 

138 self.box, data_image.array, order=1, y=self.box.y.arange, x=self.box.x.arange, unit=u.nJy 

139 ) 

140 self.assertEqual(cheby7.bounds, self.box) 

141 self.assertEqual(cheby7.unit, u.nJy) 

142 np.testing.assert_array_almost_equal(cheby7.coefficients, self.cheby.coefficients) 

143 # Fit with x and y labeling every data point, not the grid. 

144 m = self.box.meshgrid() 

145 cheby8 = ChebyshevField.fit(self.box, data_image.array, order=1, y=m.y, x=m.x) 

146 self.assertEqual(cheby8.bounds, self.box) 

147 np.testing.assert_array_almost_equal(cheby8.coefficients, self.cheby.coefficients) 

148 # Fit with x and y labeling every data point, not the grid, as well as 

149 # weights. 

150 cheby9 = ChebyshevField.fit( 

151 self.box, 

152 data_image.array, 

153 order=1, 

154 y=m.y, 

155 x=m.x, 

156 weight=self.rng.uniform(0.8, 1.2, size=data_image.array.shape), 

157 ) 

158 self.assertEqual(cheby9.bounds, self.box) 

159 np.testing.assert_array_almost_equal(cheby9.coefficients, self.cheby.coefficients) 

160 # Fit with one data points replaced by NaN, which should be ignored 

161 # (the fit is still overconstrained). 

162 new_data = data_image.array.copy() 

163 new_data[5, 7] = np.nan 

164 cheby10 = ChebyshevField.fit(self.box, new_data, order=1, y=self.box.y.arange, x=self.box.x.arange) 

165 self.assertEqual(cheby10.bounds, self.box) 

166 np.testing.assert_array_almost_equal(cheby10.coefficients, self.cheby.coefficients) 

167 

168 def test_chebyshev_evaluation_consistency(self) -> None: 

169 self.check_evaluation_consistency(self.cheby) 

170 

171 def test_chebyshev_units(self) -> None: 

172 self.check_units(self.cheby) 

173 

174 def test_spline_knot_evaluation(self) -> None: 

175 """Test that SplineField.__call__ evaluates to the input data points 

176 at their positions. 

177 """ 

178 xv, yv = np.meshgrid(self.spline.x, self.spline.y) 

179 result = self.spline(x=xv, y=yv) 

180 np.testing.assert_array_almost_equal(result, self.spline.data) 

181 

182 def test_spline_evaluation_consistency(self) -> None: 

183 self.check_evaluation_consistency(self.spline) 

184 

185 def test_spline_units(self) -> None: 

186 self.check_units(self.spline) 

187 

188 def test_product_evaluation(self) -> None: 

189 """Test ProductField.__call__ against direct calls to its operands.""" 

190 xv, yv = self.box.meshgrid(n=3) 

191 z = self.product(x=xv, y=yv) 

192 np.testing.assert_array_equal(z, self.cheby(x=xv, y=yv) * self.spline(x=xv, y=yv)) 

193 

194 def test_product_evaluation_consistency(self) -> None: 

195 self.check_evaluation_consistency(self.product) 

196 

197 def test_product_units(self) -> None: 

198 self.check_units(self.product) 

199 self.assertEqual(ProductField([self.cheby, self.spline * u.nJy]).unit, u.nJy) 

200 self.assertEqual(ProductField([self.cheby * u.nJy, self.spline]).unit, u.nJy) 

201 self.assertEqual( 

202 ProductField([self.cheby * u.nJy, self.spline / u.arcsec**2]).unit, u.nJy / u.arcsec**2 

203 ) 

204 

205 def test_sum_evaluation(self) -> None: 

206 """Test sumField.__call__ against direct calls to its operands.""" 

207 xv, yv = self.box.meshgrid(n=3) 

208 z = self.sum(x=xv, y=yv) 

209 np.testing.assert_array_equal(z, self.cheby(x=xv, y=yv) + self.spline(x=xv, y=yv)) 

210 

211 def test_sum_evaluation_consistency(self) -> None: 

212 self.check_evaluation_consistency(self.sum) 

213 

214 def test_sum_units(self) -> None: 

215 self.check_units(self.sum) 

216 with self.assertRaises(u.UnitConversionError): 

217 SumField([self.cheby, self.spline * u.nJy]) 

218 with self.assertRaises(u.UnitConversionError): 

219 SumField([self.cheby * u.nJy, self.spline]) 

220 # Test a SumField where the operands have different but compatible 

221 # units. 

222 mixed = SumField([self.cheby * u.rad, self.spline * u.deg]) 

223 self.assertEqual(mixed.unit, u.rad) 

224 small_box = Box.factory[10:12, 2:5] 

225 cheby_render = self.cheby.render(small_box) 

226 spline_render = self.spline.render(small_box) 

227 mixed_render = mixed.render(small_box) 

228 self.assertEqual(mixed_render.unit, u.rad) 

229 np.testing.assert_array_almost_equal( 

230 mixed_render.array, cheby_render.array + (spline_render.array * np.pi / 180.0) 

231 ) 

232 self.check_evaluation_consistency(mixed) 

233 

234 def check_evaluation_consistency(self, field: BaseField) -> None: 

235 """Test that `BaseField.__call__` and `BaseField.render` agree on a 

236 concrete field. 

237 """ 

238 image_1 = field.render() 

239 p = field.bounds.bbox.meshgrid() 

240 image_2 = Image(field(x=p.x, y=p.y), bbox=field.bounds.bbox, unit=field.unit) 

241 assert_images_equal(self, image_1, image_2) 

242 scaled_field = field * 2.0 

243 image_3 = scaled_field.render() 

244 image_3.array *= 0.5 

245 assert_images_equal(self, image_1, image_3) 

246 image_4 = field.render(Box.factory[9:11, -3:1]) 

247 assert_images_equal(self, image_4, image_1[image_4.bbox]) 

248 

249 def check_units(self, field: BaseField) -> None: 

250 """Test that the methods of a `BaseField` implementation correctly 

251 propogate and check units. 

252 """ 

253 self.assertIsNone(field.unit) 

254 with_units_1 = field * u.nJy 

255 self.assertEqual(with_units_1.unit, u.nJy) 

256 self.assertEqual(with_units_1(x=np.array([0.0]), y=np.array([10.0]), quantity=True).unit, u.nJy) 

257 image_1 = with_units_1.render(bbox=Box.factory[10:12, 0:3]) 

258 self.assertEqual(image_1.unit, u.nJy) 

259 self.assertEqual((with_units_1 * 2.0).unit, u.nJy) 

260 self.assertEqual((with_units_1 / u.arcsec**2).unit, u.nJy / u.arcsec**2) 

261 

262 

263@unittest.skipUnless(HAVE_LEGACY, "This test requires lsst.afw.math to be importable.") 

264class FieldLegacyTestCase(unittest.TestCase): 

265 """Test the Field classes against legacy implementations. 

266 

267 This includes many tests for correct evaluation, since the legacy types 

268 serve as our reference implementation. 

269 """ 

270 

271 def setUp(self) -> None: 

272 self.rng = np.random.default_rng(10) 

273 self.box = Box.factory[6:32, -7:26] 

274 # This Chebyshev coefficient matrix is unusual in that it has nonzero 

275 # entries for the whole matrix, not just the p + q < order triangle, 

276 # and it has different orders for y and x. 

277 # But we want to make sure those round-trip, too. 

278 self.cheby_coeffs = self.rng.random((6, 3)) 

279 self.legacy_cheby = LegacyChebyshevBoundedField(self.box.to_legacy(), self.cheby_coeffs) 

280 cheby2 = LegacyChebyshevBoundedField(self.box.to_legacy(), self.rng.standard_normal(size=(2, 2))) 

281 self.legacy_product = LegacyProductBoundedField([self.legacy_cheby, cheby2]) 

282 

283 def test_chebyshev_roundtrip(self) -> None: 

284 """Test converting ChebyshevField from and to legacy, and serializing 

285 it in between. 

286 """ 

287 cheby = field_from_legacy(self.legacy_cheby) 

288 assert isinstance(cheby, ChebyshevField) 

289 compare_field_to_legacy(self, cheby, self.legacy_cheby, subimage_bbox=Box.factory[8:12, -3:2]) 

290 with RoundtripFits(self, cheby) as roundtrip: 

291 pass 

292 compare_field_to_legacy( 

293 self, roundtrip.result, self.legacy_cheby, subimage_bbox=Box.factory[8:12, -3:2] 

294 ) 

295 compare_field_to_legacy( 

296 self, roundtrip.result, cheby.to_legacy(), subimage_bbox=Box.factory[8:12, -3:2] 

297 ) 

298 

299 def test_product_roundtrip(self) -> None: 

300 """Test converting ProductField from and to legacy, and serializing 

301 it in between. 

302 """ 

303 product = field_from_legacy(self.legacy_product) 

304 assert isinstance(product, ProductField) 

305 compare_field_to_legacy(self, product, self.legacy_product, subimage_bbox=Box.factory[8:12, -3:2]) 

306 with RoundtripFits(self, product) as roundtrip: 

307 pass 

308 compare_field_to_legacy( 

309 self, roundtrip.result, self.legacy_product, subimage_bbox=Box.factory[8:12, -3:2] 

310 ) 

311 compare_field_to_legacy( 

312 self, roundtrip.result, product.to_legacy(), subimage_bbox=Box.factory[8:12, -3:2] 

313 ) 

314 

315 def test_spline_simple(self) -> None: 

316 """Test SplineField against `lsst.afw.math.BackgroundMI`, when there 

317 is no missing data. 

318 """ 

319 from lsst.afw.image import MaskedImageF 

320 from lsst.afw.math import BackgroundMI 

321 

322 bins = MaskedImageF(Box.factory[0:5, 0:6].to_legacy()) 

323 bins.image.array[:, :] = self.rng.standard_normal(bins.image.array.shape) 

324 bins.variance.array[::] = 1.0 

325 legacy_bg = BackgroundMI(self.box.to_legacy(), bins) 

326 spline = field_from_legacy_background(legacy_bg) 

327 render_bbox = self.box.padded(-3) 

328 assert_images_equal( 

329 self, 

330 spline.render(render_bbox), 

331 Image.from_legacy( 

332 legacy_bg.getImageF( 

333 render_bbox.to_legacy(), legacy_bg.getBackgroundControl().getInterpStyle() 

334 ) 

335 ), 

336 rtol=1e-7, 

337 ) 

338 

339 def test_spline_one_nan(self) -> None: 

340 """Test SplineField against `lsst.afw.math.BackgroundMI`, when there 

341 is missing data. 

342 """ 

343 from lsst.afw.image import MaskedImageF 

344 from lsst.afw.math import BackgroundMI 

345 

346 bins = MaskedImageF(Box.factory[0:7, 0:6].to_legacy()) 

347 bins.image.array[:, :] = self.rng.standard_normal(bins.image.array.shape) 

348 bins.image.array[3, 2] = np.nan 

349 bins.variance.array[::] = 1.0 

350 legacy_bg = BackgroundMI(self.box.to_legacy(), bins) 

351 spline = field_from_legacy_background(legacy_bg) 

352 render_bbox = self.box.padded(-3) 

353 assert_images_equal( 

354 self, 

355 spline.render(render_bbox), 

356 Image.from_legacy( 

357 legacy_bg.getImageF( 

358 render_bbox.to_legacy(), 

359 legacy_bg.getBackgroundControl().getInterpStyle(), 

360 ) 

361 ), 

362 rtol=1e-7, 

363 ) 

364 

365 def test_chebyshev1_function2(self) -> None: 

366 legacy_func2a = LegacyChebyshev1Function2D(4, LegacyBox2D(self.box.to_legacy())) 

367 legacy_func2a.setParameters(self.rng.standard_normal(legacy_func2a.getNParameters())) 

368 field = ChebyshevField.from_legacy_function2(legacy_func2a) 

369 legacy_func2b = field.to_legacy_function2() 

370 self.assertEqual(field.bounds, self.box) 

371 xy_array = self.box.meshgrid(4) 

372 z_array = field(x=xy_array.x, y=xy_array.y) 

373 for z, x, y in zip(z_array.flat, xy_array.x.flat, xy_array.y.flat): 

374 assert_close(self, legacy_func2a(x, y), z) 

375 assert_close(self, legacy_func2b(x, y), z) 

376 

377 

378@unittest.skipIf(DATA_DIR is None, "This test requires the TESTDATA_IMAGES_DIR envvar to be set.") 

379@unittest.skipUnless(HAVE_LEGACY, "This test requires lsst.afw.math to be importable.") 

380class FieldLegacyDataTestCase(unittest.TestCase): 

381 """Tests for using Field classes using legacy datasets.""" 

382 

383 def test_visit_background(self) -> None: 

384 assert DATA_DIR is not None, "Guaranteed by decorator." 

385 filename = os.path.join(DATA_DIR, "dp2", "legacy", "visit_image_background.fits") 

386 legacy_bg_list = LegacyBackgroundList.readFits(filename) 

387 bg_field = field_from_legacy_background(legacy_bg_list) 

388 assert_images_equal(self, bg_field.render(), Image.from_legacy(legacy_bg_list.getImage()), rtol=1e-6) 

389 

390 

391if __name__ == "__main__": 

392 unittest.main()