Coverage for tests / test_CompensatedTophatFlux.py: 15%

1# This file is part of meas_base. 

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 

23 

24import numpy as np 

25 

26import lsst.geom 

27import lsst.afw.geom 

28import lsst.meas.base 

29import lsst.utils.tests 

30from lsst.meas.base.tests import AlgorithmTestCase 

31 

32 

33class CompensatedTophatFluxTestCase(AlgorithmTestCase, lsst.utils.tests.TestCase): 

34 def setUp(self): 

35 self.bbox = lsst.geom.Box2I(lsst.geom.Point2I(-20, -30), 

36 lsst.geom.Extent2I(1000, 1000)) 

37 

38 self.dataset = lsst.meas.base.tests.TestDataset(self.bbox) 

39 self.psf_size = 2.0 

40 self.dataset.psfShape = lsst.afw.geom.Quadrupole(self.psf_size**2., self.psf_size**2., 0.0) 

41 

42 # We want a set of point sources at various flux levels. 

43 self.dataset.addSource(10000.0, lsst.geom.Point2D(50.1, 49.8)) 

44 self.dataset.addSource(20000.0, lsst.geom.Point2D(100.5, 100.4)) 

45 self.dataset.addSource(30000.0, lsst.geom.Point2D(150.4, 149.6)) 

46 self.dataset.addSource(40000.0, lsst.geom.Point2D(200.2, 200.3)) 

47 self.dataset.addSource(50000.0, lsst.geom.Point2D(250.3, 250.1)) 

48 self.dataset.addSource(60000.0, lsst.geom.Point2D(300.4, 300.2)) 

49 self.dataset.addSource(70000.0, lsst.geom.Point2D(350.5, 350.6)) 

50 self.dataset.addSource(80000.0, lsst.geom.Point2D(400.6, 400.0)) 

51 self.dataset.addSource(90000.0, lsst.geom.Point2D(450.0, 450.0)) 

52 self.dataset.addSource(100000.0, lsst.geom.Point2D(500.7, 500.8)) 

53 

54 # Small test for Monte Carlo 

55 self.bbox_single = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), 

56 lsst.geom.Extent2I(101, 101)) 

57 self.dataset_single = lsst.meas.base.tests.TestDataset(self.bbox_single) 

58 self.dataset_single.psfShape = self.dataset.psfShape 

59 

60 self.dataset_single.addSource(100000.0, lsst.geom.Point2D(50.0, 50.0)) 

61 

62 self.all_apertures = (12, 14) 

63 self.all_scales = ((1.2, 1.7), (1.5, 2.0)) 

64 

65 def tearDown(self): 

66 del self.bbox 

67 del self.dataset 

68 del self.bbox_single 

69 del self.dataset_single 

70 

71 def makeAlgorithm(self, config=None): 

72 """Construct an algorithm and return both it and its schema. 

73 """ 

74 schema = lsst.meas.base.tests.TestDataset.makeMinimalSchema() 

75 if config is None: 

76 config = lsst.meas.base.SingleFrameCompensatedTophatFluxConfig() 

77 algorithm = lsst.meas.base.SingleFrameCompensatedTophatFluxPlugin( 

78 config, 

79 "base_CompensatedTophatFlux", 

80 schema, 

81 None 

82 ) 

83 return algorithm, schema 

84 

85 def testCompensatedTophatPlugin(self): 

86 """Test for correct instFlux given known position and shape. 

87 """ 

88 # In the z-band, HSC images have a noise of about 40.0 ADU, and a background 

89 # offset of ~ -0.6 ADU/pixel. This determines our test levels. 

90 for aperture in self.all_apertures: 

91 for inner_scale, outer_scale in self.all_scales: 

92 for dc_offset in (0.0, -1.0, 1.0): 

93 config = self.makeSingleFrameMeasurementConfig("base_CompensatedTophatFlux") 

94 config.algorithms["base_CompensatedTophatFlux"].apertures = [aperture] 

95 config.algorithms["base_CompensatedTophatFlux"].inner_scale = inner_scale 

96 config.algorithms["base_CompensatedTophatFlux"].outer_scale = outer_scale 

97 

98 task = self.makeSingleFrameMeasurementTask(config=config) 

99 exposure, catalog = self.dataset.realize(40.0, task.schema, randomSeed=0) 

100 exposure.image.array += dc_offset 

101 task.run(catalog, exposure) 

102 

103 filter_flux = catalog[f"base_CompensatedTophatFlux_{aperture}_instFlux"] 

104 filter_err = catalog[f"base_CompensatedTophatFlux_{aperture}_instFluxErr"] 

105 truth_flux = catalog["truth_instFlux"] 

106 

107 np.testing.assert_array_less(truth_flux, filter_flux + 3*filter_err) 

