Coverage for tests/test_matchSourceInjected.py: 98%

1# 

2# This file is part of ap_pipe. 

3# 

4# Developed for the LSST Data Management System. 

5# This product includes software developed by the LSST Project 

6# (http://www.lsst.org). 

7# See the COPYRIGHT file at the top-level directory of this distribution 

8# for details of code ownership. 

9# 

10# This program is free software: you can redistribute it and/or modify 

11# it under the terms of the GNU General Public License as published by 

12# the Free Software Foundation, either version 3 of the License, or 

13# (at your option) any later version. 

14# 

15# This program is distributed in the hope that it will be useful, 

16# but WITHOUT ANY WARRANTY; without even the implied warranty of 

17# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 

18# GNU General Public License for more details. 

19# 

20# You should have received a copy of the GNU General Public License 

21# along with this program. If not, see <http://www.gnu.org/licenses/>. 

22# 

23 

24import numpy as np 

25import unittest 

26 

27from astropy.table import Table, vstack 

28from lsst.afw.table import SourceCatalog, SourceTable 

29from lsst.afw.image import ExposureF 

30from lsst.afw.geom import makeSkyWcs 

31from lsst.geom import SpherePoint, degrees, Point2D, Extent2I 

32import lsst.utils.tests 

33 

34from lsst.pipe.tasks.matchDiffimSourceInjected import ( 

35 MatchInjectedToDiaSourceTask, 

36 MatchInjectedToAssocDiaSourceTask, 

37 MatchInjectedToDiaSourceConfig, 

38 MatchInjectedToAssocDiaSourceConfig, 

39) 

40 

41 

42class BaseTestMatchInjected(lsst.utils.tests.TestCase): 

43 def setUp(self): 

44 rng = np.random.RandomState(6) 

45 # 0.35 arcsec = 0.5/np.sqrt(2) arcsec (half of the diagonal of a 0.5 arcsec square) 

46 offsetFactor = 1./3600 * 0.35 

47 

48 # Create a mock injected catalog 

49 self.injectedCat = Table( 

50 { 

51 "injection_id": [1, 2, 3, 5, 6, 7, 12, 21, 31, 49], 

52 "injection_flag": np.repeat(0, 10), 

53 # random positions with 10 arcmin size 

54 # ra centered at 0, dec centered at -30 

55 "ra": (1/6.) * rng.random(size=10), 

56 "dec": -30 + (1/6.) * rng.random(size=10), 

57 "mag": 20.25 - 0.5 * rng.random(size=10), # random magnitudes 

58 "source_type": np.repeat("DeltaFunction", 10) 

59 } 

60 ) 

61 

62 # Create a mock diaSources catalog 

63 schema = SourceTable.makeMinimalSchema() 

64 self.diaSources = SourceCatalog(schema) 

65 for i in range(5): 

66 record = self.diaSources.addNew() 

67 record.setId(100 + i) 

68 record.setCoord( 

69 # Random posisions of diaSources 

70 SpherePoint((1/6.) * rng.random(), -30 + (1/6.) * rng.random(), degrees) 

71 ) 

72 

73 for i, (ra, dec) in enumerate(self.injectedCat[['ra', 'dec']][:8]): 

74 record = self.diaSources.addNew() 

75 sign = rng.choice([-1, 1], size=2) 

76 record.setId(i + 200) 

77 record.setCoord( 

78 SpherePoint( 

79 ra+sign[0]*rng.random()*offsetFactor, 

80 dec+sign[1]*rng.random()*offsetFactor, 

81 degrees 

82 ) 

83 ) 

84 

85 # Create a mock difference image 

86 self.diffIm = ExposureF(Extent2I(4096, 4096)) 

87 crpix = Point2D(0, 0) 

88 crval = SpherePoint(0, -30, degrees) 

89 cdMatrix = np.array([[5.19513851e-05, 2.81124812e-07], 

90 [3.25186974e-07, 5.19112119e-05]]) 

91 wcs = makeSkyWcs(crpix, crval, cdMatrix) 

92 self.diffIm.setWcs(wcs) 

93 

94 # add a fake source outside of the image 

95 self.injectedCatForTrimming = self.injectedCat.copy() 

96 self.injectedCatForTrimming.add_row( 

97 { 

98 'injection_id': 50, 

99 'injection_flag': 0, 

100 'ra': 360. - 0.1/6., 

101 'dec': -30 - 0.1/6., 

102 'mag': 20.0, 

103 'source_type': 'DeltaFunction' 

104 } 

105 ) 

106 

107 # Variable sources as written by the catalog generator: parent and twin 

108 # are both flagged isVariable and link to each other through twin_id. 

109 # Row 0 is a parent with a matching diaSource, row 9 has none. 

110 variableCat = self.injectedCat.copy() 

111 variableCat["twin_id"] = np.zeros(len(variableCat), dtype=np.int64) 

112 variableCat["isVariable"] = np.zeros(len(variableCat), dtype=bool) 

113 variableCat["isVisitSource"] = np.ones(len(variableCat), dtype=bool) 

114 variableCat["isTemplateSource"] = np.zeros(len(variableCat), dtype=bool) 

115 parents = [0, 9] 

116 variableCat["isVariable"][parents] = True 

117 twins = variableCat[parents].copy() 

118 twins["injection_id"] = [100, 101] 

119 twins["twin_id"] = variableCat["injection_id"][parents] 

120 twins["isVisitSource"] = False 

121 twins["isTemplateSource"] = True 

122 variableCat["twin_id"][parents] = [100, 101] 

123 self.injectedCatVariable = vstack([variableCat, twins]) 

124 

125 # only 4 injected sources are associated 

