Coverage for tests/test_compareDiaSources.py: 98%
130 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-23 09:35 +0000
1# This file is part of analysis_ap.
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/>.
22import os
23import tempfile
24import unittest
25from unittest import mock
27import astropy.table
28import numpy as np
30import lsst.afw.table as afwTable
31import lsst.daf.butler.tests as butlerTests
32import lsst.geom
33import lsst.utils.tests
34from lsst.analysis.ap.compareDiaSources import compare_dia_sources, plot_dia_source_comparison
36RA0, DEC0 = 10.0, -5.0
37ARCSEC = 1/3600
40def _position(slot, offset=0.0):
41 """Sky position of a test source: slots are 10 arcsec apart in RA, and
42 ``offset`` moves the source north, in arcseconds.
43 """
44 return RA0 + 10*slot*ARCSEC/np.cos(np.radians(DEC0)), DEC0 + offset*ARCSEC
47def _table(rows):
48 """Make a diaSource table from (id, visit, detector, slot, offset,
49 isNegative, psfFlux, pixelFlags_bad) rows.
50 """
51 columns = list(zip(*rows))
52 ra, dec = zip(*(_position(slot, offset) for slot, offset in zip(columns[3], columns[4])))
53 n = len(rows)
54 return astropy.table.Table({
55 "diaSourceId": np.array(columns[0], dtype=np.int64),
56 "visit": np.array(columns[1], dtype=np.int64),
57 "detector": np.array(columns[2], dtype=np.int64),
58 "ra": np.array(ra),
59 "dec": np.array(dec),
60 "x": np.full(n, 50.0) + 5*np.array(columns[4]),
61 "y": np.full(n, 60.0),
62 "isNegative": np.array(columns[5], dtype=bool),
63 "psfFlux": np.array(columns[6], dtype=float),
64 "psfFluxErr": np.full(n, 10.0),
65 "snr": np.array(columns[6], dtype=float)/10,
66 "reliability": np.full(n, 0.9),
67 "pixelFlags_bad": np.array(columns[7], dtype=bool),
68 "band": np.full(n, "g"),
69 })
72def _catalog(sources):
73 """Make a catalog of detections from (id, slot, reliability) entries."""
74 schema = afwTable.SourceTable.makeMinimalSchema()
75 negative = schema.addField("is_negative", type="Flag", doc="Negative detection.")
76 reliability = schema.addField("reliability", type="F", doc="Reliability score.")
77 catalog = afwTable.SourceCatalog(schema)
78 for sourceId, slot, score in sources:
79 record = catalog.addNew()
80 record.setId(sourceId)
81 record.setCoord(lsst.geom.SpherePoint(*_position(slot), lsst.geom.degrees))
82 record.set(negative, False)
83 record.set(reliability, score)
84 return catalog
87def _make_repo(path):
88 butler = butlerTests.makeTestRepo(path, {"instrument": ["notACam"], "visit": [1, 2],
89 "detector": [1, 2]})
90 butler.registry.registerRun("run")
91 butlerTests.addDatasetType(butler, "dia_source_visit", {"instrument", "visit"}, "ArrowAstropy")
92 butlerTests.addDatasetType(butler, "dia_source_detector", {"instrument", "visit", "detector"},
93 "ArrowAstropy")
94 for name in ("dia_source_low_reliability", "rejected_dia_source", "dia_source_unfiltered"):
95 butlerTests.addDatasetType(butler, name, {"instrument", "visit", "detector"}, "SourceCatalog")
96 butlerTests.addDatasetType(butler, "long_trailed_source_detector", {"instrument", "visit", "detector"},
97 "ArrowAstropy")
98 return butler.clone(collections="run", run="run")
101def _write_runs(root):
102 """Write the reference and new runs and return their butlers.
104 Visit 1, matched pairs: 101-201 (0.05 arcsec apart), 102-202 (on
105 different detectors, with a flux and flag change), 103-203 (of opposite
106 signs) and 104-204. 105 is a second detection of the 104/204 source, so
107 it is unmatched. 205-209 are new only, and the reference run's
108 catalogs explain why. Visit 2: the new run wrote only detector tables,
109 for detector 1, so 112 on detector 2 has no data in it.
110 """
111 ref = _make_repo(os.path.join(root, "ref"))
112 new = _make_repo(os.path.join(root, "new"))
113 put = dict(instrument="notACam")
114 ref.put(_table([(101, 1, 1, 0, 0.0, False, 100.0, False),
115 (102, 1, 1, 1, 0.0, False, 200.0, False),
116 (103, 1, 2, 2, 0.0, True, -50.0, False),
117 (104, 1, 1, 3, 0.0, False, 300.0, False),
118 (105, 1, 1, 3, 0.3, False, 30.0, False)]),
119 "dia_source_visit", visit=1, **put)
120 ref.put(_table([(111, 2, 1, 10, 0.0, False, 100.0, False),
121 (112, 2, 2, 11, 0.0, False, 100.0, False)]),
122 "dia_source_visit", visit=2, **put)
123 ref.put(_catalog([(901, 4, 0.2)]), "dia_source_low_reliability", visit=1, detector=1, **put)
124 ref.put(_catalog([(902, 5, np.nan)]), "rejected_dia_source", visit=1, detector=1, **put)
125 ra, dec = np.radians(_position(8))
126 ref.put(astropy.table.Table({"id": [903], "coord_ra": [ra], "coord_dec": [dec], "is_negative": [False]}),
127 "long_trailed_source_detector", visit=1, detector=1, **put)
128 # Removed detections keep their IDs from the detection catalog.
129 ref.put(_catalog([(910, 0, np.nan), (911, 1, np.nan), (913, 3, np.nan), (901, 4, np.nan),
130 (902, 5, np.nan), (906, 6, np.nan), (903, 8, np.nan)]),
131 "dia_source_unfiltered", visit=1, detector=1, **put)
133 new.put(_table([(201, 1, 1, 0, 0.05, False, 100.0, False),
134 (202, 1, 2, 1, 0.0, False, 230.0, True),
135 (203, 1, 2, 2, 0.0, False, 50.0, False),
136 (204, 1, 1, 3, 0.1, False, 300.0, False),
137 (205, 1, 1, 4, 0.0, False, 60.0, False),
138 (206, 1, 1, 5, 0.0, False, 60.0, False),
139 (207, 1, 1, 6, 0.0, False, 60.0, False),
140 (208, 1, 1, 7, 0.0, False, 60.0, False),
141 (209, 1, 1, 8, 0.0, False, 60.0, False)]),
142 "dia_source_visit", visit=1, **put)
143 new.put(_table([(211, 2, 1, 10, 0.0, False, 100.0, False)]),
144 "dia_source_detector", visit=2, detector=1, **put)
145 return ref, new
148class CompareDiaSourcesTestCase(lsst.utils.tests.TestCase):
149 @classmethod
150 def setUpClass(cls):
151 cls.root = tempfile.TemporaryDirectory()
152 cls.ref, cls.new = _write_runs(cls.root.name)
153 cls.result = compare_dia_sources(cls.ref, cls.new)
155 @classmethod
156 def tearDownClass(cls):
157 cls.root.cleanup()
159 def testSummary(self):
160 summary = self.result.summary.set_index("visit")
161 self.assertEqual(list(summary.index), [1, 2])
162 self.assertEqual(summary.loc[1, ["n_ref", "n_new", "n_matched", "n_ref_only", "n_new_only"]].tolist(),
163 [5, 9, 4, 1, 5])
164 self.assertEqual(summary.loc[2, ["n_ref", "n_new", "n_matched", "n_ref_only", "n_new_only"]].tolist(),
165 [2, 1, 1, 1, 0])
166 self.assertEqual(summary.loc[1, "table_new"], "dia_source_visit")
167 self.assertEqual(summary.loc[2, "table_new"], "dia_source_detector")
168 for visit in (1, 2):
169 expected = self.ref.registry.expandDataId(instrument="notACam", visit=visit)["band"]
170 self.assertEqual(summary.loc[visit, "band"], expected)
172 def testMatched(self):
173 matched = self.result.matched.set_index("diaSourceId_ref")
174 self.assertEqual(dict(matched["diaSourceId_new"]), {101: 201, 102: 202, 103: 203, 104: 204, 111: 211})
175 self.assertAlmostEqual(matched.loc[101, "separation_arcsec"], 0.05, places=6)
176 self.assertAlmostEqual(matched.loc[101, "dx"], 0.25)
177 self.assertTrue(np.isnan(matched.loc[102, "dx"]))
178 self.assertAlmostEqual(matched.loc[102, "delta_psfFlux_sigma"], 30/np.sqrt(200))
179 self.assertEqual(matched.loc[102, "flags_changed"], "pixelFlags_bad")
180 self.assertEqual(matched.loc[101, "flags_changed"], "")
181 self.assertEqual(matched.loc[103, "flags_changed"], "isNegative")
183 def testReasons(self):
184 ref_only = self.result.ref_only.set_index("diaSourceId")
185 self.assertEqual(dict(ref_only["reason"]),
186 {105: "matched_elsewhere", 112: "no_data"})
187 self.assertEqual(ref_only.loc[105, "other_diaSourceId"], 204)
188 new_only = self.result.new_only.set_index("diaSourceId")
189 self.assertEqual(dict(new_only["reason"]),
190 {205: "low_reliability", 206: "rejected",
191 207: "removed_later", 208: "not_detected", 209: "long_trail"})
192 self.assertEqual(new_only.loc[209, "other_diaSourceId"], 903)
193 self.assertEqual(new_only.loc[205, "other_diaSourceId"], 901)
194 self.assertAlmostEqual(new_only.loc[205, "other_reliability"], 0.2, places=6)
195 self.assertEqual(set(new_only["instrument"]), {"notACam"})
197 def testFlagChanges(self):
198 self.assertEqual(self.result.flag_changes.to_dict("records"),
199 [{"flag": "isNegative", "n_ref_only": 1, "n_new_only": 0},
200 {"flag": "pixelFlags_bad", "n_ref_only": 0, "n_new_only": 1}])
202 def testReasonCounts(self):
203 counts = self.result.reason_counts()
204 self.assertEqual(list(counts.index), ["matched_elsewhere", "low_reliability", "rejected",
205 "long_trail", "removed_later", "not_detected", "no_data"])
206 self.assertEqual(counts.loc["matched_elsewhere"].tolist(), [1, 0])
208 def testOptions(self):
209 result = compare_dia_sources(self.ref, self.new, visits=[2], explain_unmatched=False)
210 self.assertEqual(list(result.summary["visit"]), [2])
211 self.assertEqual(list(result.ref_only["reason"]), [None])
213 def testQueries(self):
214 """Catalogs are found with one query per dataset type, and detector
215 tables are only queried for visits without a visit table.
216 """
217 with mock.patch.object(self.ref, "query_datasets", wraps=self.ref.query_datasets) as query_ref, \
218 mock.patch.object(self.new, "query_datasets", wraps=self.new.query_datasets) as query_new, \
219 mock.patch.object(self.ref, "find_dataset", side_effect=AssertionError), \
220 mock.patch.object(self.new, "find_dataset", side_effect=AssertionError):
221 compare_dia_sources(self.ref, self.new, visits=[1, 2])
222 queried_ref = [call.args[0] for call in query_ref.call_args_list]
223 queried_new = [call.args[0] for call in query_new.call_args_list]
224 # Both runs have a visit table for visit 1, but only the reference
225 # run has one for visit 2.
226 self.assertNotIn("dia_source_detector", queried_ref)
227 self.assertEqual(queried_new.count("dia_source_detector"), 1)
228 self.assertEqual(query_new.call_args_list[queried_new.index("dia_source_detector")].kwargs["bind"],
229 {"visits": [2]})
230 for dataset in ("dia_source_low_reliability", "rejected_dia_source", "long_trailed_source_detector",
231 "dia_source_unfiltered"):
232 self.assertEqual(queried_ref.count(dataset), 1)
233 self.assertEqual(queried_new.count(dataset), 1)
235 def testPlot(self):
236 figure = plot_dia_source_comparison(self.result, labels=("main", "ticket"))
237 self.assertEqual(len(figure.axes), 4)
239 def testNoSources(self):
240 with tempfile.TemporaryDirectory() as root:
241 empty = _make_repo(os.path.join(root, "empty"))
242 with self.assertRaises(LookupError):
243 compare_dia_sources(empty, empty)
246class MemoryTester(lsst.utils.tests.MemoryTestCase):
247 pass
250def setup_module(module):
251 lsst.utils.tests.init()
254if __name__ == "__main__": 254 ↛ 255line 254 didn't jump to line 255 because the condition on line 254 was never true
255 lsst.utils.tests.init()
256 unittest.main()