diff --git a/src/retina_analytics/association.py b/src/retina_analytics/association.py index 6b5752c..84af0f9 100644 --- a/src/retina_analytics/association.py +++ b/src/retina_analytics/association.py @@ -214,6 +214,39 @@ def _bistatic_delay_at(target_enu, tx_enu, rx_enu=(0, 0, 0), d_bl=None): # wrong, so they abstain too. _MIN_BISECTOR = 0.15 +# Below this implied speed the heading of an implied velocity is noise, so the +# velocity-conflict test compares magnitudes only. Well under the 120 m/s +# floor of the commercial envelope, so no real target is ever exempted by it. +_MIN_HEADING_SPEED_MS = 30.0 + +# How close two pairings' grid positions must be to be merged into one solver +# input. One association grid step, not two: the old 6.0 reached a full cell +# beyond a true cluster's own spread, and two aircraft 5 km apart merged into a +# single candidate the position solver then answered by splitting the +# difference. Swept on the offline bench (15-node ring, 6 seeds, --mode track +# --cv-fit deferred), and monotone across the whole range — every metric +# improves as the radius shrinks, with the real-track count identical at each +# point: +# +# published contam. foreign nodes/solve ghost by track +# 6.0 32.0% 0.76 55.6% +# 4.5 28.2% 0.63 54.3% +# 3.0 25.1% 0.52 50.0% +# +# 3.0 is the floor of what was swept rather than a measured optimum; the trend +# says a smaller radius is worth probing, but below the grid step two pairings +# of one aircraft start landing in cells that can no longer reach each other. +_MERGE_DIST_KM = 3.0 + +# Velocity-conflict thresholds for pair-level exclusivity (see +# velocity_conflict). Two pairings sharing a track whose implied velocities +# differ by more than 80 m/s in magnitude, or 40° in heading, cannot both +# describe it. Both sit well outside what the inference's own error can +# produce on a true pairing (median 4° of heading error, measured) and well +# inside the separation a cross-aircraft pairing shows. +_PAIR_VEL_DV_MS = 80.0 +_PAIR_VEL_DTHETA_DEG = 40.0 + def _bisector_from_ranges(target_enu, tx_enu, rx_enu, d_tx, d_rx): """b = u_tx + u_rx at the target, from ranges that are already in hand.""" @@ -265,6 +298,52 @@ def implied_horizontal_speed(m_a, b_a, m_b, b_b): return None if v is None else math.hypot(v[0], v[1]) +def _km_between(p, q) -> float: + """Ground distance between two candidates' positions, in km. + + Flat-earth, the same approximation format_track_pairs_for_solver's + vectorised distance matrix uses — over a merge radius of a few km the + difference from a great circle is centimetres. + """ + return math.hypot( + (p.lat - q.lat) * KM_PER_DEG_LAT, + (p.lon - q.lon) * km_per_deg_lon(0.5 * (p.lat + q.lat)), + ) + + +def velocity_conflict(v_a, v_b, dv_ms: float, dtheta_deg: float) -> bool: + """Whether two implied velocities are too different to be one aircraft. + + Deliberately asymmetric in what it proves. Agreement proves nothing — + two Doppler projections always admit *some* level-flight velocity (see + implied_horizontal_velocity), so a crossed pairing produces one as + readily as a true one, which is why the implied *speed* measured 0% power + as a standalone rejection. Disagreement between two hypotheses that + share a track is a different question: one track is one aircraft and an + aircraft has one velocity, so if the two claims about it cannot both be + true, at most one pairing is. + + Returns False whenever either side has no inference to offer — "no + information" is not evidence of conflict. + """ + if v_a is None or v_b is None: + return False + s_a, s_b = math.hypot(v_a[0], v_a[1]), math.hypot(v_b[0], v_b[1]) + if abs(s_a - s_b) > dv_ms: + return True + # A near-zero implied speed has no meaningful heading — the direction is + # then all noise — so heading is only asked once both sides have enough + # speed for the angle to mean something. + if min(s_a, s_b) < _MIN_HEADING_SPEED_MS: + return False + ang = abs( + math.degrees( + math.atan2(v_a[0] * v_b[1] - v_a[1] * v_b[0], v_a[0] * v_b[0] + v_a[1] * v_b[1]), + ) + ) + return ang > dtheta_deg + + # ── Node Pair Configuration ───────────────────────────────────────────────── @@ -453,6 +532,12 @@ class TrackPairCandidate: doppler_b: float snr_a: float snr_b: float + # Sample time (epoch seconds) of each of those latest measurements. The + # two nodes sample at independent phases, so these differ by up to a frame + # interval and the solver — which treats its measurement set as one epoch — + # needs them to dead-reckon each delay onto a common time. + t_s_a: float + t_s_b: float # Fitted trajectory at the most recent epoch. lat: float lon: float @@ -464,6 +549,18 @@ class TrackPairCandidate: chi2_per_dof: float | None = None dof: int = 0 n_epochs: int = 0 + # |predicted - measured| delay summed over both nodes at the grid point the + # coarse gate matched, in µs. This is the ordering _pair_tracks already + # sorts its hypotheses by; carried on the candidate because the deferred + # path (cv_fit=None — production) has no chi2, so when two pairings claim + # the same track this is the only ranking in hand. It is a tie-break, not + # a test: at n=2 it is ~0 for a crossed pairing too. + grid_resid_us: float = 0.0 + # (v_east, v_north) in m/s that the two Dopplers imply at that grid point + # for level flight, or None where the geometry could not support the + # inference. On the deferred path vel_east_ms/vel_north_ms hold the same + # numbers; this field keeps "no inference" distinguishable from "zero". + implied_vel: tuple | None = None # The measurements the constant-velocity fit needs, in # fit_constant_velocity's input shape. Carried so the fit can run # somewhere other than here: it is an ~86 ms LM solve and submit_tracks @@ -929,6 +1026,19 @@ def _has_receiver_position(config: dict) -> bool: return not (_coord(config, "rx_lat") == 0.0 and _coord(config, "rx_lon") == 0.0) +def _worlds_compatible(world_a, world_b) -> bool: + """Whether two world tags may be associated with each other. + + Fail-open in both directions, which is the rule the ADS-B seed gate already + applies: an unknown world on either side — no provider, a node the provider + does not know, an untagged state — is compatible with everything, and only + two *known* and different tags are refused. A tag that is present but empty + is not a world; it reads as unknown rather than as a third world that + matches nothing. + """ + return not world_a or not world_b or world_a == world_b + + def _merge_epochs_multi(histories: list[tuple[str, list[dict]]]) -> list: """N-node generalisation of _merge_epochs, same output shape. @@ -1017,6 +1127,11 @@ def __init__( adsb_seed_max_dr_age_s: float = ADSB_SEED_MAX_DR_AGE_S, adsb_provider=None, node_world_provider=None, + merge_dist_km: float = _MERGE_DIST_KM, + pair_vel_exclusive: bool = True, + merge_vel_consistent: bool = True, + pair_vel_dv_ms: float = _PAIR_VEL_DV_MS, + pair_vel_dtheta_deg: float = _PAIR_VEL_DTHETA_DEG, ): self.delay_gate_us = delay_gate_us self.doppler_gate_hz = doppler_gate_hz @@ -1067,6 +1182,22 @@ def __init__( # stage 2 for why this exists and why greedy-on-chi2 is safe where two # earlier exclusivity schemes were not. self.cv_exclusive = cv_exclusive + # Pair-level exclusivity for the DEFERRED path only (cv_fit is None, + # i.e. production), where cv_exclusive above has nothing to rank on and + # every hypothesis the coarse grid passed is emitted. See + # _drop_velocity_conflicts. Gated so a caller that does supply a fit + # keeps exactly the chi2 arbitration it had. + self.pair_vel_exclusive = pair_vel_exclusive + self.pair_vel_dv_ms = pair_vel_dv_ms + self.pair_vel_dtheta_deg = pair_vel_dtheta_deg + # How far apart two pairings may sit and still be merged into one + # solver input — see _MERGE_DIST_KM — and whether proximity alone + # earns the merge or their implied velocities must agree too. Both + # are separable from pair_vel_exclusive above because they act at a + # different stage (clustering, not hypothesis pruning) and were + # measured separately on the bench. + self.merge_dist_km = merge_dist_km + self.merge_vel_consistent = merge_vel_consistent # node_id -> (bearing -> observed limit km | None), or None for no # constraint. Injected the same way cv_fit is, so this library keeps # knowing nothing about NodeAnalyticsManager. @@ -1089,6 +1220,12 @@ def __init__( self.track_pairs_unfitted: int = 0 # too few epochs, or the fit failed self.track_pairs_superseded: int = 0 # lost their tracks to a better fit self.track_pairs_deferred: int = 0 # rounds cut short by a budget (once per round) + # Position clusters that turned out to hold two different tracks of one + # node and were split into one solver input each — see + # format_track_pairs_for_solver. Counted once per cluster split, not + # once per input emitted, so it reads as "how often two aircraft were + # about to be solved as one". + self.cluster_splits: int = 0 # Adjacency index: node_id → set of neighbor node_ids that share a real # overlap zone (delay_pairs is non-empty). Built during registration so # submit_frame can iterate O(K) neighbors instead of O(N) all nodes. @@ -1230,6 +1367,15 @@ def __init__( # hex collision, a mis-tag) whatever its residuals say. Optional and # fail-open in both directions: no provider, or an unknown world on # either side, changes nothing. + # + # The same reasoning applies one level lower, to the overlap zones + # bottom-up pairing is drawn from: two nodes in different worlds can + # never see the same echo, so a grid between them only ever pairs a + # synthetic tracklet with a real one. On the test fleet that is 400 + # of 1653 zones (50 synthetic × 8 hardware nodes), 39 of them with a + # real overlap, and the real node ids duly turned up in a third of + # the synthetic dark solver records. register_node and + # rebuild_zones_for therefore skip cross-world pairs outright. self.node_world_provider = node_world_provider # Counters — plain += like the claiming ones above. self.adsb_seed_rounds: int = 0 @@ -1239,6 +1385,12 @@ def __init__( self.adsb_seed_world_rejects: int = 0 self.adsb_tracklets_excluded: int = 0 self.adsb_inputs_emitted: int = 0 + # Node pairs registration or a rebuild declined to build a grid for + # because the two nodes are in different worlds. Counted per pair + # considered, so it rises by O(N) on every registration in a mixed + # fleet — the useful reading is that it is non-zero at all, and + # against overlap_zones, how much of the pair space it is removing. + self.assoc_world_skipped_pairs: int = 0 def register_node(self, node_id: str, config: dict): """Register a node and pre-compute overlap zones with all existing nodes. @@ -1250,7 +1402,9 @@ def register_node(self, node_id: str, config: dict): as long as their geometry (RX/TX position) hasn't changed. A node whose config carries no receiver position is registered but takes - no part in overlap — see _has_receiver_position. + no part in overlap — see _has_receiver_position. A pair whose two nodes + are in known and different worlds gets no zone either — see + node_world_provider and _worlds_compatible. """ positioned = _has_receiver_position(config) rx_alt_km = (config.get("rx_alt_ft") or 0) * 0.3048 / 1000.0 @@ -1310,10 +1464,19 @@ def register_node(self, node_id: str, config: dict): # Pre-compute overlap zones with existing nodes (serialised to avoid # RuntimeError: dictionary changed size during iteration when multiple # nodes register concurrently from a thread-pool executor). + my_world = self._node_world(node_id) for existing_id, existing_geo in list(self.node_geometries.items()): if not self._is_positioned(existing_id): continue pair_key = tuple(sorted([node_id, existing_id])) + if not _worlds_compatible(my_world, self._node_world(existing_id)): + # Re-registration can be what moved this node between + # worlds, so drop rather than merely skip: a grid built + # while the pair was compatible must not survive the + # change that made it cross-world. + self._drop_pair(pair_key) + self.assoc_world_skipped_pairs += 1 + continue zone = compute_overlap_zone( geo if pair_key[0] == node_id else existing_geo, existing_geo if pair_key[0] == node_id else geo, @@ -1348,6 +1511,7 @@ def _reset_for_tests(self) -> None: "track_pairs_unfitted", "track_pairs_superseded", "track_pairs_deferred", + "cluster_splits", "claim_rounds", "claims_matched", "claim_conflicts", @@ -1360,6 +1524,7 @@ def _reset_for_tests(self) -> None: "adsb_seed_world_rejects", "adsb_tracklets_excluded", "adsb_inputs_emitted", + "assoc_world_skipped_pairs", ): setattr(self, name, 0) @@ -1367,6 +1532,27 @@ def _is_positioned(self, node_id: str) -> bool: """Whether a registered node has a receiver position to pair against.""" return _has_receiver_position(self.node_configs.get(node_id, {})) + def _node_world(self, node_id: str): + """This node's world, or None when nothing can say. + + No provider means no world gate at all, which is what every caller + that never injects one gets. + """ + if self.node_world_provider is None: + return None + return self.node_world_provider(node_id) + + def _drop_pair(self, pair_key: tuple[str, str]) -> None: + """Remove one pair's overlap zone and both adjacency entries. + + Caller holds _register_lock. Both nodes stay registered and stay + paired with everyone else; only this edge goes. + """ + a_id, b_id = pair_key + self.overlap_zones.pop(pair_key, None) + self._neighbors.get(a_id, set()).discard(b_id) + self._neighbors.get(b_id, set()).discard(a_id) + def _drop_zones_for(self, node_id: str) -> int: """Remove every overlap zone and adjacency entry naming this node. @@ -1427,10 +1613,18 @@ def rebuild_zones_for(self, node_id: str) -> int: geo.fov = self.fov_provider(node_id) rebuilt = 0 with self._register_lock: + my_world = self._node_world(node_id) for other_id, other_geo in list(self.node_geometries.items()): if other_id == node_id or not self._is_positioned(other_id): continue pair_key = tuple(sorted([node_id, other_id])) + if not _worlds_compatible(my_world, self._node_world(other_id)): + # Dropped, not skipped: a rebuild is exactly where a node + # that has since changed world sheds the grids it built + # against the world it left. + self._drop_pair(pair_key) + self.assoc_world_skipped_pairs += 1 + continue a, b = (geo, other_geo) if pair_key[0] == node_id else (other_geo, geo) zone = compute_overlap_zone( a, @@ -1518,9 +1712,7 @@ def _adsb_seed_round( if st is None or st.get("lat") is None or st.get("lon") is None: self.adsb_seed_no_state += 1 continue - st_world = st.get("world") - node_world = world_of.get(nid) - if st_world is not None and node_world is not None and st_world != node_world: + if not _worlds_compatible(st.get("world"), world_of.get(nid)): self.adsb_seed_world_rejects += 1 continue dt = float(last["t_s"]) - st.get("timestamp_ms", 0) / 1000.0 @@ -1616,6 +1808,7 @@ def _adsb_seed_round( "delay_us": float(m["last"]["delay_us"]), "doppler_hz": float(m["last"]["doppler_hz"]), "snr": float(m["last"].get("snr", 0.0)), + "t_s": float(m["last"]["t_s"]), } for m in matches ], @@ -1825,6 +2018,7 @@ def _claim_round( "delay_us": float(m["_last"]["delay_us"]), "doppler_hz": float(m["_last"]["doppler_hz"]), "snr": float(m["_last"].get("snr", 0.0)), + "t_s": float(m["_last"]["t_s"]), } for m in matches ], @@ -2111,15 +2305,17 @@ def _pair_tracks( if self.cv_fit is not None else _b.get("pairs", self._MAX_PAIRS_PER_ROUND) ) - ordered = sorted( - matches.items(), - key=lambda kv: ( - abs(zone._np_pred_a[kv[1]] - float(usable_a[kv[0][0]]["history"][-1]["delay_us"])) - + abs(zone._np_pred_b[kv[1]] - float(usable_b[kv[0][1]]["history"][-1]["delay_us"])) - ), - )[: max(limit, 0)] - - for (i_a, i_b), best_g in ordered: + + def _grid_resid_us(kv) -> float: + (i_a, i_b), g = kv + return float( + abs(zone._np_pred_a[g] - float(usable_a[i_a]["history"][-1]["delay_us"])) + + abs(zone._np_pred_b[g] - float(usable_b[i_b]["history"][-1]["delay_us"])) + ) + + ordered = sorted(((kv, _grid_resid_us(kv)) for kv in matches.items()), key=lambda x: x[1])[: max(limit, 0)] + + for ((i_a, i_b), best_g), grid_resid in ordered: ta, tb = usable_a[i_a], usable_b[i_b] hist_a, hist_b = ta["history"], tb["history"] last_a, last_b = hist_a[-1], hist_b[-1] @@ -2226,6 +2422,8 @@ def _pair_tracks( doppler_b=float(last_b["doppler_hz"]), snr_a=float(last_a.get("snr", 0.0)), snr_b=float(last_b.get("snr", 0.0)), + t_s_a=float(last_a["t_s"]), + t_s_b=float(last_b["t_s"]), lat=lat, lon=lon, alt_km=alt_km, @@ -2235,6 +2433,8 @@ def _pair_tracks( dof=dof, n_epochs=len(epochs), epochs=deferred_epochs, + grid_resid_us=grid_resid, + implied_vel=(vel_seed["vel_east_ms"], vel_seed["vel_north_ms"]) if vel_seed else None, ) (fitted if chi2_per_dof is not None else held).append(cand) @@ -2276,6 +2476,8 @@ def _pair_tracks( results.append(c) else: self.track_pairs_rejected += 1 + if self.cv_fit is None and self.pair_vel_exclusive: + held = self._drop_velocity_conflicts(held) for c in held: # A held pairing has no score, so it cannot claim anything — and a # track already explained by a scored winner does not get to seed a @@ -2286,6 +2488,50 @@ def _pair_tracks( results.append(c) return results + def _drop_velocity_conflicts(self, held: list[TrackPairCandidate]) -> list[TrackPairCandidate]: + """Prune pairings claiming a track another pairing explains differently. + + Stage 2 above arbitrates on chi2, which production never has: it runs + cv_fit=None so nothing is fitted here and every hypothesis the coarse + delay grid passed is emitted, several of them for the same track in a + crowded zone. Those land at similar positions and + format_track_pairs_for_solver merges them, so the solver is handed one + candidate describing two aircraft. + + The evidence already in hand is the Doppler-implied level-flight + velocity at each pairing's matched grid point — the same quantity the + fit is seeded from. It is not usable as a standalone test (see the + seed comment above: 0% power against real cross pairings), and neither + is the coarse delay residual, which an earlier assignment attempt + ranked on and lost recall to because at n=2 it is ~0 for a crossed + pairing too. The *comparison* between two hypotheses about one track + is what neither of them is alone: one track is one aircraft, so two + claims about its velocity that cannot both be true mean at most one + pairing is. Where they agree, nothing has been learned and both + survive — which is why this cannot cost the recall the residual + assignment cost. + + Ties in the residual are broken on ids so a round is deterministic. + """ + winner: dict[tuple[str, str], TrackPairCandidate] = {} + kept: list[TrackPairCandidate] = [] + for c in sorted( + held, + key=lambda p: (p.grid_resid_us, p.node_a_id, p.track_a_id, p.node_b_id, p.track_b_id), + ): + keys = ((c.node_a_id, c.track_a_id), (c.node_b_id, c.track_b_id)) + if any( + velocity_conflict(winner[k].implied_vel, c.implied_vel, self.pair_vel_dv_ms, self.pair_vel_dtheta_deg) + for k in keys + if k in winner + ): + self.track_pairs_superseded += 1 + continue + for k in keys: + winner.setdefault(k, c) + kept.append(c) + return kept + def format_track_pairs_for_solver(self, pairs: list[TrackPairCandidate]) -> list[dict]: """Cluster track pairs by fitted position into multinode solver inputs. @@ -2297,7 +2543,6 @@ def format_track_pairs_for_solver(self, pairs: list[TrackPairCandidate]) -> list if not pairs: return [] - _MERGE_DIST_KM = 6.0 n = len(pairs) parent = list(range(n)) @@ -2312,7 +2557,19 @@ def _find(x: int) -> int: km_per_lat = KM_PER_DEG_LAT km_per_lon = km_per_deg_lon(float(np.mean(lats))) dist_sq = ((lats[:, None] - lats) * km_per_lat) ** 2 + ((lons[:, None] - lons) * km_per_lon) ** 2 - rows, cols = np.where((dist_sq < _MERGE_DIST_KM**2) & (np.arange(n)[:, None] < np.arange(n))) + # Proximity alone was the whole merge criterion, and it is not enough: + # two aircraft crossing within the merge radius are as close together + # as one aircraft's own pairings are, so the union welded them into a + # single candidate. Requiring the two pairings' Doppler-implied + # velocities to be compatible as well separates the crossing case, + # where the headings differ by definition, from the co-located case, + # where they do not. Same abstention rule as everywhere else: a + # pairing with no usable inference blocks nothing. + rows, cols = np.where( + (dist_sq < self.merge_dist_km**2) + & ~self._velocity_conflict_matrix(pairs) + & (np.arange(n)[:, None] < np.arange(n)) + ) for i, j in zip(rows.tolist(), cols.tolist()): parent[_find(i)] = _find(j) @@ -2321,60 +2578,160 @@ def _find(x: int) -> int: groups[_find(i)].append(p) solver_inputs = [] - for group in groups.values(): - by_node: dict[str, dict] = {} - for p in group: - for nid, d, f, s in ( - (p.node_a_id, p.delay_a, p.doppler_a, p.snr_a), - (p.node_b_id, p.delay_b, p.doppler_b, p.snr_b), - ): - if nid not in by_node or s > by_node[nid]["snr"]: - by_node[nid] = {"node_id": nid, "delay_us": d, "doppler_hz": f, "snr": s} - - # Per-node track id sets, so a downstream trim (dropping one - # contaminated node's measurement and re-solving) can also drop - # exactly that node's source tracks from the published solve's - # provenance instead of carrying the whole cluster's track_ids. - track_ids_by_node: dict[str, set] = defaultdict(set) - for p in group: - track_ids_by_node[p.node_a_id].add(p.track_a_id) - track_ids_by_node[p.node_b_id].add(p.track_b_id) - - fitted = [p for p in group if p.chi2_per_dof is not None] - # Worst fit in the cluster, not the best: a cluster is published as - # one target, so a pairing that failed to justify itself should not - # be laundered by a well-fitted neighbour sharing its position. - worst_chi2 = max((p.chi2_per_dof for p in fitted), default=None) - - solver_inputs.append( - { - "initial_guess": { - "lat": sum(p.lat for p in group) / len(group), - "lon": sum(p.lon for p in group) / len(group), - "alt_km": sum(p.alt_km for p in group) / len(group), - }, - "initial_velocity": { - "vel_east_ms": sum(p.vel_east_ms for p in group) / len(group), - "vel_north_ms": sum(p.vel_north_ms for p in group) / len(group), - }, - "measurements": list(by_node.values()), - "n_nodes": len(by_node), - "timestamp_ms": group[0].timestamp_ms, - "adsb_hex": None, - "chi2_per_dof": worst_chi2, - "n_epochs": min(p.n_epochs for p in group), - # Present only when the fit has not run yet: the solver worker - # runs it there, on its own threads and behind its own queue, - # instead of on the frame path. Taken from the pairing with the - # most history, which is the best-conditioned in the cluster. - "cv_epochs": max((p.epochs for p in group if p.epochs), key=len, default=None), - "track_pair_ids": sorted({(p.track_a_id, p.track_b_id) for p in group})[:1], - "track_ids": sorted({p.track_a_id for p in group} | {p.track_b_id for p in group}), - "track_ids_by_node": {nid: sorted(ids) for nid, ids in track_ids_by_node.items()}, - } - ) + for merged in groups.values(): + subs = self._partition_cluster(merged) + if len(subs) > 1: + self.cluster_splits += 1 + solver_inputs.extend(self._solver_input(g) for g in subs) return solver_inputs + def _velocity_conflict_matrix(self, pairs: list[TrackPairCandidate]) -> np.ndarray: + """(n, n) mask: True where two pairings' implied velocities disagree. + + The vectorised form of velocity_conflict — same thresholds, same + abstention when either side has no inference — because this runs + against every pairing pair in the round. + """ + n = len(pairs) + has_v = np.array([p.implied_vel is not None for p in pairs]) + if not self.merge_vel_consistent or not has_v.any(): + return np.zeros((n, n), dtype=bool) + ve = np.array([p.implied_vel[0] if p.implied_vel else 0.0 for p in pairs], dtype=np.float64) + vn = np.array([p.implied_vel[1] if p.implied_vel else 0.0 for p in pairs], dtype=np.float64) + speed = np.hypot(ve, vn) + d_speed = np.abs(speed[:, None] - speed) + ang = np.abs( + np.degrees( + np.arctan2(ve[:, None] * vn - vn[:, None] * ve, ve[:, None] * ve + vn[:, None] * vn), + ) + ) + heading_meaningful = np.minimum(speed[:, None], speed) >= _MIN_HEADING_SPEED_MS + bad = (d_speed > self.pair_vel_dv_ms) | ((ang > self.pair_vel_dtheta_deg) & heading_meaningful) + return bad & has_v[:, None] & has_v + + def _partition_cluster(self, group: list[TrackPairCandidate]) -> list[list[TrackPairCandidate]]: + """Break one connected position cluster into consistent solver inputs. + + The union-find above is transitive, so "within 4.5 km" chains: three + aircraft strung out over 12 km arrive here as one group, and before + this partition ran the whole chain became a single solver input. On + the 15-node bench scene the baseline emitted inputs spanning up to 23 + single-node tracks — an aircraft has one track per node, so a 23-track + input at 10 nodes is describing at least three aircraft as one. + + Two ways a cluster can be describing more than one aircraft, and this + rejects both: + + - A node appears with two different tracks. A node's tracker gives one + track per aircraft, so that is two aircraft outright. The old answer + — keep the node's louder track — is the worst available: the quiet + track's aircraft silently borrows a measurement belonging to the loud + one's, the solver gets a candidate no position explains, and the trim + that follows drops legitimate nodes about as often as contaminated + ones (33 of 69, measured live). + + - The cluster is wider than the merge distance, or holds pairings whose + implied velocities contradict each other. Membership is therefore + tested against EVERY pairing already in the sub-cluster, not just one + of them, which is what bounds the diameter at merge_dist_km instead + of letting it grow with the chain. + + Every resulting sub-cluster is emitted. This stage cannot tell which + aircraft is real — the solver's gates and the resolve slot decide that + downstream — and suppressing the runners-up here is how the SNR pick + went wrong in the first place. + + Best-residual-first, so the strongest hypothesis founds the first + sub-cluster rather than whichever pairing the zone iteration reached + first. + """ + if len(group) == 1: + return [group] + subs: list[list[TrackPairCandidate]] = [] + sub_nodes: list[dict[str, str]] = [] + for p in sorted( + group, + key=lambda c: (c.grid_resid_us, c.node_a_id, c.track_a_id, c.node_b_id, c.track_b_id), + ): + own = ((p.node_a_id, p.track_a_id), (p.node_b_id, p.track_b_id)) + for sub, nodes in zip(subs, sub_nodes): + if any(nodes.get(nid, tid) != tid for nid, tid in own): + continue + if any(_km_between(p, q) >= self.merge_dist_km for q in sub): + continue + if self.merge_vel_consistent and any( + velocity_conflict(p.implied_vel, q.implied_vel, self.pair_vel_dv_ms, self.pair_vel_dtheta_deg) + for q in sub + ): + continue + sub.append(p) + nodes.update(own) + break + else: + subs.append([p]) + sub_nodes.append(dict(own)) + return subs + + def _solver_input(self, group: list[TrackPairCandidate]) -> dict: + """One node-consistent cluster, in the shape the solver worker takes.""" + by_node: dict[str, dict] = {} + for p in group: + for nid, d, f, s, t in ( + (p.node_a_id, p.delay_a, p.doppler_a, p.snr_a, p.t_s_a), + (p.node_b_id, p.delay_b, p.doppler_b, p.snr_b, p.t_s_b), + ): + # First writer wins, and that is not a choice between rivals: + # _split_node_conflicts guarantees every pairing here agrees on + # which track each node contributed, and a track's measurement + # is its own history[-1], so the repeats are the same numbers. + # t_s travels with the measurement it belongs to (see + # TrackPairCandidate.t_s_a): the consumer aligns the cluster's + # delays onto one epoch. + by_node.setdefault(nid, {"node_id": nid, "delay_us": d, "doppler_hz": f, "snr": s, "t_s": t}) + + # Per-node track id sets, so a downstream trim (dropping one + # contaminated node's measurement and re-solving) can also drop + # exactly that node's source tracks from the published solve's + # provenance instead of carrying the whole cluster's track_ids. + track_ids_by_node: dict[str, set] = defaultdict(set) + for p in group: + track_ids_by_node[p.node_a_id].add(p.track_a_id) + track_ids_by_node[p.node_b_id].add(p.track_b_id) + + # Worst fit in the cluster, not the best: a cluster is published as one + # target, so a pairing that failed to justify itself should not be + # laundered by a well-fitted neighbour sharing its position. None + # whenever no pairing here was fitted, which on the deferred path + # (cv_fit=None — production) is always: the fit runs in the solver + # worker, so this field says "not scored yet" rather than "scored 0". + worst_chi2 = max((p.chi2_per_dof for p in group if p.chi2_per_dof is not None), default=None) + + return { + "initial_guess": { + "lat": sum(p.lat for p in group) / len(group), + "lon": sum(p.lon for p in group) / len(group), + "alt_km": sum(p.alt_km for p in group) / len(group), + }, + "initial_velocity": { + "vel_east_ms": sum(p.vel_east_ms for p in group) / len(group), + "vel_north_ms": sum(p.vel_north_ms for p in group) / len(group), + }, + "measurements": list(by_node.values()), + "n_nodes": len(by_node), + "timestamp_ms": group[0].timestamp_ms, + "adsb_hex": None, + "chi2_per_dof": worst_chi2, + "n_epochs": min(p.n_epochs for p in group), + # Present only when the fit has not run yet: the solver worker runs + # it there, on its own threads and behind its own queue, instead of + # on the frame path. Taken from the pairing with the most history, + # which is the best-conditioned in the cluster. + "cv_epochs": max((p.epochs for p in group if p.epochs), key=len, default=None), + "track_pair_ids": sorted({(p.track_a_id, p.track_b_id) for p in group})[:1], + "track_ids": sorted({p.track_a_id for p in group} | {p.track_b_id for p in group}), + "track_ids_by_node": {nid: sorted(ids) for nid, ids in track_ids_by_node.items()}, + } + def get_overlap_summary(self) -> list[dict]: """Return summary of all overlap zones.""" summaries = [] diff --git a/tests/test_track_association.py b/tests/test_track_association.py index b0ccaa3..90fd587 100644 --- a/tests/test_track_association.py +++ b/tests/test_track_association.py @@ -15,6 +15,7 @@ InterNodeAssociator, TrackPairCandidate, _merge_epochs, + velocity_conflict, ) from retina_analytics.constants import KM_PER_DEG_LAT @@ -120,6 +121,38 @@ def _assoc(cv_fit=None, **kw): return a +def _candidate(track_a_id, track_b_id, **kw): + """A TrackPairCandidate at a default position, for the clustering tests. + + Clustering is judged on position, per-node track ids and implied velocity + alone, so these need no zone geometry — building them by hand keeps a + two-aircraft scene readable and exactly reproducible. + """ + fields = dict( + timestamp_ms=1000, + node_a_id="site-a", + node_b_id="site-b", + delay_a=30.0, + delay_b=40.0, + doppler_a=5.0, + doppler_b=-5.0, + snr_a=15.0, + snr_b=15.0, + t_s_a=1.0, + t_s_b=1.0, + lat=34.88, + lon=-82.35, + alt_km=7.0, + vel_east_ms=180.0, + vel_north_ms=-90.0, + implied_vel=(180.0, -90.0), + dof=14, + n_epochs=6, + ) + fields.update(kw) + return TrackPairCandidate(track_a_id=track_a_id, track_b_id=track_b_id, **fields) + + def _cv_fit(): """The real fit, imported lazily so the rest of the module runs without it.""" pytest.importorskip("retina_geolocator") @@ -372,37 +405,240 @@ def test_cluster_reports_its_worst_fit(self): The cluster is published as one target, so its quality is the quality of the weakest pairing holding it together. + + The two pairings here share node-a track a1, which is what makes them + one cluster rather than two: same node, same track, two neighbours — + the 3-node case, not a conflict. It used to be a2/b2 against a1/b1, + which _split_node_conflicts now (correctly) separates into two targets, + so the worst-fit rule would never have been reached. """ - base = dict( - timestamp_ms=1000, - node_a_id="site-a", - node_b_id="site-b", - delay_a=30.0, - delay_b=40.0, - doppler_a=5.0, - doppler_b=-5.0, - snr_a=15.0, - snr_b=15.0, - lat=34.88, - lon=-82.35, - alt_km=7.0, - vel_east_ms=180.0, - vel_north_ms=-90.0, - dof=14, - n_epochs=6, - ) pairs = [ - TrackPairCandidate(track_a_id="a1", track_b_id="b1", chi2_per_dof=0.4, **base), - TrackPairCandidate(track_a_id="a2", track_b_id="b2", chi2_per_dof=9.9, **base), + _candidate("a1", "b1", chi2_per_dof=0.4), + _candidate("a1", "b2", node_b_id="site-c", chi2_per_dof=9.9), ] inputs = InterNodeAssociator().format_track_pairs_for_solver(pairs) assert len(inputs) == 1 assert inputs[0]["chi2_per_dof"] == 9.9 + def test_unfitted_cluster_reports_no_chi2(self): + """Production defers the fit, so the field must say so rather than 0. + + With cv_fit=None nothing here is scored; a numeric chi2 would be the + solver worker's gate reading a quality nobody measured. + """ + inputs = InterNodeAssociator().format_track_pairs_for_solver([_candidate("a1", "b1")]) + assert inputs[0]["chi2_per_dof"] is None + def test_empty_input(self): assert InterNodeAssociator().format_track_pairs_for_solver([]) == [] +class TestClusterPartition: + """A cluster is one aircraft only if nothing in it says otherwise. + + Two things say otherwise: a node appearing with two different tracks (its + tracker gives one track per aircraft, so that is two aircraft), and a + cluster wider than the merge distance (the union-find that builds it is + transitive, so "within 4.5 km" chains across an arbitrary distance). + Before this partition ran, 58-65% of dark candidates measured on the live + fleet carried a node that could not see the aircraft they were published + as. + """ + + def test_a_chain_is_not_one_target(self): + """Three pairings 2 km apart in a line span 4 km — two targets, not one. + + The union edge joins 1-2 and 2-3, so union-find hands the whole chain + over as one group; membership is tested against every pairing already + in a sub-cluster, which is what bounds the diameter. + """ + a = InterNodeAssociator() + step = 2.0 / KM_PER_DEG_LAT + inputs = a.format_track_pairs_for_solver( + [ + _candidate("a1", "b1", lat=34.88), + _candidate("a2", "b2", lat=34.88 + step, node_a_id="site-c", node_b_id="site-d"), + _candidate("a3", "b3", lat=34.88 + 2 * step, node_a_id="site-e", node_b_id="site-f"), + ] + ) + assert len(inputs) == 2 + assert a.cluster_splits == 1 + widest = max(inputs, key=lambda s: s["n_nodes"]) + assert widest["n_nodes"] == 4 + + def test_two_aircraft_split_instead_of_being_merged(self): + """Both nodes see both aircraft, 2 km apart — one input each.""" + a = InterNodeAssociator() + inputs = a.format_track_pairs_for_solver( + [ + _candidate("aP", "bP", lat=34.88, snr_a=20.0, snr_b=20.0), + _candidate("aQ", "bQ", lat=34.88 + 2.0 / KM_PER_DEG_LAT, snr_a=6.0, snr_b=6.0), + ] + ) + assert len(inputs) == 2 + assert a.cluster_splits == 1 + by_node = [s_in["track_ids_by_node"] for s_in in inputs] + assert {"site-a": ["aP"], "site-b": ["bP"]} in by_node + assert {"site-a": ["aQ"], "site-b": ["bQ"]} in by_node + # No SNR pick: the quiet aircraft keeps its own measurements rather + # than borrowing the loud one's. + assert sorted(m["snr"] for s_in in inputs for m in s_in["measurements"]) == [6.0, 6.0, 20.0, 20.0] + + def test_a_false_pairing_between_them_stands_alone(self): + """The cross pairing is emitted too, but never welded onto a true one. + + This stage cannot tell which of three hypotheses is real — the solver's + gates and the resolve slot decide that — so the requirement is only + that no emitted input mixes two aircraft. + """ + a = InterNodeAssociator() + inputs = a.format_track_pairs_for_solver( + [ + _candidate("aP", "bP", lat=34.88, grid_resid_us=0.1), + _candidate("aQ", "bQ", lat=34.88 + 2.0 / KM_PER_DEG_LAT, grid_resid_us=0.1), + # The cross pairing: node a's aircraft P against node b's Q, + # landing between the two true clusters and inside the merge + # radius of both. + _candidate("aP", "bQ", lat=34.88 + 1.0 / KM_PER_DEG_LAT, grid_resid_us=0.9), + ] + ) + assert len(inputs) == 3 + for s_in in inputs: + assert all(len(ids) == 1 for ids in s_in["track_ids_by_node"].values()) + assert {"site-a": ["aP"], "site-b": ["bP"]} in [s["track_ids_by_node"] for s in inputs] + assert {"site-a": ["aQ"], "site-b": ["bQ"]} in [s["track_ids_by_node"] for s in inputs] + + def test_one_aircraft_on_three_nodes_is_not_split(self): + """The legitimate multi-node cluster is untouched. + + Two pairings sharing node-a track a1 are the same aircraft seen by + three nodes — one track per node, no conflict — and must still merge + into a single 3-node solver input. + """ + a = InterNodeAssociator() + inputs = a.format_track_pairs_for_solver( + [ + _candidate("a1", "b1"), + _candidate("a1", "c1", node_b_id="site-c"), + ] + ) + assert len(inputs) == 1 + assert inputs[0]["n_nodes"] == 3 + assert a.cluster_splits == 0 + + def test_velocity_disagreement_blocks_the_merge(self): + """Two coincident pairings heading opposite ways are two targets. + + Proximity alone made this one cluster, which is the crossing case the + old 6 km radius could not distinguish from one aircraft's own + neighbouring pairings. + """ + a = InterNodeAssociator() + inputs = a.format_track_pairs_for_solver( + [ + _candidate("a1", "b1", implied_vel=(200.0, 0.0)), + _candidate("a2", "b2", implied_vel=(-200.0, 0.0), lat=34.881), + ] + ) + assert len(inputs) == 2 + # Split by the merge edge, not by the conflict split, so nothing was + # ever a single cluster to begin with. + assert a.cluster_splits == 0 + + +class TestVelocityExclusivity: + """Deferred-path exclusivity: production has no chi2 to rank on. + + cv_fit is None live, so stage 2's chi2 selection never runs and every + hypothesis the coarse grid passed is emitted — including the several that + claim the same track. The Doppler-implied velocity cannot reject a + pairing on its own (0% power, measured), but two claims about one track's + velocity that cannot both be true mean at most one pairing is. + """ + + def _run(self, **kw): + a = _assoc(**kw) + a.submit_tracks( + "site-a", + [{"track_id": "aP", "history": _history(_NODE_A, 34.88, -82.35, 7.0, 180.0, -90.0, anchor="end")}], + 1000, + ) + return a, a.submit_tracks( + "site-b", + [ + {"track_id": "bP", "history": _history(_NODE_B, 34.88, -82.35, 7.0, 180.0, -90.0, anchor="end")}, + {"track_id": "bQ", "history": _history(_NODE_B, 34.88, -82.35, 7.0, -150.0, 170.0, anchor="end")}, + ], + 2000, + ) + + def test_the_contradicting_hypothesis_is_dropped(self): + """Both pairings pass the coarse gate at the same grid point and the + same delay residual — the residual is ~0 for a crossed pairing at n=2, + which is why an assignment on it alone cost recall. The implied + headings differ by 75°, and that is decisive.""" + a, pairs = self._run() + assert [(p.track_a_id, p.track_b_id) for p in pairs] == [("aP", "bP")] + assert a.track_pairs_superseded == 1 + + def test_off_by_flag_restores_both(self): + a, pairs = self._run(pair_vel_exclusive=False) + assert len(pairs) == 2 + assert a.track_pairs_superseded == 0 + + def test_agreement_keeps_both(self): + """Where the two hypotheses agree, nothing has been learned. + + A three-node target legitimately pairs one track against two + neighbours; this must never become an excuse to drop one of them, and + it is the reason the test is on disagreement rather than on rank. + """ + a = _assoc() + a.submit_tracks( + "site-a", + [{"track_id": "aP", "history": _history(_NODE_A, 34.88, -82.35, 7.0, 180.0, -90.0, anchor="end")}], + 1000, + ) + pairs = a.submit_tracks( + "site-b", + [ + {"track_id": "bP", "history": _history(_NODE_B, 34.88, -82.35, 7.0, 180.0, -90.0, anchor="end")}, + {"track_id": "bP2", "history": _history(_NODE_B, 34.88, -82.35, 7.0, 182.0, -88.0, anchor="end")}, + ], + 2000, + ) + assert len(pairs) == 2 + assert a.track_pairs_superseded == 0 + + def test_a_supplied_fit_keeps_the_chi2_path(self): + """Gated on cv_fit is None, so an inline-fitting caller is untouched.""" + a, pairs = self._run(cv_fit=_cv_fit(), cv_min_span_s=2.0, cv_min_epochs=2) + assert len(pairs) == 1 + assert pairs[0].chi2_per_dof is not None + + +class TestVelocityConflict: + def test_speed_alone_can_decide(self): + assert velocity_conflict((200.0, 0.0), (60.0, 0.0), 80.0, 40.0) + assert not velocity_conflict((200.0, 0.0), (150.0, 0.0), 80.0, 40.0) + + def test_heading_alone_can_decide(self): + assert velocity_conflict((200.0, 0.0), (0.0, 200.0), 80.0, 40.0) + assert not velocity_conflict((200.0, 0.0), (190.0, 40.0), 80.0, 40.0) + + def test_missing_inference_never_conflicts(self): + """No information is not evidence — the abstention this shares with + implied_horizontal_velocity, which returns None on geometry that + cannot support the inference at all.""" + assert not velocity_conflict(None, (200.0, 0.0), 80.0, 40.0) + assert not velocity_conflict((200.0, 0.0), None, 80.0, 40.0) + + def test_near_zero_speeds_compare_magnitude_only(self): + """Below 30 m/s the heading is noise, so opposite directions at a + crawl are not called a conflict.""" + assert not velocity_conflict((5.0, 0.0), (-5.0, 0.0), 80.0, 40.0) + + # Third node at the same receiver — a triple-illuminator site — for exercising # cross-node-pair sharing, which hypothesis selection must never forbid. _NODE_C = { diff --git a/tests/test_world_overlap_gate.py b/tests/test_world_overlap_gate.py new file mode 100644 index 0000000..9cf8694 --- /dev/null +++ b/tests/test_world_overlap_gate.py @@ -0,0 +1,122 @@ +"""Bottom-up pairing does not cross worlds. + +The ADS-B gate added a world check to seeding and claiming, but the overlap +zones pairing draws its candidates from were still built between every pair of +positioned nodes. A synthetic fleet sharing a footprint with real hardware — +the test droplet's 50 simulated nodes over the same city as 8 receivers — +therefore kept a grid for every sim/real pair, and real node ids turned up +inside synthetic dark solves. Two nodes in different worlds can never see the +same echo, so those grids have no true pairing to find. + +Same fail-open rule as the seed gate: only two known, different worlds are +refused. +""" + +from retina_analytics.association import InterNodeAssociator + +# Two receivers a few km apart on the same illuminator: a real, non-empty +# overlap, so "no zone" below is the world gate and not the geometry. +_CFG_A = dict(rx_lat=34.85, rx_lon=-82.40, tx_lat=34.90, tx_lon=-82.30, max_range_km=50, max_bistatic_range_km=50) +_CFG_B = dict(rx_lat=34.86, rx_lon=-82.36, tx_lat=34.90, tx_lon=-82.30, max_range_km=50, max_bistatic_range_km=50) + +_PAIR = ("a", "b") + + +def _register_pair(**kwargs) -> InterNodeAssociator: + assoc = InterNodeAssociator(grid_step_km=5.0, **kwargs) + assoc.register_node("a", dict(_CFG_A)) + assoc.register_node("b", dict(_CFG_B)) + return assoc + + +class TestOverlapZoneWorldGate: + def test_cross_world_pair_gets_no_zone(self): + worlds = {"a": "sim", "b": "real"} + assoc = _register_pair(node_world_provider=worlds.get) + + assert _PAIR not in assoc.overlap_zones + assert assoc.overlap_zones == {} + assert assoc._neighbors.get("a", set()) == set() + assert assoc._neighbors.get("b", set()) == set() + assert assoc.assoc_world_skipped_pairs == 1 + + def test_same_world_pair_is_unchanged(self): + assoc = _register_pair(node_world_provider=lambda nid: "sim") + + assert assoc.overlap_zones[_PAIR].delay_pairs + assert assoc._neighbors["a"] == {"b"} + assert assoc.assoc_world_skipped_pairs == 0 + + def test_one_untagged_node_still_pairs(self): + """Fail-open: an unknown world is compatible with every world, which is + what a node the provider has not seen yet has.""" + worlds = {"a": "real"} + assoc = _register_pair(node_world_provider=worlds.get) + + assert assoc.overlap_zones[_PAIR].delay_pairs + assert assoc._neighbors["a"] == {"b"} + assert assoc.assoc_world_skipped_pairs == 0 + + def test_an_empty_world_tag_is_not_a_world(self): + worlds = {"a": "", "b": "real"} + assoc = _register_pair(node_world_provider=worlds.get) + + assert assoc.overlap_zones[_PAIR].delay_pairs + assert assoc.assoc_world_skipped_pairs == 0 + + def test_no_provider_pairs_everything(self): + assoc = _register_pair() + + assert assoc.overlap_zones[_PAIR].delay_pairs + assert assoc._neighbors["a"] == {"b"} + assert assoc.assoc_world_skipped_pairs == 0 + + def test_a_node_that_changes_world_loses_its_zone_on_rebuild(self): + """The gate has to be able to take a grid away, not only decline to + build one: registration order and a late handshake both mean a pair can + be compatible when it is first built and cross-world afterwards.""" + worlds = {"a": "sim", "b": "sim"} + assoc = _register_pair(node_world_provider=worlds.get) + assert assoc.overlap_zones[_PAIR].delay_pairs + + worlds["b"] = "real" + rebuilt = assoc.rebuild_zones_for("b") + + assert rebuilt == 0 + assert _PAIR not in assoc.overlap_zones + assert assoc._neighbors.get("a", set()) == set() + assert assoc.assoc_world_skipped_pairs == 1 + + def test_rebuild_keeps_same_world_zones(self): + assoc = _register_pair(node_world_provider=lambda nid: "sim") + + rebuilt = assoc.rebuild_zones_for("b") + + assert rebuilt == 1 + assert assoc.overlap_zones[_PAIR].delay_pairs + assert assoc._neighbors["a"] == {"b"} + assert assoc.assoc_world_skipped_pairs == 0 + + def test_a_re_registration_that_changes_world_drops_the_zone(self): + """register_node's own drop path, which rebuild_zones_for cannot cover: + a reconnecting node re-registers with a moved receiver.""" + worlds = {"a": "sim", "b": "sim"} + assoc = _register_pair(node_world_provider=worlds.get) + assert _PAIR in assoc.overlap_zones + + worlds["b"] = "real" + moved = dict(_CFG_B, rx_lat=34.87) + assoc.register_node("b", moved) + + assert _PAIR not in assoc.overlap_zones + assert assoc._neighbors.get("a", set()) == set() + assert assoc.assoc_world_skipped_pairs == 1 + + def test_the_counter_resets_with_the_others(self): + worlds = {"a": "sim", "b": "real"} + assoc = _register_pair(node_world_provider=worlds.get) + assert assoc.assoc_world_skipped_pairs == 1 + + assoc._reset_for_tests() + + assert assoc.assoc_world_skipped_pairs == 0