Coverage for python/lsst/meas/astrom/matcher_probabilistic.py: 79%
271 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-07-25 08:55 +0000
1# This file is part of meas_astrom.
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__all__ = [
23 "CatalogExtras", "ComparableCatalog", "ConvertCatalogCoordinatesConfig", "MatchProbabilisticConfig",
24 "MatcherProbabilistic",
25]
27from dataclasses import dataclass
28import logging
29import time
30from typing import Callable
31import warnings
33import astropy.table
34import lsst.pex.config as pexConfig
35import numpy as np
36from scipy.spatial import cKDTree
37from smatch.matcher import Matcher
39logger_default = logging.getLogger(__name__)
42def _mul_column(column: np.array, value: float):
43 if value is not None and value != 1: 43 ↛ 44line 43 didn't jump to line 44 because the condition on line 43 was never true
44 column *= value
45 return column
48def _radec_to_xyz(ra, dec):
49 """Convert input ra/dec coordinates to spherical unit vectors.
51 Parameters
52 ----------
53 ra, dec: `numpy.ndarray`
54 Arrays of right ascension/declination in degrees.
56 Returns
57 -------
58 vectors : `numpy.ndarray`, (N, 3)
59 Output unit vectors.
60 """
61 if ra.size != dec.size:
62 raise ValueError('ra and dec must be same size')
63 ras = np.radians(ra)
64 decs = np.radians(dec)
65 vectors = np.empty((ras.size, 3))
67 sin_dec = np.sin(np.pi / 2 - decs)
68 vectors[:, 0] = sin_dec * np.cos(ras)
69 vectors[:, 1] = sin_dec * np.sin(ras)
70 vectors[:, 2] = np.cos(np.pi / 2 - decs)
72 return vectors
75@dataclass
76class CatalogExtras:
77 """Store frequently-reference (meta)data relevant for matching a catalog.
79 Parameters
80 ----------
81 catalog : `astropy.table.Table`
82 A Table to store extra information for.
83 select : `numpy.array`
84 A numpy boolean array of the same length as catalog to be used for
85 target selection.
86 """
88 n: int
89 indices: np.array
90 select: np.array
92 coordinate_factor: float = None
94 def __init__(
95 self,
96 catalog: astropy.table.Table,
97 select: np.array = None,
98 coordinate_factor: float = None,
99 ):
100 self.n = len(catalog)
101 self.select = np.ones(self.n, dtype=bool) if select is None else select
102 self.indices = np.flatnonzero(select) if select is not None else np.arange(self.n)
103 self.coordinate_factor = coordinate_factor
106@dataclass(frozen=True)
107class ComparableCatalog:
108 """A catalog with sources with coordinate columns in some standard format/units.
110 catalog : `astropy.table.Table`
111 A catalog with comparable coordinate columns.
112 column_coord1 : `str`
113 The first spatial coordinate column name.
114 column_coord2 : `str`
115 The second spatial coordinate column name.
116 coord1 : `numpy.array`
117 The first spatial coordinate values.
118 coord2 : `numpy.array`
119 The second spatial coordinate values.
120 extras : `CatalogExtras`
121 Extra cached (meta)data for the `catalog`.
122 """
124 catalog: astropy.table.Table
125 column_coord1: str
126 column_coord2: str
127 coord1: np.array
128 coord2: np.array
129 extras: CatalogExtras
132class ConvertCatalogCoordinatesConfig(pexConfig.Config):
133 """Configuration for the MatchProbabilistic matcher."""
135 column_ref_coord1 = pexConfig.Field[str](
136 default='ra',
137 doc='The reference table column for the first spatial coordinate (usually x or ra).',
138 )
139 column_ref_coord2 = pexConfig.Field[str](
140 default='dec',
141 doc='The reference table column for the second spatial coordinate (usually y or dec).'
142 'Units must match column_ref_coord1.',
143 )
144 column_target_coord1 = pexConfig.Field[str](
145 default='coord_ra',
146 doc='The target table column for the first spatial coordinate (usually x or ra).'
147 'Units must match column_ref_coord1.',
148 )
149 column_target_coord2 = pexConfig.Field[str](
150 default='coord_dec',
151 doc='The target table column for the second spatial coordinate (usually y or dec).'
152 'Units must match column_ref_coord2.',
153 )
154 coords_spherical = pexConfig.Field[bool](
155 default=True,
156 doc='Whether column_*_coord[12] are spherical coordinates (ra/dec) or not (pixel x/y).',
157 )
158 coords_ref_factor = pexConfig.Field[float](
159 default=1.0,
160 doc='Multiplicative factor for reference catalog coordinates.'
161 'If coords_spherical is true, this must be the number of degrees per unit increment of '
162 'column_ref_coord[12]. Otherwise, it must convert the coordinate to the same units'
163 ' as the target coordinates.',
164 )
165 coords_target_factor = pexConfig.Field[float](
166 default=1.0,
167 doc='Multiplicative factor for target catalog coordinates.'
168 'If coords_spherical is true, this must be the number of degrees per unit increment of '
169 'column_target_coord[12]. Otherwise, it must convert the coordinate to the same units'
170 ' as the reference coordinates.',
171 )
172 coords_ref_to_convert = pexConfig.DictField[str, str](
173 default=None,
174 optional=True,
175 dictCheck=lambda x: len(x) == 2,
176 doc='Dict mapping sky coordinate columns to be converted to pixel columns.',
177 )
178 mag_zeropoint_ref = pexConfig.Field[float](
179 default=31.4,
180 doc='Magnitude zeropoint for reference catalog.',
181 )
182 return_converted_coords = pexConfig.Field[float](
183 default=True,
184 doc='Whether to return converted coordinates for matching or only write them.',
185 )
187 def format_catalogs(
188 self,
189 catalog_ref: astropy.table.Table,
190 catalog_target: astropy.table.Table,
191 select_ref: np.array = None,
192 select_target: np.array = None,
193 radec_to_xy_func: Callable = None,
194 **kwargs,
195 ):
196 """Format matched catalogs that may require coordinate conversions.
198 Parameters
199 ----------
200 catalog_ref : `astropy.table.Table`
201 A reference catalog for comparison to `catalog_target`.
202 catalog_target : `astropy.table.Table`
203 A target catalog with measurements for comparison to `catalog_ref`.
204 select_ref : `numpy.ndarray`, (Nref,)
205 A boolean array of len `catalog_ref`, True for valid match candidates.
206 select_target : `numpy.ndarray`, (Ntarget,)
207 A boolean array of len `catalog_target`, True for valid match candidates.
208 radec_to_xy_func : `typing.Callable`
209 Function taking equal-length ra, dec arrays and returning an ndarray of
210 - ``x``: current parameter (`float`).
211 - ``extra_args``: additional arguments (`dict`).
212 kwargs
213 Additional keyword arguments to pass to radec_to_xy_func.
215 Returns
216 -------
217 compcat_ref, compcat_target : `ComparableCatalog`
218 Comparable catalogs corresponding to the input reference and target.
219 """
220 convert_ref = self.coords_ref_to_convert
221 if convert_ref and not callable(radec_to_xy_func): 221 ↛ 222line 221 didn't jump to line 222 because the condition on line 221 was never true
222 raise TypeError('radec_to_xy_func must be callable if converting ref coords')
224 # Set up objects with frequently-used attributes like selection bool array
225 extras_ref, extras_target = (
226 CatalogExtras(catalog, select=select, coordinate_factor=coord_factor)
227 for catalog, select, coord_factor in zip(
228 (catalog_ref, catalog_target),
229 (select_ref, select_target),
230 (self.coords_ref_factor, self.coords_target_factor),
231 )
232 )
234 compcats = []
236 # Retrieve coordinates and multiply them by scaling factors
237 for catalog, extras, (column1, column2), convert in (
238 (catalog_ref, extras_ref, (self.column_ref_coord1, self.column_ref_coord2), convert_ref),
239 (catalog_target, extras_target, (self.column_target_coord1, self.column_target_coord2), False),
240 ):
241 coord1, coord2 = (
242 _mul_column(catalog[column], extras.coordinate_factor)
243 for column in (column1, column2)
244 )
245 if convert: 245 ↛ 246line 245 didn't jump to line 246 because the condition on line 245 was never true
246 xy_ref = radec_to_xy_func(coord1, coord2, self.coords_ref_factor, **kwargs)
247 for idx_coord, column_out in enumerate(self.coords_ref_to_convert.values()):
248 coord = np.array([xy[idx_coord] for xy in xy_ref])
249 catalog[column_out] = coord
250 if convert_ref: 250 ↛ 251line 250 didn't jump to line 251 because the condition on line 250 was never true
251 column1, column2 = self.coords_ref_to_convert.values()
252 if self.return_converted_coords:
253 coord1, coord2 = catalog[column1], catalog[column2]
255 compcats.append(ComparableCatalog(
256 catalog=catalog, column_coord1=column1, column_coord2=column2,
257 coord1=coord1, coord2=coord2, extras=extras,
258 ))
260 return compcats[0], compcats[1]
263class MatchProbabilisticConfig(pexConfig.Config):
264 """Configuration for the MatchProbabilistic matcher."""
266 column_ref_order = pexConfig.Field(
267 dtype=str,
268 default=None,
269 optional=True,
270 doc='Name of column in reference catalog specifying order for matching'
271 ' Derived from columns_ref_flux if not set.',
272 )
274 @property
275 def columns_in_ref(self) -> set[str]:
276 return set(self.columns_ordered_in_ref)
278 @property
279 def columns_ordered_in_ref(self) -> dict[str, None]:
280 columns_all = [
281 self.coord_format.column_ref_coord1,
282 self.coord_format.column_ref_coord2,
283 ]
284 for columns in (
285 self.columns_ref_flux,
286 self.columns_ref_meas,
287 self.columns_ref_select_false,
288 self.columns_ref_select_true,
289 self.columns_ref_copy,
290 ):
291 columns_all.extend(columns)
292 if self.column_ref_order: 292 ↛ 293line 292 didn't jump to line 293 because the condition on line 292 was never true
293 columns_all.append(self.column_ref_order)
295 return {k: None for k in columns_all}
297 @property
298 def columns_in_target(self) -> set[str]:
299 return set(self.columns_ordered_in_target)
301 @property
302 def columns_ordered_in_target(self) -> dict[str, None]:
303 columns_all = [
304 self.coord_format.column_target_coord1,
305 self.coord_format.column_target_coord2,
306 ]
307 for columns in (
308 self.columns_target_meas,
309 self.columns_target_err,
310 self.columns_target_select_false,
311 self.columns_target_select_true,
312 self.columns_target_copy,
313 ):
314 columns_all.extend(columns)
315 return {k: None for k in columns_all}
317 columns_ref_copy = pexConfig.ListField(
318 dtype=str,
319 default=[],
320 listCheck=lambda x: len(set(x)) == len(x),
321 optional=True,
322 doc='Reference table columns to copy unchanged into both match tables',
323 )
324 columns_ref_flux = pexConfig.ListField(
325 dtype=str,
326 default=[],
327 optional=True,
328 doc="List of reference flux columns to nansum total magnitudes from if column_order is None",
329 )
330 columns_ref_meas = pexConfig.ListField(
331 dtype=str,
332 doc='The reference table columns to compute match likelihoods from '
333 '(usually centroids and fluxes/magnitudes)',
334 )
335 columns_ref_select_true = pexConfig.ListField(
336 dtype=str,
337 default=tuple(),
338 doc='Reference table columns to require to be True for selecting sources',
339 )
340 columns_ref_select_false = pexConfig.ListField(
341 dtype=str,
342 default=tuple(),
343 doc='Reference table columns to require to be False for selecting sources',
344 )
345 columns_target_copy = pexConfig.ListField(
346 dtype=str,
347 default=[],
348 listCheck=lambda x: len(set(x)) == len(x),
349 optional=True,
350 doc='Target table columns to copy unchanged into both match tables',
351 )
352 columns_target_meas = pexConfig.ListField(
353 dtype=str,
354 doc='Target table columns with measurements corresponding to columns_ref_meas',
355 )
356 columns_target_err = pexConfig.ListField(
357 dtype=str,
358 doc='Target table columns with standard errors (sigma) corresponding to columns_ref_meas',
359 )
360 columns_target_select_true = pexConfig.ListField(
361 dtype=str,
362 default=[],
363 doc='Target table columns to require to be True for selecting sources',
364 )
365 columns_target_select_false = pexConfig.ListField(
366 dtype=str,
367 default=[],
368 doc='Target table columns to require to be False for selecting sources',
369 )
370 coord_format = pexConfig.ConfigField(
371 dtype=ConvertCatalogCoordinatesConfig,
372 doc="Configuration for coordinate conversion",
373 )
374 mag_brightest_ref = pexConfig.Field(
375 dtype=float,
376 default=-np.inf,
377 doc='Bright magnitude cutoff for selecting reference sources to match.'
378 ' Ignored if column_ref_order is None.'
379 )
380 mag_faintest_ref = pexConfig.Field(
381 dtype=float,
382 default=np.inf,
383 doc='Faint magnitude cutoff for selecting reference sources to match.'
384 ' Ignored if column_ref_order is None.'
385 )
386 match_dist_max = pexConfig.Field(
387 dtype=float,
388 default=0.5,
389 doc='Maximum match distance. Units must be arcseconds if coords_spherical, '
390 'or else match those of column_*_coord[12] multiplied by coords_*_factor.',
391 )
392 match_n_max = pexConfig.Field(
393 dtype=int,
394 default=10,
395 optional=True,
396 doc='Maximum number of spatial matches to consider (in ascending distance order).',
397 check=lambda x: x >= 1,
398 )
399 match_n_finite_min = pexConfig.Field(
400 dtype=int,
401 default=2,
402 optional=True,
403 doc='Minimum number of columns with a finite value to measure match likelihood',
404 check=lambda x: x >= 2,
405 )
406 order_ascending = pexConfig.Field(
407 dtype=bool,
408 default=False,
409 optional=True,
410 doc='Whether to order reference match candidates in ascending order of column_ref_order '
411 '(should be False if the column is a flux and True if it is a magnitude.',
412 )
414 def validate(self):
415 super().validate()
416 n_ref_meas = len(self.columns_ref_meas)
417 n_target_meas = len(self.columns_target_meas)
418 n_target_err = len(self.columns_target_err)
419 match_n_finite_min = self.match_n_finite_min
420 errors = []
421 if n_target_meas != n_ref_meas:
422 errors.append(f"{len(self.columns_target_meas)=} != {len(self.columns_ref_meas)=}")
423 if n_target_err != n_ref_meas:
424 errors.append(f"{len(self.columns_target_err)=} != {len(self.columns_ref_meas)=}")
425 if not (n_ref_meas >= match_n_finite_min): 425 ↛ 426line 425 didn't jump to line 426 because the condition on line 425 was never true
426 errors.append(
427 f"{len(self.columns_ref_meas)=} !>= {self.match_n_finite_min=}, no matches possible"
428 )
429 if self.column_ref_order is None:
430 if len(self.columns_ref_flux) == 0: 430 ↛ 431line 430 didn't jump to line 431 because the condition on line 430 was never true
431 errors.append("Must specify one of columns_ref_flux or column_ref_order")
432 if errors:
433 raise ValueError("\n".join(errors))
436def default_value(dtype):
437 if dtype is str: 437 ↛ 438line 437 didn't jump to line 438 because the condition on line 437 was never true
438 return ''
439 elif np.issubdtype(dtype, np.signedinteger): 439 ↛ 441line 439 didn't jump to line 441 because the condition on line 439 was always true
440 return np.inf
441 elif np.issubdtype(dtype, np.unsignedinteger):
442 return -np.inf
443 return None
446class MatcherProbabilistic:
447 """A probabilistic, greedy catalog matcher.
449 Parameters
450 ----------
451 config: `MatchProbabilisticConfig`
452 A configuration instance.
453 """
455 config: MatchProbabilisticConfig
457 def __init__(
458 self,
459 config: MatchProbabilisticConfig,
460 ):
461 self.config = config
463 def match(
464 self,
465 catalog_ref: astropy.table.Table,
466 catalog_target: astropy.table.Table,
467 select_ref: np.array = None,
468 select_target: np.array = None,
469 logger: logging.Logger = None,
470 logging_n_rows: int = None,
471 **kwargs
472 ):
473 """Match catalogs.
475 Parameters
476 ----------
477 catalog_ref : `astropy.table.Table`
478 A reference catalog to match in order of a given column (i.e. greedily).
479 catalog_target : `astropy.table.Table`
480 A target catalog for matching sources from `catalog_ref`. Must contain measurements with errors.
481 select_ref : `numpy.array`
482 A boolean array of the same length as `catalog_ref` selecting the sources that can be matched.
483 select_target : `numpy.array`
484 A boolean array of the same length as `catalog_target` selecting the sources that can be matched.
485 logger : `logging.Logger`
486 A Logger for logging.
487 logging_n_rows : `int`
488 The number of sources to match before printing a log message.
489 kwargs
490 Additional keyword arguments to pass to `format_catalogs`.
492 Returns
493 -------
494 catalog_out_ref : `astropy.table.Table`
495 A catalog of identical length to `catalog_ref`, containing match information for rows selected by
496 `select_ref` (including the matching row index in `catalog_target`).
497 catalog_out_target : `astropy.table.Table`
498 A catalog of identical length to `catalog_target`, containing the indices of matching rows in
499 `catalog_ref`.
500 exceptions : `dict` [`int`, `Exception`]
501 A dictionary keyed by `catalog_target` row number of the first exception caught when matching.
502 """
503 if logger is None: 503 ↛ 504line 503 didn't jump to line 504 because the condition on line 503 was never true
504 logger = logger_default
506 t_init = time.process_time()
507 config = self.config
509 # Transform any coordinates, if required
510 # Note: The returned objects contain the original catalogs, as well as
511 # transformed coordinates, and the selection of sources for matching.
512 # These might be identical to the arrays passed as kwargs, but that
513 # depends on config settings.
514 # For the rest of this function, the selection arrays will be used,
515 # but the indices of the original, unfiltered catalog will also be
516 # output, so some further indexing steps are needed.
517 with warnings.catch_warnings():
518 # We already issued a deprecation warning; no need to repeat it.
519 warnings.filterwarnings(action="ignore", category=FutureWarning)
520 ref, target = config.coord_format.format_catalogs(
521 catalog_ref=catalog_ref, catalog_target=catalog_target,
522 select_ref=select_ref, select_target=select_target,
523 **kwargs
524 )
526 # If no order is specified, take nansum of all flux columns for a 'total flux'
527 # Note: it won't actually be a total flux if bands overlap significantly
528 # (or it might define a filter with >100% efficiency
529 column_order = (
530 catalog_ref[config.column_ref_order][ref.extras.select]
531 if config.column_ref_order is not None else
532 np.nansum([catalog_ref[col][ref.extras.select] for col in config.columns_ref_flux], axis=0)
533 )
534 order = np.argsort(column_order if config.order_ascending else -column_order)
536 n_ref_select = len(ref.extras.indices)
538 coords_spherical = config.coord_format.coords_spherical
539 coords_ref, coords_target = (
540 (cat.coord1[cat.extras.select], cat.coord2[cat.extras.select])
541 for cat in (ref, target)
542 )
544 # Generate K-d tree to compute distances
545 logger.info('Generating cKDTree with match_n_max=%d', config.match_n_max)
547 if coords_spherical: 547 ↛ 558line 547 didn't jump to line 558 because the condition on line 547 was always true
548 match_dist_max = config.match_dist_max/3600.
549 with Matcher(coords_target[0], coords_target[1]) as matcher:
550 idxs_target_select = matcher.query_knn(
551 coords_ref[0], coords_ref[1],
552 distance_upper_bound=match_dist_max,
553 k=config.match_n_max,
554 )
555 # Call scipy for non-spherical case
556 # The spherical case won't trigger, but the implementation is left for comparison, if needed
557 else:
558 match_dist_max = np.radians(config.match_dist_max/3600.)
559 # Convert ra/dec sky coordinates to spherical vectors for accurate distances
560 func_convert = _radec_to_xyz if coords_spherical else np.vstack
561 vec_ref, vec_target = (
562 func_convert(coords[0], coords[1])
563 for coords in (coords_ref, coords_target)
564 )
565 tree_obj = cKDTree(vec_target)
566 _, idxs_target_select = tree_obj.query(
567 vec_ref,
568 distance_upper_bound=match_dist_max,
569 k=config.match_n_max,
570 )
572 n_target_select = len(target.extras.indices)
573 n_matches = np.sum(idxs_target_select != n_target_select, axis=1)
574 n_matched_max = np.sum(n_matches == config.match_n_max)
575 if n_matched_max > 0: 575 ↛ 576line 575 didn't jump to line 576 because the condition on line 575 was never true
576 logger.warning(
577 '%d/%d (%.2f%%) selected true objects have n_matches=n_match_max(%d)',
578 n_matched_max, n_ref_select, 100.*n_matched_max/n_ref_select, config.match_n_max
579 )
581 # Pre-allocate outputs
582 target_row_match = np.full(target.extras.n, np.iinfo(np.int64).min, dtype=np.int64)
583 ref_candidate_match = np.zeros(ref.extras.n, dtype=bool)
584 ref_row_match = np.full(ref.extras.n, np.iinfo(np.int64).min, dtype=np.int64)
585 ref_match_count = np.zeros(ref.extras.n, dtype=np.int32)
586 ref_match_meas_finite = np.zeros(ref.extras.n, dtype=np.int32)
587 ref_chisq = np.full(ref.extras.n, np.nan, dtype=float)
589 # Need the original reference row indices for output
590 idx_orig_ref, idx_orig_target = (np.argwhere(cat.extras.select)[:, 0] for cat in (ref, target))
592 # Retrieve required columns, including any converted ones (default to original column name)
593 columns_convert = config.coord_format.coords_ref_to_convert
594 if columns_convert is None: 594 ↛ 596line 594 didn't jump to line 596 because the condition on line 594 was always true
595 columns_convert = {}
596 data_ref = np.array([
597 ref.catalog[columns_convert.get(column, column)][ref.extras.indices[order]]
598 for column in config.columns_ref_meas
599 ])
600 data_target = np.array([
601 target.catalog[col][target.extras.select] for col in config.columns_target_meas
602 ])
603 errors_target = np.array([
604 target.catalog[col][target.extras.select] for col in config.columns_target_err
605 ])
607 exceptions = {}
608 # The kdTree uses len(inputs) as a sentinel value for no match
609 matched_target = {n_target_select, }
610 index_ref = idx_orig_ref[order]
611 # Fill in the candidate column
612 ref_candidate_match[index_ref] = True
614 # Count this as the time when disambiguation begins
615 t_begin = time.process_time()
617 # Exclude unmatched sources
618 matched_ref = idxs_target_select[order, 0] != n_target_select
619 order = order[matched_ref]
620 idx_first = idxs_target_select[order, 0]
621 data_ref_matched = data_ref[:, matched_ref]
622 chi_0 = (data_target[:, idx_first] - data_ref_matched)/errors_target[:, idx_first]
623 chi_finite_0 = np.isfinite(chi_0)
624 n_finite_0 = np.sum(chi_finite_0, axis=0)
625 chi_0[~chi_finite_0] = 0
626 chisq_sum_0 = np.sum(chi_0*chi_0, axis=0)
627 n_meas = len(config.columns_ref_meas)
628 n_ambiguous = 0
630 logger.info('Disambiguating %d/%d matches/targets', len(order), len(ref.catalog))
631 for index_n, index_row_select in enumerate(order):
632 index_row = idx_orig_ref[index_row_select]
633 found = idxs_target_select[index_row_select, :]
634 # Unambiguous match, short-circuit some evaluations
635 if (found[1] == n_target_select) and (found[0] not in matched_target):
636 n_finite = n_finite_0[index_n]
637 if not (n_finite >= config.match_n_finite_min): 637 ↛ 638line 637 didn't jump to line 638 because the condition on line 637 was never true
638 continue
639 idx_chisq_min = 0
640 n_matched = 1
641 chisq_sum = chisq_sum_0[index_n]
642 else:
643 # Select match candidates from nearby sources not already matched
644 # Note: set lookup is apparently fast enough that this is a few percent faster than:
645 # found = [x for x in found[found != n_target_select] if x not in matched_target]
646 # ... at least for ~1M sources
647 found = [x for x in found if x not in matched_target]
648 n_found = len(found)
649 if n_found == 0: 649 ↛ 650line 649 didn't jump to line 650 because the condition on line 649 was never true
650 continue
651 # This is an ndarray of n_found rows x len(data_ref/target) columns
652 chi = (
653 data_target[:, found] - data_ref_matched[:, index_n].reshape((n_meas, 1))
654 )/errors_target[:, found]
655 finite = np.isfinite(chi)
656 n_finite = np.sum(finite, axis=0)
657 # Require some number of finite chi_sq to match
658 chisq_good = n_finite >= config.match_n_finite_min
659 if not any(chisq_good): 659 ↛ 660line 659 didn't jump to line 660 because the condition on line 659 was never true
660 continue
661 try:
662 chisq_sum = np.zeros(n_found, dtype=float)
663 chisq_sum[chisq_good] = np.nansum(chi[:, chisq_good] ** 2, axis=0)
664 idx_chisq_min = np.nanargmin(chisq_sum / n_finite)
665 n_finite = n_finite[idx_chisq_min]
666 n_matched = len(chisq_good)
667 chisq_sum = chisq_sum[idx_chisq_min]
668 n_ambiguous += 1
669 except Exception as error:
670 # Can't foresee any exceptions, but they shouldn't prevent
671 # matching subsequent sources
672 exceptions[index_row] = error
673 ref_match_meas_finite[index_row] = n_finite
674 ref_match_count[index_row] = n_matched
675 ref_chisq[index_row] = chisq_sum
676 idx_match_select = found[idx_chisq_min]
677 row_target = target.extras.indices[idx_match_select]
678 ref_row_match[index_row] = row_target
680 target_row_match[row_target] = index_row
681 matched_target.add(idx_match_select)
683 if logging_n_rows and ((index_n + 1) % logging_n_rows == 0):
684 t_elapsed = time.process_time() - t_begin
685 logger.info(
686 'Processed %d/%d in %.2fs at sort value=%.3f',
687 index_n + 1, n_ref_select, t_elapsed, column_order[order[index_n]],
688 )
690 data_ref = {
691 'match_candidate': ref_candidate_match,
692 'match_row': ref_row_match,
693 'match_count': ref_match_count,
694 'match_chi2': ref_chisq,
695 'match_n_chi2_finite': ref_match_meas_finite,
696 }
697 data_target = {
698 'match_candidate': target.extras.select if target.extras.select is not None else (
699 np.ones(target.extras.n, dtype=bool)),
700 'match_row': target_row_match,
701 }
703 for (columns, out_original, out_matched, in_original, in_matched, matches, name_cat) in (
704 (
705 self.config.columns_ref_copy,
706 data_ref,
707 data_target,
708 ref,
709 target,
710 target_row_match,
711 'target',
712 ),
713 (
714 self.config.columns_target_copy,
715 data_target,
716 data_ref,
717 target,
718 ref,
719 ref_row_match,
720 'reference',
721 ),
722 ):
723 matched = matches >= 0
724 idx_matched = matches[matched]
725 logger.info('Matched %d/%d %s sources', np.sum(matched), len(matched), name_cat)
727 for column in columns:
728 values = in_original.catalog[column]
729 out_original[column] = values
730 dtype = in_original.catalog[column].dtype
732 # Pandas object columns can have mixed types - check for that
733 if dtype is object: 733 ↛ 734line 733 didn't jump to line 734 because the condition on line 733 was never true
734 types = list(set((type(x) for x in values)))
735 if len(types) != 1:
736 raise RuntimeError(f'Column {column} dtype={dtype} has multiple types={types}')
737 dtype = types[0]
739 value_fill = default_value(dtype)
741 # Without this, the dtype would be '<U1' for an empty Unicode string
742 if dtype is str: 742 ↛ 743line 742 didn't jump to line 743 because the condition on line 742 was never true
743 dtype = f'<U{max(len(x) for x in values)}'
745 column_match = np.full(in_matched.extras.n, value_fill, dtype=dtype)
746 column_match[matched] = in_original.catalog[column][idx_matched]
747 out_matched[f'match_{column}'] = column_match
749 logger.info(
750 'Completed match disambiguating in %.2fs (total %.2fs) with %d disambiguated',
751 time.process_time() - t_begin,
752 time.process_time() - t_init,
753 n_ambiguous,
754 )
756 catalog_out_ref = astropy.table.Table(data_ref)
757 for column, description in {
758 'match_candidate': 'Whether the object was selected as a candidate for matching',
759 'match_row': 'The index of the best matched row in the target table, if any',
760 'match_count': 'The number of candidate matching target objects, i.e. those within the match'
761 ' distance but excluding objects already matched to a ref object',
762 'match_chi2': 'The sum of all finite reduced chi-squared values over all match columns',
763 'match_n_chi2_finite': 'The number of match columns with finite chi-squared',
764 }.items():
765 catalog_out_ref[column].description = description
766 catalog_out_target = astropy.table.Table(data_target)
768 return catalog_out_ref, catalog_out_target, exceptions