126 self.assocDiaSources = Table( 

127 { 

128 "diaSourceId": [101, 102, 103, 201, 202, 205, 207], 

129 "band": np.repeat("r", 7), 

130 "visit": np.repeat(410001, 7), 

131 "detector": np.repeat(0, 7), 

132 "diaObjectId": np.arange(7), 

133 } 

134 ) 

135 

136 

137class TestMatchInjectedToDiaSourceTask(BaseTestMatchInjected): 

138 

139 def test_run(self): 

140 config = MatchInjectedToDiaSourceConfig() 

141 config.matchDistanceArcseconds = 0.5 

142 config.doMatchVisit = False 

143 config.doForcedMeasurement = False 

144 

145 task = MatchInjectedToDiaSourceTask(config=config) 

146 

147 result = task.run(self.injectedCat, self.diffIm, self.diaSources) 

148 self.assertEqual(len(result.matchDiaSources), len(self.injectedCat)) 

149 self.assertEqual(np.sum(result.matchDiaSources['diaSourceId'] > 0), 8) 

150 self.assertEqual(np.sum(result.matchDiaSources['dist_diaSrc'] > 0), 8) 

151 self.assertEqual( 

152 np.sum(np.abs(result.matchDiaSources['dist_diaSrc']) < config.matchDistanceArcseconds), 8 

153 ) 

154 

155 def test_run_trimming(self): 

156 config = MatchInjectedToDiaSourceConfig() 

157 config.matchDistanceArcseconds = 0.5 

158 config.doMatchVisit = True 

159 config.doForcedMeasurement = False 

160 

161 task = MatchInjectedToDiaSourceTask(config=config) 

162 result = task.run(self.injectedCatForTrimming, self.diffIm, self.diaSources) 

163 

164 self.assertEqual(len(result.matchDiaSources), len(self.injectedCatForTrimming) - 1) 

165 self.assertEqual(np.sum(result.matchDiaSources['diaSourceId'] > 0), 8) 

166 self.assertEqual(np.sum(result.matchDiaSources['dist_diaSrc'] > 0), 8) 

167 self.assertEqual( 

168 np.sum(np.abs(result.matchDiaSources['dist_diaSrc']) < config.matchDistanceArcseconds), 8 

169 ) 

170 

171 def test_splitVariables(self): 

172 config = MatchInjectedToDiaSourceConfig() 

173 task = MatchInjectedToDiaSourceTask(config=config) 

174 

175 initial, doubles = task._splitVariables(self.injectedCatVariable) 

176 

177 # Parents stay with the static sources, only the twins are split out. 

178 self.assertEqual(len(initial), 10) 

179 self.assertEqual(sorted(doubles["injection_id"]), [100, 101]) 

180 self.assertIn(1, initial["injection_id"]) 

181 self.assertIn(49, initial["injection_id"]) 

182 

183 def test_run_variables(self): 

184 config = MatchInjectedToDiaSourceConfig() 

185 config.matchDistanceArcseconds = 0.5 

186 config.doMatchVisit = False 

187 config.doForcedMeasurement = False 

188 

189 task = MatchInjectedToDiaSourceTask(config=config) 

190 result = task.run(self.injectedCatVariable, self.diffIm, self.diaSources) 

191 matched = result.matchDiaSources 

192 

193 self.assertEqual(len(matched), len(self.injectedCatVariable)) 

194 byId = {row["injection_id"]: row for row in matched} 

195 

196 # The twin of a matched parent inherits its diaSourceId. 

197 self.assertGreater(byId[1]["diaSourceId"], 0) 

198 self.assertEqual(byId[100]["diaSourceId"], byId[1]["diaSourceId"]) 

199 self.assertEqual(byId[100]["dist_diaSrc"], byId[1]["dist_diaSrc"]) 

200 # The twin of an unmatched parent stays unmatched. 

201 self.assertEqual(byId[49]["diaSourceId"], 0) 

202 self.assertEqual(byId[101]["diaSourceId"], 0) 

203 self.assertEqual(byId[101]["dist_diaSrc"], -1) 

204 # 8 static matches plus the one matched twin. 

205 self.assertEqual(np.sum(matched["diaSourceId"] > 0), 9) 

206 

207 def test_getVectors(self): 

208 config = MatchInjectedToDiaSourceConfig() 

209 config.matchDistanceArcseconds = 0.5 

210 config.doMatchVisit = False 

211 config.doForcedMeasurement = False 

212 

213 task = MatchInjectedToDiaSourceTask(config=config) 

214 

215 ras = np.radians(self.injectedCat['ra']) 

216 decs = np.radians(self.injectedCat['dec']) 

217 vectors = task._getVectors(ras, decs) 

218 self.assertEqual(vectors.shape, (10, 3)) 

219 

220 

221class TestMatchInjectedToAssocDiaSourceTask(BaseTestMatchInjected): 

222 

223 def test_run(self): 

224 config = MatchInjectedToDiaSourceConfig() 

225 config.matchDistanceArcseconds = 0.5 

226 config.doMatchVisit = False 

227 config.doForcedMeasurement = False 

228 

229 task = MatchInjectedToDiaSourceTask(config=config) 

230 

231 result = task.run(self.injectedCat, self.diffIm, self.diaSources) 

232 

233 configAssoc = MatchInjectedToAssocDiaSourceConfig() 

234 taskAssoc = MatchInjectedToAssocDiaSourceTask(config=configAssoc) 

235 resultAssoc = taskAssoc.run(self.assocDiaSources, result.matchDiaSources) 

236 self.assertEqual(len(resultAssoc.matchAssocDiaSources), len(self.injectedCat)) 

237 self.assertEqual(np.sum(resultAssoc.matchAssocDiaSources['isAssocDiaSource']), 4) 

238 

239 

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

241 unittest.main()