Coverage for tests/test_compare.py: 96%
80 statements
« prev ^ index » next coverage.py v7.16.1, created at 2026-09-22 10:41 +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 TestFluxResiduals(lsst.utils.tests.TestCase):
111 """`flux_residuals` joins matched pairs with catalog 2 to compute
112 flux differences in units of sigma."""
114 def test_zero_residuals_when_identical(self):
115 srcs1 = _make_frame([
116 {"diaSourceId": 1, "ra": 10.0, "dec": 0.0, "visit": 1, "detector": 0,
117 "psfFlux": 100.0, "psfFluxErr": 5.0},
118 {"diaSourceId": 2, "ra": 10.001, "dec": 0.0, "visit": 1, "detector": 0,
119 "psfFlux": 200.0, "psfFluxErr": 10.0},
120 ])
121 srcs2 = _make_frame([
122 {"diaSourceId": 11, "ra": 10.0, "dec": 0.0, "visit": 1, "detector": 0,
123 "psfFlux": 100.0, "psfFluxErr": 5.0},
124 {"diaSourceId": 12, "ra": 10.001, "dec": 0.0, "visit": 1, "detector": 0,
125 "psfFlux": 200.0, "psfFluxErr": 10.0},
126 ])
127 matched, _, _ = match_catalogs(srcs1, srcs2, radius=1*u.arcsec)
128 residuals = flux_residuals(matched, srcs2)
129 np.testing.assert_array_equal(residuals["delta_flux"], [0.0, 0.0])
130 np.testing.assert_array_equal(residuals["delta_flux_sigma"], [0.0, 0.0])
132 def test_known_offset(self):
133 # Set up a 3-sigma flux offset on a single matched pair.
134 srcs1 = _make_frame([
135 {"diaSourceId": 1, "ra": 10.0, "dec": 0.0, "visit": 1, "detector": 0,
136 "psfFlux": 130.0, "psfFluxErr": 5.0},
137 ])
138 srcs2 = _make_frame([
139 {"diaSourceId": 11, "ra": 10.0, "dec": 0.0, "visit": 1, "detector": 0,
140 "psfFlux": 100.0, "psfFluxErr": 5.0},
141 ])
142 matched, _, _ = match_catalogs(srcs1, srcs2, radius=1*u.arcsec)
143 residuals = flux_residuals(matched, srcs2)
144 self.assertEqual(float(residuals["delta_flux"].iloc[0]), 30.0)
145 # Combined sigma is sqrt(5^2 + 5^2); 30/sqrt(50) ~ 4.243.
146 self.assertAlmostEqual(
147 float(residuals["delta_flux_sigma"].iloc[0]),
148 30.0 / np.sqrt(50.0),
149 places=10,
150 )
153class TestMatchToTruth(lsst.utils.tests.TestCase):
154 """`match_to_truth` reports purity and completeness against a truth
155 catalog by matching in both directions."""
157 def setUp(self):
158 # Three detected sources: two near truth, one in empty sky.
159 self.srcs = _make_frame([
160 {"diaSourceId": 1, "ra": 10.0, "dec": 0.0},
161 {"diaSourceId": 2, "ra": 11.0, "dec": 0.0},
162 {"diaSourceId": 3, "ra": 99.0, "dec": -50.0}, # bogus, not real
163 ])
164 # Three truth sources: two recovered, one missed.
165 self.truth = _make_frame([
166 {"injection_id": 100, "ra": 10.0, "dec": 0.0},
167 {"injection_id": 101, "ra": 11.0, "dec": 0.0},
168 {"injection_id": 102, "ra": 50.0, "dec": -10.0}, # not recovered
169 ])
171 def test_purity_and_completeness(self):
172 result = match_to_truth(self.srcs, self.truth, radius=1*u.arcsec)
173 # 2 of 3 detections are real; 2 of 3 truths are detected.
174 self.assertAlmostEqual(result["purity"], 2/3, places=10)
175 self.assertAlmostEqual(result["completeness"], 2/3, places=10)
176 # The bogus detection has is_real=False and a NA partner.
177 bogus = result["srcs"].set_index("diaSourceId").loc[3]
178 self.assertFalse(bool(bogus["is_real"]))
179 self.assertTrue(pd.isna(bogus["injection_id_match"]))
180 # The unrecovered truth has detected=False.
181 miss = result["truth"].set_index("injection_id").loc[102]
182 self.assertFalse(bool(miss["detected"]))
185class TestMemory(lsst.utils.tests.MemoryTestCase):
186 pass
189def setup_module(module):
190 lsst.utils.tests.init()
193if __name__ == "__main__": 193 ↛ 194line 193 didn't jump to line 194 because the condition on line 193 was never true
194 lsst.utils.tests.init()
195 unittest.main()