Coverage for tests/test_association_task.py: 98%

1# This file is part of ap_association. 

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 numpy as np 

23import pandas as pd 

24import unittest 

25 

26import lsst.geom as geom 

27import lsst.utils.tests 

28from lsst.ap.association import AssociationTask 

29from lsst.ap.association.association import AssociationConfig 

30 

31 

32class TestAssociationTask(unittest.TestCase): 

33 

34 def setUp(self): 

35 """Create sets of diaSources and diaObjects. 

36 """ 

37 rng = np.random.default_rng(1234) 

38 self.nObjects = 5 

39 scatter = 0.1/3600 

40 self.diaObjects = pd.DataFrame(data=[ 

41 {"ra": 0.04*(idx + 1), "dec": 0.04*(idx + 1), 

42 "diaObjectId": idx + 1} 

43 for idx in range(self.nObjects)]) 

44 self.diaObjects.set_index("diaObjectId", drop=False, inplace=True) 

45 self.nSources = 5 

46 self.diaSources = pd.DataFrame(data=[ 

47 {"ra": 0.04*idx + scatter*rng.uniform(-1, 1), 

48 "dec": 0.04*idx + scatter*rng.uniform(-1, 1), 

49 "diaSourceId": idx + 1 + self.nObjects, "diaObjectId": 0, "trailLength": 5.5*idx, 

50 "flags": 0} 

51 for idx in range(self.nSources)]) 

52 self.diaSourceZeroScatter = pd.DataFrame(data=[ 

53 {"ra": 0.04*idx, 

54 "dec": 0.04*idx, 

55 "diaSourceId": idx + 1 + self.nObjects, "diaObjectId": 0, "trailLength": 5.5*idx, 

56 "flags": 0} 

57 for idx in range(self.nSources)]) 

58 

59 def test_run(self): 

60 """Test the full task by associating a set of diaSources to 

61 existing diaObjects. 

62 """ 

63 config = AssociationTask.ConfigClass() 

64 assocTask = AssociationTask(config=config) 

65 results = assocTask.run(self.diaSources, self.diaObjects) 

66 

67 self.assertEqual(results.nUpdatedDiaObjects, len(self.diaObjects) - 1) 

68 self.assertEqual(results.nUnassociatedDiaObjects, 1) 

69 self.assertEqual(len(results.matchedDiaSources), 

70 len(self.diaObjects) - 1) 

71 self.assertEqual(len(results.unAssocDiaSources), 1) 

72 np.testing.assert_array_equal(results.matchedDiaSources["diaObjectId"].values, [1, 2, 3, 4]) 

73 np.testing.assert_array_equal(results.unAssocDiaSources["diaObjectId"].values, [0]) 

74 

75 def test_run_no_existing_objects(self): 

76 """Test the run method with a completely empty database. 

77 """ 

78 assocTask = AssociationTask() 

79 results = assocTask.run( 

80 self.diaSources, 

81 pd.DataFrame(columns=["ra", "dec", "diaObjectId", "trailLength"])) 

82 self.assertEqual(results.nUpdatedDiaObjects, 0) 

83 self.assertEqual(results.nUnassociatedDiaObjects, 0) 

84 self.assertEqual(len(results.matchedDiaSources), 0) 

85 self.assertTrue(np.all(results.unAssocDiaSources["diaObjectId"] == 0)) 

86 

87 def test_associate_sources(self): 

88 """Test the performance of the associate_sources method in 

89 AssociationTask. 

90 """ 

91 assoc_task = AssociationTask() 

92 assoc_result = assoc_task.associate_sources( 

93 self.diaObjects, self.diaSources) 

94 

95 for test_obj_id, expected_obj_id in zip( 

96 assoc_result.diaSources["diaObjectId"].to_numpy(), 

97 [0, 1, 2, 3, 4]): 

98 self.assertEqual(test_obj_id, expected_obj_id) 

99 np.testing.assert_array_equal(assoc_result.diaSources["diaObjectId"].values, [0, 1, 2, 3, 4]) 

100 

101 def test_score_and_match(self): 

102 """Test association between a set of sources and an existing 

103 DIAObjectCollection. 

104 """ 

105 

106 assoc_task = AssociationTask() 

107 score_struct = assoc_task.score(self.diaObjects, 

108 self.diaSourceZeroScatter, 

109 1.0 * geom.arcseconds) 

110 # Source 0 has no nearby DIAObject (closest is ~144" away). 

111 self.assertNotIn(0, score_struct.src_idx) 