108 np.testing.assert_array_less(filter_flux - 3*filter_err, truth_flux) 

109 

110 if dc_offset == 0.0: 

111 # Use the no-offset run as a comparison for offset runs. 

112 flux_0 = filter_flux 

113 else: 

114 # Note: this tolerance is determined empirically, but this is 

115 # larger than preferable. 

116 self.assertFloatsAlmostEqual(filter_flux, flux_0, rtol=5e-3) 

117 

118 # The ratio of the filter flux to the truth flux should be consistent. 

119 # I'm not sure how to scale this with the error, so this is a loose 

120 # tolerance now. 

121 ratio = filter_flux / truth_flux 

122 self.assertLess(np.std(ratio), 0.04) 

123 

124 def testCompensatedTophatPluginFailure(self): 

125 """Test that the correct flag is set on failures.""" 

126 config = self.makeSingleFrameMeasurementConfig("base_CompensatedTophatFlux") 

127 config.algorithms["base_CompensatedTophatFlux"].apertures = [5, 15] 

128 config.algorithms["base_CompensatedTophatFlux"].inner_scale = 1.5 

129 config.algorithms["base_CompensatedTophatFlux"].outer_scale = 2.0 

130 

131 task = self.makeSingleFrameMeasurementTask(config=config) 

132 exposure, catalog = self.dataset.realize(40.0, task.schema, randomSeed=0) 

133 # Modify two objects to trigger the 2 failure conditions. 

134 catalog[0]["slot_Centroid_x"] = -20.0 

135 catalog[1]["slot_Centroid_x"] = -10.5 

136 task.run(catalog, exposure) 

137 self.assertTrue(catalog["base_CompensatedTophatFlux_5_flag"][0]) 

138 self.assertTrue(catalog["base_CompensatedTophatFlux_5_flag_bounds"][0]) 

139 self.assertTrue(catalog["base_CompensatedTophatFlux_15_flag"][0]) 

140 self.assertTrue(catalog["base_CompensatedTophatFlux_15_flag_bounds"][0]) 

141 self.assertTrue(catalog["base_CompensatedTophatFlux_5_flag"][1]) 

142 self.assertTrue(catalog["base_CompensatedTophatFlux_5_flag_bounds"][1]) 

143 

144 def testMonteCarlo(self): 

145 """Test an ideal simulation, with no noise. 

146 

147 Demonstrate that: 

148 

149 - We get exactly the right answer, and 

150 - The reported uncertainty agrees with a Monte Carlo test of the noise. 

151 """ 

152 nSamples = 500 

153 

154 for aperture in self.all_apertures: 

155 for inner_scale, outer_scale in self.all_scales: 

156 config = lsst.meas.base.SingleFrameCompensatedTophatFluxConfig() 

157 config.apertures = [aperture] 

158 config.inner_scale = inner_scale 

159 config.outer_scale = outer_scale 

160 

161 algorithm, schema = self.makeAlgorithm(config=config) 

162 

163 # Make a noiseless catalog. 

164 exposure, catalog = self.dataset_single.realize(1E-8, schema, randomSeed=1) 

165 

166 # Only use the high-flux source for the error tests. 

167 record = catalog[0] 

168 algorithm.measure(record, exposure) 

169 inst_flux = record[f"base_CompensatedTophatFlux_{aperture}_instFlux"] 

170 

171 for noise in (0.001, 0.01, 0.1): 

172 fluxes = np.zeros(nSamples) 

173 errs = np.zeros_like(fluxes) 

174 

175 for repeat in range(nSamples): 

176 # By using ``repeat`` to seed the RNG, we get results which 

177 # fall within the tolerances defined below. If we allow this 

178 # test to be truly random, passing becomes RNG-dependent. 

179 exposure_samp, catalog_samp = self.dataset_single.realize( 

180 noise*inst_flux, 

181 schema, 

182 randomSeed=repeat, 

183 ) 

184 record_samp = catalog_samp[0] 

185 algorithm.measure(record_samp, exposure_samp) 

186 fluxes[repeat] = record_samp[f"base_CompensatedTophatFlux_{aperture}_instFlux"] 

187 errs[repeat] = record_samp[f"base_CompensatedTophatFlux_{aperture}_instFluxErr"] 

188 

189 err_mean = np.mean(errs) 

190 flux_std = np.std(fluxes) 

191 self.assertFloatsAlmostEqual(err_mean, flux_std, rtol=0.03) 

192 

193 

194class TestMemory(lsst.utils.tests.MemoryTestCase): 

195 pass 

196 

197 

198def setup_module(module): 

199 lsst.utils.tests.init() 

200 

201 

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

203 lsst.utils.tests.init() 

204 unittest.main()