Coverage for tests/test_compare.py: 97%
102 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 unittest
24import astropy.units as u
25import numpy as np
26import pandas as pd
28import lsst.utils.tests
29from lsst.analysis.ap.compare import (
30 match_catalogs, flux_residuals, match_to_truth,
31)
34def _make_frame(rows, **defaults):
35 """Build a small DiaSource-like DataFrame from a list of partial dicts."""
36 full = []
37 for row in rows:
38 merged = dict(defaults)
39 merged.update(row)
40 full.append(merged)
41 return pd.DataFrame(full)
44class TestMatchCatalogs(lsst.utils.tests.TestCase):
45 """Spatial cross-matching is the building block for compare_sources;
46 these tests pin down its behavior on small hand-crafted catalogs."""
48 def setUp(self):
49 # Three nearby points, two on (visit=1, detector=0), one on
50 # (visit=1, detector=1). All separations are tiny (~0.1 arcsec).
51 self.srcs1 = _make_frame([
52 {"diaSourceId": 1, "ra": 10.0, "dec": -5.0, "visit": 1, "detector": 0},
53 {"diaSourceId": 2, "ra": 10.001, "dec": -5.0, "visit": 1, "detector": 0},
54 {"diaSourceId": 3, "ra": 11.0, "dec": -5.0, "visit": 1, "detector": 1},
55 ])
56 # srcs2: id 11 matches id 1 (same coord). id 12 is far from id 2.
57 # No row in srcs2 for (visit=1, detector=1), so id 3 must be unique1.
58 self.srcs2 = _make_frame([
59 {"diaSourceId": 11, "ra": 10.0, "dec": -5.0, "visit": 1, "detector": 0},
60 {"diaSourceId": 12, "ra": 10.5, "dec": -5.0, "visit": 1, "detector": 0},
61 ])
63 def test_basic_match(self):
64 matched, unique1, unique2 = match_catalogs(
65 self.srcs1, self.srcs2, radius=1*u.arcsec)
66 # id=1 in srcs1 matches id=11 in srcs2.
67 self.assertEqual(list(matched["diaSourceId"]), [1])
68 self.assertEqual(list(matched["diaSourceId_2"]), [11])
69 # The match is essentially zero arcsec apart.
70 self.assertLess(float(matched["xmatch_dist_arcsec"].iloc[0]), 1e-6)
71 # ids 2, 3 from srcs1 are unique; ids 12 from srcs2 is unique.
72 self.assertEqual(set(unique1["diaSourceId"]), {2, 3})
73 self.assertEqual(set(unique2["diaSourceId"]), {12})
75 def test_radius_units(self):
76 # Bare-float radius is interpreted as arcseconds.
77 matched_q, _, _ = match_catalogs(self.srcs1, self.srcs2,
78 radius=1*u.arcsec)
79 matched_f, _, _ = match_catalogs(self.srcs1, self.srcs2, radius=1.0)
80 self.assertEqual(len(matched_q), len(matched_f))
82 def test_no_grouping(self):
83 # With on=(), sources match across visit/detector boundaries.
84 srcs1 = _make_frame([
85 {"diaSourceId": 1, "ra": 10.0, "dec": 0.0, "visit": 1, "detector": 0},
86 ])
87 srcs2 = _make_frame([
88 {"diaSourceId": 2, "ra": 10.0, "dec": 0.0, "visit": 99, "detector": 99},
89 ])
90 # With grouping, no match (different visits/detectors).
91 matched, _, _ = match_catalogs(srcs1, srcs2, radius=1*u.arcsec)
92 self.assertEqual(len(matched), 0)
93 # Without grouping, the spatial match succeeds.
94 matched, _, _ = match_catalogs(srcs1, srcs2, radius=1*u.arcsec, on=())
95 self.assertEqual(len(matched), 1)
97 def test_empty_inputs(self):
98 empty = self.srcs1.iloc[0:0]
99 matched, u1, u2 = match_catalogs(empty, self.srcs2, radius=1*u.arcsec)
100 self.assertEqual(len(matched), 0)
101 self.assertEqual(len(u1), 0)
102 self.assertEqual(set(u2["diaSourceId"]), {11, 12})
104 matched, u1, u2 = match_catalogs(self.srcs1, empty, radius=1*u.arcsec)
105 self.assertEqual(len(matched), 0)
106 self.assertEqual(set(u1["diaSourceId"]), {1, 2, 3})
107 self.assertEqual(len(u2), 0)
110class TestOneToOneMatch(lsst.utils.tests.TestCase):
111 """Two sources in catalog 1 compete for the one source in catalog 2
112 that is nearest to both."""
114 def setUp(self):
115 arcsec = 1/3600
116 self.srcs1 = _make_frame([
117 {"diaSourceId": 1, "ra": 10.0, "dec": -5.0 + 0.3*arcsec},
118 {"diaSourceId": 2, "ra": 10.0, "dec": -5.0 - 0.1*arcsec},
119 ], visit=1, detector=0)
120 self.srcs2 = _make_frame([
121 {"diaSourceId": 11, "ra": 10.0, "dec": -5.0},
122 {"diaSourceId": 12, "ra": 10.0, "dec": -5.0 + 0.7*arcsec},
123 ], visit=1, detector=0)
125 def test_nearest_neighbor_shares_partner(self):
126 matched, _, unique2 = match_catalogs(self.srcs1, self.srcs2, radius=1*u.arcsec)
127 self.assertEqual(list(matched["diaSourceId_2"]), [11, 11])
128 self.assertEqual(list(unique2["diaSourceId"]), [12])
130 def test_one_to_one(self):
131 matched, unique1, unique2 = match_catalogs(self.srcs1, self.srcs2, radius=1*u.arcsec,
132 one_to_one=True)
133 # The closest pair (2, 11) is taken first, leaving 1 with 12.
134 self.assertEqual(dict(zip(matched["diaSourceId"], matched["diaSourceId_2"])), {1: 12, 2: 11})
135 self.assertEqual(matched["diaSourceId_2"].dtype, self.srcs2["diaSourceId"].dtype)
136 np.testing.assert_allclose(matched.set_index("diaSourceId")["xmatch_dist_arcsec"],
137 [0.4, 0.1], atol=1e-6)
138 self.assertEqual(len(unique1), 0)
139 self.assertEqual(len(unique2), 0)
141 def test_one_to_one_radius(self):
142 matched, unique1, unique2 = match_catalogs(self.srcs1, self.srcs2, radius=0.2*u.arcsec,
143 one_to_one=True)
144 self.assertEqual(list(matched["diaSourceId"]), [2])
145 self.assertEqual(list(unique1["diaSourceId"]), [1])
146 self.assertEqual(list(unique2["diaSourceId"]), [12])
147 self.assertNotIn("xmatch_dist_arcsec", unique1.columns)
150class TestFluxResiduals(lsst.utils.tests.TestCase):
151 """`flux_residuals` joins matched pairs with catalog 2 to compute
152 flux differences in units of sigma."""
154 def test_zero_residuals_when_identical(self):
155 srcs1 = _make_frame([
156 {"diaSourceId": 1, "ra": 10.0, "dec": 0.0, "visit": 1, "detector": 0,
157 "psfFlux": 100.0, "psfFluxErr": 5.0},
158 {"diaSourceId": 2, "ra": 10.001, "dec": 0.0, "visit": 1, "detector": 0,
159 "psfFlux": 200.0, "psfFluxErr": 10.0},
160 ])
161 srcs2 = _make_frame([
162 {"diaSourceId": 11, "ra": 10.0, "dec": 0.0, "visit": 1, "detector": 0,
163 "psfFlux": 100.0, "psfFluxErr": 5.0},
164 {"diaSourceId": 12, "ra": 10.001, "dec": 0.0, "visit": 1, "detector": 0,
165 "psfFlux": 200.0, "psfFluxErr": 10.0},
166 ])
167 matched, _, _ = match_catalogs(srcs1, srcs2, radius=1*u.arcsec)
168 residuals = flux_residuals(matched, srcs2)
169 np.testing.assert_array_equal(residuals["delta_flux"], [0.0, 0.0])
170 np.testing.assert_array_equal(residuals["delta_flux_sigma"], [0.0, 0.0])
172 def test_known_offset(self):
173 # Set up a 3-sigma flux offset on a single matched pair.
174 srcs1 = _make_frame([
175 {"diaSourceId": 1, "ra": 10.0, "dec": 0.0, "visit": 1, "detector": 0,
176 "psfFlux": 130.0, "psfFluxErr": 5.0},
177 ])
178 srcs2 = _make_frame([
179 {"diaSourceId": 11, "ra": 10.0, "dec": 0.0, "visit": 1, "detector": 0,
180 "psfFlux": 100.0, "psfFluxErr": 5.0},
181 ])
182 matched, _, _ = match_catalogs(srcs1, srcs2, radius=1*u.arcsec)
183 residuals = flux_residuals(matched, srcs2)
184 self.assertEqual(float(residuals["delta_flux"].iloc[0]), 30.0)
185 # Combined sigma is sqrt(5^2 + 5^2); 30/sqrt(50) ~ 4.243.
186 self.assertAlmostEqual(
187 float(residuals["delta_flux_sigma"].iloc[0]),
188 30.0 / np.sqrt(50.0),
189 places=10,
190 )
193class TestMatchToTruth(lsst.utils.tests.TestCase):
194 """`match_to_truth` reports purity and completeness against a truth
195 catalog by matching in both directions."""
197 def setUp(self):
198 # Three detected sources: two near truth, one in empty sky.
199 self.srcs = _make_frame([
200 {"diaSourceId": 1, "ra": 10.0, "dec": 0.0},
201 {"diaSourceId": 2, "ra": 11.0, "dec": 0.0},
202 {"diaSourceId": 3, "ra": 99.0, "dec": -50.0}, # bogus, not real
203 ])
204 # Three truth sources: two recovered, one missed.
205 self.truth = _make_frame([
206 {"injection_id": 100, "ra": 10.0, "dec": 0.0},
207 {"injection_id": 101, "ra": 11.0, "dec": 0.0},
208 {"injection_id": 102, "ra": 50.0, "dec": -10.0}, # not recovered
209 ])
211 def test_purity_and_completeness(self):
212 result = match_to_truth(self.srcs, self.truth, radius=1*u.arcsec)
213 # 2 of 3 detections are real; 2 of 3 truths are detected.
214 self.assertAlmostEqual(result["purity"], 2/3, places=10)
215 self.assertAlmostEqual(result["completeness"], 2/3, places=10)
216 # The bogus detection has is_real=False and a NA partner.
217 bogus = result["srcs"].set_index("diaSourceId").loc[3]
218 self.assertFalse(bool(bogus["is_real"]))
219 self.assertTrue(pd.isna(bogus["injection_id_match"]))
220 # The unrecovered truth has detected=False.
221 miss = result["truth"].set_index("injection_id").loc[102]
222 self.assertFalse(bool(miss["detected"]))
225class TestMemory(lsst.utils.tests.MemoryTestCase):
226 pass
229def setup_module(module):
230 lsst.utils.tests.init()
233if __name__ == "__main__": 233 ↛ 234line 233 didn't jump to line 234 because the condition on line 233 was never true
234 lsst.utils.tests.init()
235 unittest.main()