112 # Sources 1..4 each have at least one candidate at ~zero score. 

113 for src_idx in range(1, len(self.diaSources)): 

114 mask = score_struct.src_idx == src_idx 

115 self.assertTrue(np.any(mask)) 

116 self.assertAlmostEqual(score_struct.scores[mask].min(), 0.0, 

117 places=10) 

118 

119 # Linear-assignment match: 4 sources match 4 distinct objects; 

120 # the fifth object is left unassociated. 

121 match_result = assoc_task.match( 

122 self.diaObjects, self.diaSources, score_struct) 

123 self.assertEqual(match_result.nUpdatedDiaObjects, 4) 

124 self.assertEqual(match_result.nUnassociatedDiaObjects, 1) 

125 

126 def test_remove_nan_dia_sources(self): 

127 """Test removing DiaSources with NaN locations. 

128 """ 

129 self.diaSources.loc[2, "ra"] = np.nan 

130 self.diaSources.loc[3, "dec"] = np.nan 

131 self.diaSources.loc[4, "ra"] = np.nan 

132 self.diaSources.loc[4, "dec"] = np.nan 

133 assoc_task = AssociationTask() 

134 out_dia_sources = assoc_task.check_dia_source_radec(self.diaSources) 

135 self.assertEqual(len(out_dia_sources), len(self.diaSources) - 3) 

136 

137 def test_score_falls_back_without_error_columns(self): 

138 """Score uses chord distance when raErr/decErr columns are absent. 

139 """ 

140 task = AssociationTask() 

141 result = task.score(self.diaObjects, self.diaSources, 

142 1.0 * geom.arcseconds) 

143 # 4 of 5 sources have at least one candidate. 

144 self.assertEqual(len(np.unique(result.src_idx)), 4) 

145 # Chord distance on the unit sphere — values are radians. 

146 self.assertTrue(np.all(result.scores < 1e-3)) 

147 

148 def test_score_falls_back_when_errors_all_nan(self): 

149 """Empty (all-NaN) raErr/decErr columns trigger the chord fallback. 

150 """ 

151 objects = self.diaObjects.copy() 

152 objects["raErr"] = np.nan 

153 objects["decErr"] = np.nan 

154 sources = self.diaSources.copy() 

155 sources["raErr"] = np.nan 

156 sources["decErr"] = np.nan 

157 task = AssociationTask() 

158 result = task.score(objects, sources, 1.0 * geom.arcseconds) 

159 self.assertEqual(len(np.unique(result.src_idx)), 4) 

160 self.assertTrue(np.all(result.scores < 1e-3)) 

161 

162 def test_chi2_accepts_within_uncertainty(self): 

163 """A 0.05" separation with 0.05" per-axis sigma yields chi^2 ~ 0.5. 

164 """ 

165 sig_deg = 0.05 / 3600.0 

166 objects = pd.DataFrame([ 

167 {"ra": 1.0, "dec": 1.0, "raErr": sig_deg, "decErr": sig_deg, 

168 "diaObjectId": 1} 

169 ]).set_index("diaObjectId", drop=False) 

170 sources = pd.DataFrame([ 

171 {"ra": 1.0, "dec": 1.0 + 0.05 / 3600.0, 

172 "raErr": sig_deg, "decErr": sig_deg, 

173 "diaSourceId": 100, "diaObjectId": 0} 

174 ]) 

175 task = AssociationTask() 

176 result = task.score(objects, sources, 1.0 * geom.arcseconds) 

177 self.assertEqual(len(result.src_idx), 1) 

178 self.assertEqual(result.src_idx[0], 0) 

179 self.assertEqual(result.obj_idx[0], 0) 

180 # var per axis = 2 * sig^2 (above the floor); dra=0, ddec=0.05". 

181 # chi^2 = (0.05)^2 / (2 * 0.05^2) = 0.5, and 

182 # NLL = 0.5*chi^2 + 0.5*ln(var_ra * var_dec). 

183 var = 2 * sig_deg ** 2 

184 expected_nll = 0.5 * 0.5 + 0.5 * np.log(var * var) 

185 self.assertAlmostEqual(result.scores[0], expected_nll, places=6) 

186 

187 def test_within_maxdist_matches_despite_high_chi2(self): 

188 """A 0.5" separation with 0.05" per-axis sigma (chi^2 ~ 50, i.e. 

189 ~7 sigma) is still matched: the chi^2 only ranks candidates, and 

190 ``maxDistArcSeconds`` is the sole association gate. 

191 """ 

