Coverage for python/lsst/analysis/ap/compareDiaSources.py: 92%
249 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-22 15:23 +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/>.
22"""Compare the diaSources of two pipeline runs without an APDB.
24The sources come from the visit-level table that
25``consolidateDiaSourceTable`` writes (``dia_source_visit`` in AP), or from the
26detector-level tables it consolidates when a run has no visit table. Sources
27are matched by position, since diaSource IDs differ between runs. A source
28found in only one run is explained by looking for it in the other run's
29earlier catalogs, which record the detections that were removed before the
30table was written.
31"""
33__all__ = [
34 "DEFAULT_UPSTREAM_DATASETS",
35 "DiaSourceComparison",
36 "compare_dia_sources",
37 "plot_dia_source_comparison",
38]
40import dataclasses
42import astropy.units as u
43import numpy as np
44import pandas as pd
45from astropy.coordinates import SkyCoord, match_coordinates_sky
47from lsst.daf.butler import MissingDatasetTypeError
49from .compare import match_catalogs
51DEFAULT_UPSTREAM_DATASETS = {
52 "low_reliability": "dia_source_low_reliability",
53 "rejected": "rejected_dia_source",
54 "long_trail": "long_trailed_source_detector",
55}
56"""Catalogs of removed detections searched for an unmatched source, keyed
57by the reason recorded when the source is found in one (`dict` [`str`,
58`str`]). The first, the most common reason, is searched by position; the
59rest by the ID of the source's detection.
60"""
62_POSITION_SEARCHED = 1
63"""Number of catalogs at the start of the upstream datasets searched by
64position before the detection catalog is read (`int`). The rest are searched
65by the ID of the source's detection.
66"""
68_COMPARED_COLUMNS = ("psfFlux", "psfFluxErr", "apFlux", "scienceFlux", "snr", "reliability")
69"""Columns reported for both sources of a matched pair."""
71_REASONS = ("matched_elsewhere", *DEFAULT_UPSTREAM_DATASETS, "removed_later",
72 "not_detected", "no_data", "unknown")
73"""Every reason the default datasets can give, in the order they are
74tested.
75"""
77_GROUP_COLUMNS = ["instrument", "visit", "detector"]
80@dataclasses.dataclass
81class DiaSourceComparison:
82 """The diaSources of two runs, matched by position.
84 Attributes
85 ----------
86 summary : `pandas.DataFrame`
87 One row per instrument and visit, with the band of the visit
88 (``band``), the number of sources in each run (``n_ref``,
89 ``n_new``), matched (``n_matched``) and in one run only
90 (``n_ref_only``, ``n_new_only``), and the dataset each run's sources
91 were read from (``table_ref``, ``table_new``; `None` if the run has
92 no sources for the visit).
93 matched : `pandas.DataFrame`
94 One row per matched pair: the IDs and detectors of both sources,
95 their separation (``separation_arcsec``), the pixel offset of the new
96 source (``dx``, ``dy``; NaN if the detectors differ), each compared
97 column with ``_ref`` and ``_new`` suffixes, the PSF flux difference in
98 units of its uncertainty (``delta_psfFlux_sigma``), and the boolean
99 columns whose values differ (``flags_changed``, comma separated;
100 ``isNegative`` if the sign changed).
101 ref_only, new_only : `pandas.DataFrame`
102 Sources found in one run only, with all their columns plus
103 ``reason``, and ``other_diaSourceId``, ``other_separation_arcsec``
104 and ``other_reliability`` for the source the reason refers to (the
105 ID is that catalog's, not a diaSourceId, for a removed detection).
106 flag_changes : `pandas.DataFrame`
107 One row per boolean column that differs in at least one pair, with
108 the number of pairs where it is set in the reference run only
109 (``n_ref_only``) and in the new run only (``n_new_only``).
110 """
112 summary: pd.DataFrame
113 matched: pd.DataFrame
114 ref_only: pd.DataFrame
115 new_only: pd.DataFrame
116 flag_changes: pd.DataFrame
118 def reason_counts(self):
119 """Count the unmatched sources of each run by reason.
121 Returns
122 -------
123 counts : `pandas.DataFrame`
124 Indexed by reason, with columns ``ref_only`` and ``new_only``.
125 """
126 counts = pd.DataFrame({
127 "ref_only": self.ref_only["reason"].value_counts(),
128 "new_only": self.new_only["reason"].value_counts(),
129 }).fillna(0).astype(int)
130 order = [reason for reason in _REASONS if reason in counts.index]
131 return counts.loc[order + sorted(set(counts.index) - set(order))]
134def _query_refs(butler, dataset, visits=None):
135 """Return the refs of a dataset type, optionally for some visits only.
137 Parameters
138 ----------
139 butler : `lsst.daf.butler.Butler`
140 Butler to query.
141 dataset : `str`
142 Dataset type, with a ``visit`` dimension.
143 visits : `dict` [`str`, `collections.abc.Iterable` [`int`]], optional
144 Visits to find, keyed by instrument; all visits if not given.
146 Returns
147 -------
148 refs : `list` [`lsst.daf.butler.DatasetRef`]
149 The refs; none if the repository does not define the dataset type.
150 """
151 try:
152 if visits is None:
153 return list(butler.query_datasets(dataset, limit=None, explain=False))
154 refs = []
155 for instrument, instrument_visits in visits.items():
156 instrument_visits = sorted({int(visit) for visit in instrument_visits})
157 if instrument_visits: 157 ↛ 155line 157 didn't jump to line 155 because the condition on line 157 was always true
158 refs.extend(butler.query_datasets(
159 dataset, data_id={"instrument": instrument}, where="visit IN (:visits)",
160 bind={"visits": instrument_visits}, limit=None, explain=False))
161 return refs
162 except MissingDatasetTypeError:
163 return []
166def _read_table(butler, ref):
167 """Read a diaSource table as a DataFrame with an ``instrument`` column."""
168 table = butler.get(ref)
169 frame = table.to_pandas() if hasattr(table, "to_pandas") else pd.DataFrame(table)
170 frame["instrument"] = ref.dataId["instrument"]
171 return frame
174def _load_run(butler, dataset, detector_dataset, visits):
175 """Read the diaSources of one run.
177 Each visit is read from its ``dataset`` table, or from its
178 ``detector_dataset`` tables if it has none.
180 Returns
181 -------
182 sources : `pandas.DataFrame`
183 The diaSources, with an ``instrument`` column.
184 tables : `dict` [`tuple`, `str`]
185 The dataset each ``(instrument, visit)`` was read from.
186 bands : `dict` [`tuple`, `str`]
187 The band of each ``(instrument, visit)``.
188 """
189 if visits is not None:
190 # A visit ID alone does not name a visit, so look for it under every
191 # instrument.
192 instruments = [record.name for record in butler.query_dimension_records("instrument",
193 explain=False)]
194 visits = {instrument: list(visits) for instrument in instruments}
195 frames = []
196 tables = {}
197 bands = {}
198 for ref in _query_refs(butler, dataset, visits):
199 frames.append(_read_table(butler, ref))
200 tables[(ref.dataId["instrument"], ref.dataId["visit"])] = dataset
201 bands[(ref.dataId["instrument"], ref.dataId["visit"])] = ref.dataId["band"]
202 if visits is not None:
203 visits = {instrument: [visit for visit in instrument_visits if (instrument, visit) not in tables]
204 for instrument, instrument_visits in visits.items()}
205 if visits is None or any(visits.values()):
206 for ref in _query_refs(butler, detector_dataset, visits):
207 key = (ref.dataId["instrument"], ref.dataId["visit"])
208 if tables.get(key, detector_dataset) == detector_dataset: 208 ↛ 206line 208 didn't jump to line 206 because the condition on line 208 was always true
209 frames.append(_read_table(butler, ref))
210 tables[key] = detector_dataset
211 bands[key] = ref.dataId["band"]
212 if not frames:
213 empty = pd.DataFrame(columns=["diaSourceId", "instrument", "visit", "detector", "ra", "dec"])
214 return empty, tables, bands
215 return pd.concat(frames, ignore_index=True), tables, bands
218def _read_catalog(butler, ref):
219 """Read the positions in a catalog of detections.
221 Parameters
222 ----------
223 butler : `lsst.daf.butler.Butler`
224 Butler to read from.
225 ref : `lsst.daf.butler.DatasetRef` or `None`
226 The catalog: an afw `~lsst.afw.table.SourceCatalog` or a table.
228 Returns
229 -------
230 catalog : `pandas.DataFrame` or `None`
231 Columns ``id``, ``ra`` and ``dec`` (degrees), and ``reliability``
232 (NaN where the catalog has none); `None` if ``ref`` is `None`.
233 """
234 if ref is None:
235 return None
236 catalog = butler.get(ref)
237 if hasattr(catalog, "schema"):
238 # Take the few columns needed from the afw catalog; converting the
239 # whole catalog to astropy costs more than reading it.
240 if not catalog.isContiguous(): 240 ↛ 241line 240 didn't jump to line 241 because the condition on line 240 was never true
241 catalog = catalog.copy(deep=True)
242 names = set(catalog.schema.getNames())
243 else:
244 names = set(catalog.colnames)
245 n = len(catalog)
246 if "coord_ra" in names: 246 ↛ 250line 246 didn't jump to line 250 because the condition on line 246 was always true
247 ra = np.degrees(np.asarray(catalog["coord_ra"], dtype=float))
248 dec = np.degrees(np.asarray(catalog["coord_dec"], dtype=float))
249 else:
250 ra = np.asarray(catalog["ra"], dtype=float)
251 dec = np.asarray(catalog["dec"], dtype=float)
252 return pd.DataFrame({
253 "id": np.asarray(catalog["id"] if "id" in names else catalog["diaSourceId"]),
254 "ra": ra,
255 "dec": dec,
256 "reliability": (np.asarray(catalog["reliability"], dtype=float) if "reliability" in names
257 else np.full(n, np.nan)),
258 })
261def _nearest(sources, catalog, radius):
262 """Find the nearest catalog entry within ``radius`` of each source.
264 Returns
265 -------
266 index : `numpy.ndarray`
267 Row of ``catalog`` for each source, or -1 if none is within
268 ``radius``.
269 separation : `numpy.ndarray`
270 Separation in arcseconds, NaN where there is no entry.
271 """
272 index = np.full(len(sources), -1)
273 separation = np.full(len(sources), np.nan)
274 if len(sources) == 0 or len(catalog) == 0: 274 ↛ 275line 274 didn't jump to line 275 because the condition on line 274 was never true
275 return index, separation
276 idx, sep, _ = match_coordinates_sky(SkyCoord(sources["ra"].to_numpy(), sources["dec"].to_numpy(),
277 unit="deg"),
278 SkyCoord(catalog["ra"].to_numpy(), catalog["dec"].to_numpy(),
279 unit="deg"))
280 sep = sep.to(u.arcsec).value
281 found = sep <= radius.to_value(u.arcsec)
282 index[found] = idx[found]
283 separation[found] = sep[found]
284 return index, separation
287def _explain_unmatched(unmatched, other, other_butler, radius, upstream_datasets,
288 detection_dataset):
289 """Find why each unmatched source has no partner in the other run.
291 Parameters
292 ----------
293 unmatched : `pandas.DataFrame`
294 Sources of one run with no partner.
295 other : `pandas.DataFrame`
296 All diaSources of the other run.
297 other_butler : `lsst.daf.butler.Butler`
298 Butler for the other run.
299 radius : `astropy.units.Quantity`
300 Search radius.
301 upstream_datasets : `dict` [`str`, `str`]
302 Catalogs of removed detections, keyed by reason; see
303 `compare_dia_sources`.
304 detection_dataset : `str`
305 Catalog of all detections.
307 Returns
308 -------
309 unmatched : `pandas.DataFrame`
310 Copy of ``unmatched`` with ``reason``, ``other_diaSourceId``,
311 ``other_separation_arcsec`` and ``other_reliability`` columns.
312 """
313 result = unmatched.copy()
314 result["reason"] = pd.Series(None, index=result.index, dtype=object)
315 result["other_diaSourceId"] = pd.Series(pd.NA, index=result.index, dtype="Int64")
316 result["other_separation_arcsec"] = np.nan
317 result["other_reliability"] = np.nan
319 def assign(rows, reason, catalog):
320 """Record ``reason`` for the rows that have a catalog entry
321 nearby, and return the rows that still have no reason.
322 """
323 index, separation = _nearest(rows, catalog, radius)
324 hit = index >= 0
325 labels = rows.index[hit]
326 result.loc[labels, "reason"] = reason
327 result.loc[labels, "other_diaSourceId"] = catalog["id"].to_numpy()[index[hit]]
328 result.loc[labels, "other_separation_arcsec"] = separation[hit]
329 result.loc[labels, "other_reliability"] = catalog["reliability"].to_numpy()[index[hit]]
330 return result.loc[rows.index][result.loc[rows.index, "reason"].isna()]
332 other_positions = pd.DataFrame({
333 "id": other["diaSourceId"].to_numpy(),
334 "ra": other["ra"].to_numpy(dtype=float),
335 "dec": other["dec"].to_numpy(dtype=float),
336 "reliability": (other["reliability"].to_numpy(dtype=float) if "reliability" in other
337 else np.full(len(other), np.nan)),
338 "instrument": other["instrument"].to_numpy(),
339 "visit": other["visit"].to_numpy(),
340 "detector": other["detector"].to_numpy(),
341 })
342 # One query per catalog for every detector with an unmatched source.
343 catalog_datasets = [*upstream_datasets.values(), detection_dataset]
344 visits = {instrument: rows["visit"].unique()
345 for instrument, rows in unmatched.groupby("instrument")}
346 catalog_refs = {}
347 for dataset in catalog_datasets:
348 for ref in _query_refs(other_butler, dataset, visits):
349 data_id = ref.dataId
350 catalog_refs[dataset, data_id["instrument"], data_id["visit"], data_id["detector"]] = ref
352 def read(dataset, key):
353 """Read a catalog of the other run for an (instrument, visit,
354 detector), or return `None` if it was not written.
355 """
356 return _read_catalog(other_butler, catalog_refs.get((dataset, *key)))
358 upstream = list(upstream_datasets.items())
359 for (instrument, visit, detector), rows in unmatched.groupby(_GROUP_COLUMNS):
360 same_visit = other_positions[(other_positions["instrument"] == instrument)
361 & (other_positions["visit"] == visit)]
362 remaining = assign(rows, "matched_elsewhere", same_visit)
363 if not len(remaining):
364 continue
365 key = (instrument, visit, detector)
366 written = False
367 for reason, dataset in upstream[:_POSITION_SEARCHED]:
368 if (catalog := read(dataset, key)) is not None:
369 written = True
370 remaining = assign(remaining, reason, catalog)
371 if not len(remaining): 371 ↛ 372line 371 didn't jump to line 372 because the condition on line 371 was never true
372 continue
373 detections = read(detection_dataset, key)
374 if detections is None:
375 # Without the detection catalog, search the rest by position.
376 for reason, dataset in upstream[_POSITION_SEARCHED:]:
377 if len(remaining) and (catalog := read(dataset, key)) is not None: 377 ↛ 378line 377 didn't jump to line 378 because the condition on line 377 was never true
378 written = True
379 remaining = assign(remaining, reason, catalog)
380 if len(remaining): 380 ↛ 383line 380 didn't jump to line 383 because the condition on line 380 was always true
381 no_data = not written and not (same_visit["detector"] == detector).any()
382 result.loc[remaining.index, "reason"] = "no_data" if no_data else "unknown"
383 continue
384 # A source absent from the detection catalog was not detected. One in
385 # it was removed by whichever catalog lists its detection ID, since
386 # removed detections keep their IDs.
387 undetected = assign(remaining, "removed_later", detections)
388 result.loc[undetected.index, "reason"] = "not_detected"
389 unresolved = remaining.index.difference(undetected.index)
390 for reason, dataset in upstream[_POSITION_SEARCHED:]:
391 if unresolved.empty: 391 ↛ 392line 391 didn't jump to line 392 because the condition on line 391 was never true
392 break
393 if (catalog := read(dataset, key)) is None: 393 ↛ 394line 393 didn't jump to line 394 because the condition on line 393 was never true
394 continue
395 reliability = catalog.set_index("id")["reliability"]
396 ids = result.loc[unresolved, "other_diaSourceId"]
397 listed = ids.isin(reliability.index).to_numpy()
398 labels = unresolved[listed]
399 result.loc[labels, "reason"] = reason
400 result.loc[labels, "other_reliability"] = reliability.loc[ids[listed].to_numpy()].to_numpy()
401 unresolved = unresolved[~listed]
402 return result
405def _flag_columns(ref, new):
406 """Return the boolean columns both tables have."""
407 return [column for column in ref.columns
408 if column in new.columns
409 and pd.api.types.is_bool_dtype(ref[column]) and pd.api.types.is_bool_dtype(new[column])]
412def _pair_table(pairs, ref, new):
413 """Build the table of matched pairs and the flag-change counts."""
414 partners = new.set_index("diaSourceId").loc[pairs["diaSourceId_2"].to_numpy()]
415 matched = pd.DataFrame({
416 "instrument": pairs["instrument"].to_numpy(),
417 "visit": pairs["visit"].to_numpy(),
418 "diaSourceId_ref": pairs["diaSourceId"].to_numpy(),
419 "diaSourceId_new": partners.index.to_numpy(),
420 "detector_ref": pairs["detector"].to_numpy(),
421 "detector_new": partners["detector"].to_numpy(),
422 "separation_arcsec": pairs["xmatch_dist_arcsec"].to_numpy(dtype=float),
423 })
424 same_detector = matched["detector_ref"] == matched["detector_new"]
425 for axis in ("x", "y"):
426 if axis in pairs and axis in partners: 426 ↛ 425line 426 didn't jump to line 425 because the condition on line 426 was always true
427 offset = partners[axis].to_numpy(dtype=float) - pairs[axis].to_numpy(dtype=float)
428 matched[f"d{axis}"] = np.where(same_detector, offset, np.nan)
429 for column in _COMPARED_COLUMNS:
430 if column in pairs and column in partners:
431 matched[f"{column}_ref"] = pairs[column].to_numpy(dtype=float)
432 matched[f"{column}_new"] = partners[column].to_numpy(dtype=float)
433 if {"psfFlux_ref", "psfFluxErr_ref", "psfFlux_new", "psfFluxErr_new"} <= set(matched.columns): 433 ↛ 439line 433 didn't jump to line 439 because the condition on line 433 was always true
434 sigma = np.hypot(matched["psfFluxErr_ref"], matched["psfFluxErr_new"])
435 with np.errstate(divide="ignore", invalid="ignore"):
436 matched["delta_psfFlux_sigma"] = np.where(
437 sigma > 0, (matched["psfFlux_new"] - matched["psfFlux_ref"])/sigma, np.nan)
439 flags = _flag_columns(ref, new)
440 set_ref = pairs[flags].to_numpy(dtype=bool)
441 set_new = partners[flags].to_numpy(dtype=bool)
442 changed = set_ref != set_new
443 matched["flags_changed"] = [",".join(np.asarray(flags)[row]) for row in changed]
444 flag_changes = pd.DataFrame({
445 "flag": flags,
446 "n_ref_only": (set_ref & ~set_new).sum(axis=0),
447 "n_new_only": (~set_ref & set_new).sum(axis=0),
448 })
449 flag_changes = flag_changes[(flag_changes["n_ref_only"] + flag_changes["n_new_only"]) > 0]
450 flag_changes = flag_changes.sort_values(["n_ref_only", "n_new_only", "flag"],
451 ascending=[False, False, True], ignore_index=True)
452 return matched, flag_changes
455def _summarize(ref, new, matched, ref_only, new_only, tables_ref, tables_new, bands):
456 """Count the sources of each visit."""
457 keys = sorted(set(tables_ref) | set(tables_new))
459 def count(frame, columns=("instrument", "visit")):
460 return frame.groupby(list(columns)).size() if len(frame) else pd.Series(dtype=int)
462 counts = {
463 "n_ref": count(ref),
464 "n_new": count(new),
465 "n_matched": count(matched),
466 "n_ref_only": count(ref_only),
467 "n_new_only": count(new_only),
468 }
469 rows = []
470 for key in keys:
471 row = {"instrument": key[0], "visit": key[1], "band": bands.get(key)}
472 row.update({name: int(series.get(key, 0)) for name, series in counts.items()})
473 row.update(table_ref=tables_ref.get(key), table_new=tables_new.get(key))
474 rows.append(row)
475 columns = ["instrument", "visit", "band", *counts, "table_ref", "table_new"]
476 return pd.DataFrame(rows, columns=columns)
479def compare_dia_sources(butler_ref, butler_new, *, dataset="dia_source_visit",
480 detector_dataset="dia_source_detector", visits=None, radius=0.5*u.arcsec,
481 explain_unmatched=True,
482 upstream_datasets=DEFAULT_UPSTREAM_DATASETS,
483 detection_dataset="dia_source_unfiltered"):
484 """Match the diaSources of two runs and explain the differences.
486 Parameters
487 ----------
488 butler_ref, butler_new : `lsst.daf.butler.Butler`
489 Butlers whose default collections hold the reference and new runs.
490 dataset : `str`, optional
491 Visit-level diaSource table; ``diaSourceTable`` in DRP.
492 detector_dataset : `str`, optional
493 Detector-level diaSource table read for a visit that has no
494 ``dataset`` table, as in a run that skipped
495 ``consolidateDiaSourceTable``.
496 visits : `collections.abc.Iterable` [`int`], optional
497 Visits to compare; all visits if not given.
498 radius : `astropy.units.Quantity`, optional
499 Largest separation of a matched pair.
500 explain_unmatched : `bool`, optional
501 Look for each unmatched source in the other run's earlier catalogs
502 to find why it has no partner. This reads those catalogs for every
503 detector with an unmatched source. If False, every ``reason`` is
504 `None`.
505 upstream_datasets : `dict` [`str`, `str`], optional
506 Catalogs of removed detections, keyed by the reason recorded when a
507 source is found in one. The first is searched by position, so list
508 the most common reason first; the rest are searched by the ID of the
509 source's detection in ``detection_dataset``.
510 detection_dataset : `str`, optional
511 Catalog of all detections, whose IDs the removed detections keep. A
512 source found only there has reason ``removed_later``; one absent from
513 it has ``not_detected``.
515 Returns
516 -------
517 comparison : `DiaSourceComparison`
518 The matched pairs, unmatched sources and summaries.
520 Raises
521 ------
522 LookupError
523 Raised if neither run has any diaSources.
525 Notes
526 -----
527 Sources are paired one to one within each visit, closest pairs first
528 (`match_catalogs` with ``one_to_one=True``), so a source that lands on a
529 neighboring detector in the other run still matches.
531 The reason for an unmatched source is the first that applies:
533 ``matched_elsewhere``
534 The other run has a source within ``radius``, which was paired with
535 a closer source (a split or merged detection).
536 a key of ``upstream_datasets``
537 The other run detected the source and removed it; by default
538 ``low_reliability`` (``other_reliability`` is its score),
539 ``rejected`` (by ``filterDiaSource``) or ``long_trail``.
540 ``removed_later``
541 The other run detected the source but it is in none of those
542 catalogs.
543 ``not_detected``
544 The other run's detection catalog has no source there.
545 ``no_data``
546 The other run has none of these catalogs for the detector, nor any
547 sources on it.
548 ``unknown``
549 The other run lacks the detection catalog but wrote others.
550 """
551 ref, tables_ref, bands_ref = _load_run(butler_ref, dataset, detector_dataset, visits)
552 new, tables_new, bands_new = _load_run(butler_new, dataset, detector_dataset, visits)
553 if len(ref) == 0 and len(new) == 0:
554 raise LookupError(f"Neither run has any {dataset!r} or {detector_dataset!r} diaSources.")
555 pairs, ref_only, new_only = match_catalogs(ref, new, radius=radius, on=("instrument", "visit"),
556 id_col="diaSourceId", one_to_one=True)
557 matched, flag_changes = _pair_table(pairs, ref, new)
558 if explain_unmatched:
559 ref_only = _explain_unmatched(ref_only, new, butler_new, radius, upstream_datasets,
560 detection_dataset)
561 new_only = _explain_unmatched(new_only, ref, butler_ref, radius, upstream_datasets,
562 detection_dataset)
563 else:
564 ref_only = ref_only.assign(reason=None)
565 new_only = new_only.assign(reason=None)
566 # A visit's band is the same in both runs; the reference run's is used
567 # where both have the visit.
568 summary = _summarize(ref, new, matched, ref_only, new_only, tables_ref, tables_new,
569 bands_new | bands_ref)
570 return DiaSourceComparison(summary=summary, matched=matched.reset_index(drop=True),
571 ref_only=ref_only.reset_index(drop=True),
572 new_only=new_only.reset_index(drop=True), flag_changes=flag_changes)
575def plot_dia_source_comparison(comparison, labels=("ref", "new")):
576 """Plot the differences between the diaSources of two runs.
578 Parameters
579 ----------
580 comparison : `DiaSourceComparison`
581 Result of `compare_dia_sources`.
582 labels : `tuple` [`str`, `str`], optional
583 Names of the reference and new runs.
585 Returns
586 -------
587 figure : `matplotlib.figure.Figure`
588 Four panels: PSF flux difference against signal to noise for
589 matched pairs, reliability in each run for matched pairs, the
590 number of unmatched sources by reason, and the sky positions of the
591 unmatched sources.
592 """
593 from matplotlib.figure import Figure
595 label_ref, label_new = labels
596 matched = comparison.matched
597 figure = Figure(figsize=(11, 9), layout="constrained")
598 (flux_ax, reliability_ax), (reason_ax, sky_ax) = figure.subplots(2, 2)
600 if "delta_psfFlux_sigma" in matched and "snr_ref" in matched: 600 ↛ 605line 600 didn't jump to line 605 because the condition on line 600 was always true
601 flux_ax.scatter(matched["snr_ref"], matched["delta_psfFlux_sigma"], s=6)
602 flux_ax.set_xscale("log")
603 for level in (-1, 0, 1):
604 flux_ax.axhline(level, color="gray", lw=0.8, ls="-" if level == 0 else ":")
605 flux_ax.set_xlabel(f"snr ({label_ref})")
606 flux_ax.set_ylabel(f"psfFlux ({label_new} - {label_ref}) / sigma")
607 flux_ax.set_title(f"{len(matched)} matched pairs")
609 if "reliability_ref" in matched: 609 ↛ 612line 609 didn't jump to line 612 because the condition on line 609 was always true
610 reliability_ax.scatter(matched["reliability_ref"], matched["reliability_new"], s=6)
611 reliability_ax.plot([0, 1], [0, 1], color="gray", lw=0.8)
612 reliability_ax.set_xlabel(f"reliability ({label_ref})")
613 reliability_ax.set_ylabel(f"reliability ({label_new})")
615 counts = comparison.reason_counts()
616 positions = np.arange(len(counts))
617 reason_ax.barh(positions - 0.2, counts["ref_only"], height=0.4, label=f"{label_ref} only")
618 reason_ax.barh(positions + 0.2, counts["new_only"], height=0.4, label=f"{label_new} only")
619 reason_ax.set_yticks(positions, [str(reason) for reason in counts.index])
620 reason_ax.invert_yaxis()
621 reason_ax.set_xlabel("unmatched sources")
622 reason_ax.legend()
624 for frame, marker, label in ((comparison.ref_only, "o", label_ref),
625 (comparison.new_only, "x", label_new)):
626 sky_ax.scatter(frame["ra"], frame["dec"], s=12, marker=marker, label=f"{label} only")
627 sky_ax.invert_xaxis()
628 sky_ax.set_xlabel("ra (deg)")
629 sky_ax.set_ylabel("dec (deg)")
630 sky_ax.legend()
631 return figure