192 sig_deg = 0.05 / 3600.0 

193 objects = pd.DataFrame([ 

194 {"ra": 1.0, "dec": 1.0, "raErr": sig_deg, "decErr": sig_deg, 

195 "diaObjectId": 1} 

196 ]).set_index("diaObjectId", drop=False) 

197 sources = pd.DataFrame([ 

198 {"ra": 1.0, "dec": 1.0 + sig_deg*10, # 10 sigma in dec 

199 "raErr": sig_deg, "decErr": sig_deg, 

200 "diaSourceId": 100, "diaObjectId": 0} 

201 ]) 

202 task = AssociationTask() 

203 result = task.score(objects, sources, 1.0 * geom.arcseconds) 

204 # The pair is kept as a candidate despite the large offset 

205 # (chi^2 ~ 50) since there is no significance cut. 

206 self.assertEqual(len(result.src_idx), 1) 

207 self.assertTrue(np.isfinite(result.scores[0])) 

208 # End-to-end the source is associated to the object. 

209 run_result = task.run(sources, objects) 

210 self.assertEqual(run_result.nUpdatedDiaObjects, 1) 

211 matched = run_result.matchedDiaSources.set_index("diaSourceId") 

212 self.assertEqual(int(matched.loc[100, "diaObjectId"]), 1) 

213 

214 def test_nll_prefers_better_localized_object(self): 

215 """With two candidate objects, the NLL cost associates to the 

216 better-localized (tighter) object even though a bare chi^2 would 

217 pick the looser, farther one, whose larger variance deflates its 

218 chi^2. 

219 """ 

220 src_sig = 0.05 / 3600.0 

221 objects = pd.DataFrame([ 

222 {"ra": 1.0, "dec": 1.0 + 0.15 / 3600.0, # tight and close 

223 "raErr": 0.05 / 3600.0, "decErr": 0.05 / 3600.0, 

224 "diaObjectId": 1}, 

225 {"ra": 1.0, "dec": 1.0 + 0.30 / 3600.0, # loose and far 

226 "raErr": 0.5 / 3600.0, "decErr": 0.5 / 3600.0, 

227 "diaObjectId": 2}, 

228 ]).set_index("diaObjectId", drop=False) 

229 sources = pd.DataFrame([ 

230 {"ra": 1.0, "dec": 1.0, 

231 "raErr": src_sig, "decErr": src_sig, 

232 "diaSourceId": 100, "diaObjectId": 0} 

233 ]) 

234 # A bare chi^2 would prefer the loose object 2 (smaller chi^2)... 

235 chi2_tight = (0.15 / 3600.0) ** 2 / (2 * (0.05 / 3600.0) ** 2) 

236 chi2_loose = (0.30 / 3600.0) ** 2 / ((0.05 / 3600.0) ** 2 

237 + (0.5 / 3600.0) ** 2) 

238 self.assertLess(chi2_loose, chi2_tight) 

239 

240 task = AssociationTask() 

241 score_struct = task.score(objects, sources, 1.0 * geom.arcseconds) 

242 # ...but the NLL cost prefers the tight object 1. obj_idx are 

243 # positional: 0 -> diaObjectId 1 (tight), 1 -> diaObjectId 2. 

244 nll = {int(obj): float(score) for obj, score 

245 in zip(score_struct.obj_idx, score_struct.scores)} 

246 self.assertLess(nll[0], nll[1]) 

247 # End-to-end the source associates to the tighter object 1. 

248 result = task.run(sources, objects) 

249 matched = result.matchedDiaSources.set_index("diaSourceId") 

250 self.assertEqual(int(matched.loc[100, "diaObjectId"]), 1) 

251 

252 def test_chi2_ra_wraparound(self): 

253 """Pairs straddling RA=0/360 are scored correctly. 

254 """ 

255 sig_deg = 0.05 / 3600.0 

256 objects = pd.DataFrame([ 

257 {"ra": 359.99999, "dec": 0.0, "raErr": sig_deg, "decErr": sig_deg, 

258 "diaObjectId": 1} 

259 ]).set_index("diaObjectId", drop=False) 

260 sources = pd.DataFrame([ 

261 {"ra": 0.00001, "dec": 0.0, 

262 "raErr": sig_deg, "decErr": sig_deg, 

263 "diaSourceId": 100, "diaObjectId": 0} 

264 ]) 

265 task = AssociationTask() 

266 result = task.score(objects, sources, 1.0 * geom.arcseconds) 

267 # Separation is 0.02" of RA across the wrap — within both cuts. 

268 self.assertEqual(len(result.src_idx), 1) 

269 self.assertEqual(result.src_idx[0], 0) 

270 self.assertEqual(result.obj_idx[0], 0) 

271 

272 def test_lap_recovers_match_lost_to_greedy(self): 

273 """LAP keeps a source matched when its nearest object is taken 

274 by a better-fitting competitor — the case that single-nearest 

275 greedy used to drop into a new DIAObject. 

276 """ 

277 sig_deg = 0.5 / 3600.0 # 0.5" per axis 

278 objects = pd.DataFrame([ 

279 {"ra": 1.0, "dec": 1.0, 

280 "raErr": sig_deg, "decErr": sig_deg, "diaObjectId": 1}, 

281 {"ra": 1.0, "dec": 1.0 + 0.99 / 3600.0, 

282 "raErr": sig_deg, "decErr": sig_deg, "diaObjectId": 2}, 

283 ]).set_index("diaObjectId", drop=False) 

284 # Source 100: closer to obj 1, but obj 2 is also a viable candidate. 

285 # Source 101: sits exactly on obj 1. 

286 sources = pd.DataFrame([ 

287 {"ra": 1.0, "dec": 1.0 + 0.4 / 3600.0, 

288 "raErr": sig_deg, "decErr": sig_deg, "diaSourceId": 100, 

289 "diaObjectId": 0}, 

290 {"ra": 1.0, "dec": 1.0, 

291 "raErr": sig_deg, "decErr": sig_deg, "diaSourceId": 101, 

292 "diaObjectId": 0}, 

293 ]) 

294 task = AssociationTask() 

295 result = task.run(sources, objects) 

296 # Both sources matched; LAP picks src 100 → obj 2, src 101 → obj 1. 

297 self.assertEqual(result.nUpdatedDiaObjects, 2) 

298 self.assertEqual(result.nUnassociatedDiaObjects, 0) 

299 matched = result.matchedDiaSources.set_index("diaSourceId") 

300 self.assertEqual(int(matched.loc[100, "diaObjectId"]), 2) 

301 self.assertEqual(int(matched.loc[101, "diaObjectId"]), 1) 

302 

303 def test_match_preserves_diaObjectId_dtype(self): 

304 """match() preserves the incoming diaObjectId dtype rather than 

305 overwriting it with uint64. This keeps pd.concat dtype-stable 

306 downstream — mixing uint64 and int64 silently promotes to 

307 float64. 

308 """ 

309 self.assertEqual(self.diaSources["diaObjectId"].dtype, np.int64) 

310 task = AssociationTask() 

311 result = task.run(self.diaSources, self.diaObjects) 

312 self.assertEqual( 

313 result.matchedDiaSources["diaObjectId"].dtype, np.int64) 

314 self.assertEqual( 

315 result.unAssocDiaSources["diaObjectId"].dtype, np.int64) 

316 

317 sources_u64 = self.diaSources.copy() 

318 sources_u64["diaObjectId"] = sources_u64["diaObjectId"].astype( 

319 np.uint64) 

320 result = task.run(sources_u64, self.diaObjects) 

321 self.assertEqual( 

322 result.matchedDiaSources["diaObjectId"].dtype, np.uint64) 

323 self.assertEqual( 

324 result.unAssocDiaSources["diaObjectId"].dtype, np.uint64) 

325 with self.assertRaises(ValueError): 

326 sources_float = self.diaSources.copy() 

327 sources_float["diaObjectId"] = sources_float["diaObjectId"].astype(float) 

328 task.run(sources_float, self.diaObjects) 

329 

330 def test_fallback_sigma_scores_missing_error_pair(self): 

331 """A pair whose errors are missing is scored using 

332 ``fallbackSigmaArcSeconds``: the substituted per-axis uncertainty 

333 sets the pair's score. The fallback affects only the score 

334 (ranking), not whether the pair is a candidate. 

335 """ 

336 sig_deg = 0.05 / 3600.0 

337 # One row in each catalog carries a real error so chi^2 mode 

338 # triggers; the rows under test (id 100 / obj 1) have NaN errors. 

339 objects = pd.DataFrame([ 

340 {"ra": 1.0, "dec": 1.0, 

341 "raErr": np.nan, "decErr": np.nan, "diaObjectId": 1}, 

342 {"ra": 2.0, "dec": 2.0, 

343 "raErr": sig_deg, "decErr": sig_deg, "diaObjectId": 2}, 

344 ]).set_index("diaObjectId", drop=False) 

345 sources = pd.DataFrame([ 

346 {"ra": 1.0, "dec": 1.0 + 0.5 / 3600.0, 

347 "raErr": np.nan, "decErr": np.nan, 

348 "diaSourceId": 100, "diaObjectId": 0}, 

349 {"ra": 2.0, "dec": 2.0, 

350 "raErr": sig_deg, "decErr": sig_deg, 

351 "diaSourceId": 200, "diaObjectId": 0}, 

352 ]) 

353 

354 def missing_pair_score(fallback): 

355 cfg = AssociationConfig() 

356 cfg.fallbackSigmaArcSeconds = fallback 

357 result = AssociationTask(config=cfg).score( 

358 objects, sources, 1.0 * geom.arcseconds) 

359 # The missing-error pair is source row 0 -> object row 0. 

360 mask = (result.src_idx == 0) & (result.obj_idx == 0) 

361 self.assertEqual(int(mask.sum()), 1) 

362 return float(result.scores[mask][0]) 

363 

364 def expected_nll(fallback, chi2): 

365 # Both rows have missing errors, so the combined per-axis 

366 # variance is 2 * fallback^2 (deg^2), above the 0.05" floor. 

367 var = 2 * (fallback / 3600.0) ** 2 

368 return 0.5 * chi2 + 0.5 * np.log(var * var) 

369 

370 # ddec = 0.5"; fallback 0.05" -> chi^2 = 50, fallback 0.2" -> 3.125. 

371 self.assertAlmostEqual( 

372 missing_pair_score(0.05), expected_nll(0.05, 50.0), places=6) 

373 self.assertAlmostEqual( 

374 missing_pair_score(0.2), expected_nll(0.2, 3.125), places=6) 

375 

376 # Even with the small fallback (large chi^2) the pair is within 

377 # maxDist, so it is still matched to its object. 

378 cfg = AssociationConfig() 

379 cfg.fallbackSigmaArcSeconds = 0.05 

380 result = AssociationTask(config=cfg).run(sources, objects) 

381 matched = result.matchedDiaSources.set_index("diaSourceId") 

382 self.assertEqual(int(matched.loc[100, "diaObjectId"]), 1) 

383 self.assertEqual(int(matched.loc[200, "diaObjectId"]), 2) 

384 

385 def test_contention_leaves_extra_source_unassociated(self): 

386 """When two sources fall within ``maxDistArcSeconds`` of a single 

387 object, the better-fitting one is matched and the other is left 

388 unassociated rather than forced onto the same object. 

389 """ 

390 sig_deg = 0.1 / 3600.0 

391 objects = pd.DataFrame([ 

392 {"ra": 1.0, "dec": 1.0, 

393 "raErr": sig_deg, "decErr": sig_deg, "diaObjectId": 1}, 

394 ]).set_index("diaObjectId", drop=False) 

395 sources = pd.DataFrame([ 

396 {"ra": 1.0, "dec": 1.0 + 0.2 / 3600.0, # closer to the object 

397 "raErr": sig_deg, "decErr": sig_deg, "diaSourceId": 100, 

398 "diaObjectId": 0}, 

399 {"ra": 1.0, "dec": 1.0 + 0.5 / 3600.0, # farther 

400 "raErr": sig_deg, "decErr": sig_deg, "diaSourceId": 101, 

401 "diaObjectId": 0}, 

402 ]) 

403 task = AssociationTask() 

404 result = task.run(sources, objects) 

405 # One object, so at most one source can match it. 

406 self.assertEqual(result.nUpdatedDiaObjects, 1) 

407 self.assertEqual(result.nUnassociatedDiaObjects, 0) 

408 self.assertEqual(len(result.matchedDiaSources), 1) 

409 self.assertEqual(len(result.unAssocDiaSources), 1) 

410 # The closer source wins the object; the other is unassociated. 

411 matched = result.matchedDiaSources.set_index("diaSourceId") 

412 self.assertEqual(int(matched.loc[100, "diaObjectId"]), 1) 

413 self.assertEqual( 

414 int(result.unAssocDiaSources["diaSourceId"].iloc[0]), 101) 

415 

416 

417class MemoryTester(lsst.utils.tests.MemoryTestCase): 

418 pass 

419 

420 

421def setup_module(module): 

422 lsst.utils.tests.init() 

423 

424 

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

426 lsst.utils.tests.init() 

427 unittest.main()