Fix destination reachability and add shared-basis evaluation metrics

CHANGELOG.md

@@ -6,6 +6,39 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
 ## [Unreleased]
+### Fixed
+- Link-state routing now treats a ground station as reachable through any
+  satellite above its horizon, choosing whichever minimises path length plus
+  GSL length. The Hypatia-derived code path collapsed visibility to the single
+  nearest satellite before routing, which made the destination a fixed
+  satellite rather than the ground station. Under sparse connectivity that
+  satellite is often in an unreachable component while another visible one is
+  reachable, so link-state reported failure on pairs that were deliverable:
+  in a Ring topology it delivered 80 of 408 deliverable pairs instead of all
+  of them.
+- Path stretch is now also reported against an algorithm-independent baseline.
+  The previous baseline was a shortest path to whichever egress satellite the
+  algorithm being measured happened to reach, so each algorithm was graded
+  against a different target and over a different subset of pairs, and neither
+  the subset nor its size was reported.
+### Added
+- Delivery accounting per snapshot (`delivery_*`): deliverable pairs, delivered
+  pairs, delivery rate, forwarding failures, and the separate causes of
+  non-delivery (no source visibility, no destination visibility, graph
+  disconnection).
+- `stretch_hop_shared` and `stretch_dist_shared`, graded against the best
+  end-to-end route to any satellite the destination can see.
+- `delivery_non_optimal_egress_rate`, the share of delivered pairs that exited
+  through a satellite other than the optimal one.
+- `scripts/run-matrix-parallel.sh` for running an evaluation matrix as parallel
+  Docker jobs with per-job wall-clock recorded.
+### Removed
+- `predictive_link_state` and `traditional_segment_routing`, neither of which
+  was used by any published result. The former was link-state evaluated on a
+  future snapshot with no update suppression, so it could only lose on every
+  axis measured and did not implement the scheme its name suggested.
+- The `--prediction-horizon-minutes` and `--segment-mode` flags, which no
+  remaining algorithm reads.
 
 ## [0.1.4] - 2026-06-18
 ### Fixed

docs/evaluation.md

@@ -4,17 +4,38 @@ LEOPath evaluation focuses on path optimality, stability, and state size under d
 
 ## Metrics
 
-- **Stretch**: hop and distance ratio vs shortest-path baseline.
+- **Delivery**: how many ground-station pairs were physically deliverable in a snapshot, and how many the algorithm actually delivered.
+- **Stretch**: hop and distance ratio against a shortest-path baseline, reported on two bases (see below).
 - **Churn**: next-hop changes between consecutive snapshots.
 - **Memory footprint**: forwarding state size per satellite.
 - **Compute time**: wall-clock time to compute routing state per time step.
 
 Why these metrics:
 
-- Stretch captures path optimality cost of reduced state.
+- Delivery gives stretch a denominator. Without it, an algorithm that fails on the hard pairs is rewarded, because only its successes reach the average.
+- Stretch captures the path-optimality cost of reduced state.
 - Churn indicates update frequency and control-plane overhead.
-- Memory footprint reflects routing table scalability.
-- Compute time serves as a practical proxy for algorithmic complexity.
+- Memory footprint reflects routing-table scalability.
+- Compute time is a rough proxy for algorithmic complexity on the host, not a measure of on-board forwarding cost.
+
+### Reachability and the stretch baseline
+
+Reachability is decided from the topology before any algorithm runs, so every algorithm is scored over the same pairs. Each ordered ground-station pair in a snapshot falls into one of five buckets, reported as `delivery_*` columns:
+
+| Bucket | Meaning |
+| --- | --- |
+| `no_src_visibility` | the source ground station sees no satellite |
+| `no_dst_visibility` | the destination ground station sees no satellite |
+| `disconnected` | no ISL path reaches any satellite the destination can see |
+| `deliverable` | a path exists, so the pair counts toward `delivery_rate` |
+| `delivered` | the algorithm got a packet there; the shortfall is `forwarding_failure` |
+
+A ground station is reachable through **any** satellite above its horizon, not only its nearest one. The baseline for stretch is therefore the best end-to-end route to any of them, which makes it identical for every algorithm. Two stretch families are written:
+
+- `stretch_hop` / `stretch_dist` grade an algorithm against a shortest path to whichever egress satellite it happened to reach. An algorithm that delivers through a poor egress still scores near 1.0, because the baseline follows it there. Kept for continuity with earlier runs.
+- `stretch_hop_shared` / `stretch_dist_shared` grade every algorithm against the same lower bound. Use these for comparisons between algorithms.
+
+`delivery_non_optimal_egress_rate` reports how often an algorithm delivered through an egress other than the optimal one, which is what separates the two families. A shortest-path algorithm scores 1.000000 on the shared basis by construction, so link-state doubles as a correctness check on the metric itself.
 
 Optional metrics to add later:
 

leopath/experiments/eval_harness.py

@@ -326,6 +326,9 @@ def run_evaluation(
                 **flatten_distribution("srv6_srh_bytes", explicit_srv6_srh_stats),
                 **flatten_distribution("stretch_hop", stretch_stats["hop"]),
                 **flatten_distribution("stretch_dist", stretch_stats["distance"]),
+                **flatten_distribution("stretch_hop_shared", stretch_stats["hop_shared"]),
+                **flatten_distribution("stretch_dist_shared", stretch_stats["distance_shared"]),
+                **{f"delivery_{key}": value for key, value in stretch_stats["delivery"].items()},
                 **{
                     f"explicit_failover_{key}": value
                     for key, value in explicit_failover_stats.items()

leopath/experiments/metrics.py

@@ -648,22 +648,55 @@ def compute_path_stretch(
     route_plans = route_plans or {}
     sat_set = set(satellite_ids)
     sat_graph = topology_graph.subgraph(satellite_ids)
-    hop_stretches = []
-    dist_stretches = []
+
+    # Legacy basis: each algorithm graded against a shortest path to whichever
+    # egress it happened to reach. Retained for continuity with earlier runs.
+    hop_stretches: list[float] = []
+    dist_stretches: list[float] = []
+    # Shared basis: every algorithm graded against the same lower bound, the
+    # best end-to-end route to any satellite the destination can see.
+    shared_hop_stretches: list[float] = []
+    shared_dist_stretches: list[float] = []
+
+    total_pairs = 0
+    no_src_visibility = 0
+    no_dst_visibility = 0
+    disconnected = 0
+    deliverable = 0
+    delivered = 0
+    non_optimal_egress = 0
+
     for src_index, src_gs_id in enumerate(ground_station_ids):
         for dst_index, dst_gs_id in enumerate(ground_station_ids):
             if src_gs_id == dst_gs_id:
                 continue
+            total_pairs += 1
+
             src_sat, src_gsl_dist = attachments[src_index]
-            if src_sat is None:
-                continue
-            if src_sat not in sat_set:
+            if src_sat is None or src_sat not in sat_set:
+                no_src_visibility += 1
                 continue
+
             destination_visibility = (
                 None
                 if ground_station_satellites_in_range is None
                 else ground_station_satellites_in_range[dst_index]
             )
+            if ground_station_satellites_in_range is not None and not destination_visibility:
+                no_dst_visibility += 1
+                continue
+
+            # Reachability is decided from the topology alone, before any
+            # algorithm is consulted, so every algorithm is scored over the
+            # same set of pairs and against the same lower bound.
+            best_hop_sat, best_hops_total, best_dist_sat, best_dist_total = _best_reachable_egress(
+                sat_graph, src_sat, src_gsl_dist, destination_visibility
+            )
+            if best_dist_total is None:
+                disconnected += 1
+                continue
+            deliverable += 1
+
             dst_sat = _resolve_routed_destination_satellite(
                 fstate,
                 topology_graph,
@@ -675,7 +708,7 @@ def compute_path_stretch(
                 destination_visibility,
             )
             if dst_sat is None:
-                continue
+                continue  # deliverable, but this algorithm failed to deliver
             dst_gsl_dist = _lookup_visible_satellite_distance(destination_visibility, dst_sat)
             if dst_gsl_dist is None:
                 _, nearest_dst_gsl_dist = attachments[dst_index]
@@ -704,13 +737,36 @@ def compute_path_stretch(
             if algo_hops is None or algo_dist is None:
                 continue
 
+            delivered += 1
+            if dst_sat != best_dist_sat:
+                non_optimal_egress += 1
+
             if opt_hops_total > 0:
                 hop_stretches.append(algo_hops / opt_hops_total)
             if opt_dist_total > 0.0:
                 dist_stretches.append(algo_dist / opt_dist_total)
+            if best_hops_total:
+                shared_hop_stretches.append(algo_hops / best_hops_total)
+            if best_dist_total > 0.0:
+                shared_dist_stretches.append(algo_dist / best_dist_total)
+
     return {
         "hop": summarize_distribution(hop_stretches),
         "distance": summarize_distribution(dist_stretches),
+        "hop_shared": summarize_distribution(shared_hop_stretches),
+        "distance_shared": summarize_distribution(shared_dist_stretches),
+        "delivery": {
+            "total_pairs": float(total_pairs),
+            "no_src_visibility": float(no_src_visibility),
+            "no_dst_visibility": float(no_dst_visibility),
+            "disconnected": float(disconnected),
+            "deliverable": float(deliverable),
+            "delivered": float(delivered),
+            "forwarding_failure": float(deliverable - delivered),
+            "delivery_rate": (delivered / deliverable) if deliverable else 0.0,
+            "non_optimal_egress": float(non_optimal_egress),
+            "non_optimal_egress_rate": ((non_optimal_egress / delivered) if delivered else 0.0),
+        },
     }
 
 
@@ -770,6 +826,49 @@ def _resolve_routed_destination_satellite(
     return None
 
 
+def _best_reachable_egress(
+    sat_graph: nx.Graph,
+    src_sat: int,
+    src_gsl_dist: float,
+    destination_visibility: list[tuple[float, int]] | None,
+) -> tuple[int | None, float | None, int | None, float | None]:
+    """Best end-to-end route to a ground station over any satellite it can see.
+
+    A ground station is reachable through any satellite currently above its
+    horizon, so the lower bound on an end-to-end path is the minimum over
+    those satellites, not the route to whichever one an algorithm happened
+    to pick. Returns the hop-optimal and distance-optimal egress separately,
+    since the two need not be the same satellite.
+
+    Returns ``(hop_egress, hop_total, dist_egress, dist_total)``; the entries
+    are ``None`` when no visible satellite is reachable from ``src_sat``.
+    """
+    if not destination_visibility:
+        return None, None, None, None
+
+    best_hops: float | None = None
+    best_hop_sat: int | None = None
+    best_dist: float | None = None
+    best_dist_sat: int | None = None
+
+    for gsl_dist, candidate_sat in destination_visibility:
+        hops, dist = _shortest_path_lengths(sat_graph, src_sat, candidate_sat)
+        if hops is None or dist is None:
+            continue
+        # Both GSL legs are part of the end-to-end path, so the baseline
+        # counts them exactly as the delivered path does.
+        total_hops = hops + 2
+        total_dist = float(dist) + float(src_gsl_dist) + float(gsl_dist)
+        if best_hops is None or total_hops < best_hops:
+            best_hops = total_hops
+            best_hop_sat = candidate_sat
+        if best_dist is None or total_dist < best_dist:
+            best_dist = total_dist
+            best_dist_sat = candidate_sat
+
+    return best_hop_sat, best_hops, best_dist_sat, best_dist
+
+
 def _shortest_path_lengths(
     sat_graph: nx.Graph,
     src_sat: int,

leopath/experiments/plot_seam_robustness.py

@@ -80,8 +80,11 @@ def plot_delivery(eval_data: Path, output_dir: Path) -> None:
     width = 0.38
     fig, ax = plt.subplots(figsize=(8, 4.5))
     bars_pivot = ax.bar(
-        [i - width / 2 for i in x], pivot, width,
-        label="Topological routing with pivot (ours)", color=PIVOT_COLOR
+        [i - width / 2 for i in x],
+        pivot,
+        width,
+        label="Topological routing with pivot (ours)",
+        color=PIVOT_COLOR,
     )
     bars_dra = ax.bar(
         [i + width / 2 for i in x], dra, width, label=r"Hop-only ($\delta_{hop}$)", color=DRA_COLOR
@@ -191,8 +194,16 @@ def pos(p, s):
     px = [pos(p, s)[0] for p, s in pivot_nodes]
     py = [pos(p, s)[1] for p, s in pivot_nodes]
     ax.plot(px, py, color=PIVOT_COLOR, lw=3, zorder=2)
-    ax.scatter([dx], [dy], marker="o", s=120, facecolors="none", edgecolors=PIVOT_COLOR,
-               zorder=5, linewidths=3)
+    ax.scatter(
+        [dx],
+        [dy],
+        marker="o",
+        s=120,
+        facecolors="none",
+        edgecolors=PIVOT_COLOR,
+        zorder=5,
+        linewidths=3,
+    )
     ax.text(
         pos(2, 2)[0],
         pos(0, 2)[1] + 0.45,

leopath/network_state/routing_algorithms/shortest_path_link_state_routing/fstate_calculation.py

@@ -19,6 +19,7 @@ def calculate_fstate_shortest_path_object_no_gs_relay(
     ground_stations: list[GroundStation],
     gsl_attachment_strategy: GSLAttachmentStrategy,
     current_time: Time,
+    ground_station_satellites_in_range: list | None = None,
 ) -> dict:
     """
     Calculates forwarding state using shortest paths over ISLs only (no GS relays).
@@ -28,23 +29,34 @@ def calculate_fstate_shortest_path_object_no_gs_relay(
         ground_stations: List of ground stations
         gsl_attachment_strategy: Strategy for selecting GSL attachments
         current_time: Current simulation time for satellite positioning
+        ground_station_satellites_in_range: Full per-GS visibility list. When
+            given, a ground station is treated as reachable through any
+            satellite currently above its horizon, and the routing chooses
+            whichever of those minimises path length plus GSL length.
+
+    Passing only a single attachment per ground station, which is what the
+    Hypatia-derived code path did unconditionally, makes the destination a
+    fixed satellite rather than the ground station itself. Under sparse
+    connectivity that satellite is frequently in an unreachable component
+    even when another visible satellite is reachable, so the algorithm
+    reports failure for pairs that are in fact deliverable.
     """
     log.debug("Calculating shortest path fstate object (no GS relay)")
 
-    # Use the GSL attachment strategy to compute visibility
-    gsl_attachments = gsl_attachment_strategy.select_attachments(
-        topology_with_isls, ground_stations, current_time
-    )
-
-    # Convert single attachments to the expected format for compatibility
-    # TODO: Refactor the routing algorithm to work directly with single attachments
-    ground_station_satellites_in_range = []
-    for gs_idx, (distance, sat_id) in enumerate(gsl_attachments):
-        if sat_id != -1:  # Valid attachment
-            ground_station_satellites_in_range.append([(distance, sat_id)])
-        else:  # No attachment found
-            ground_station_satellites_in_range.append([])
-            log.warning(f"Ground station {gs_idx} has no satellite attachment")
+    if ground_station_satellites_in_range is None:
+        # Fallback: collapse to the single nearest attachment per ground
+        # station. Kept only so the function remains usable without a
+        # precomputed visibility list; it understates reachability.
+        gsl_attachments = gsl_attachment_strategy.select_attachments(
+            topology_with_isls, ground_stations, current_time
+        )
+        ground_station_satellites_in_range = []
+        for gs_idx, (distance, sat_id) in enumerate(gsl_attachments):
+            if sat_id != -1:  # Valid attachment
+                ground_station_satellites_in_range.append([(distance, sat_id)])
+            else:  # No attachment found
+                ground_station_satellites_in_range.append([])
+                log.warning(f"Ground station {gs_idx} has no satellite attachment")
 
     full_graph = topology_with_isls.graph
     sat_neighbor_to_if = topology_with_isls.sat_neighbor_to_if

leopath/network_state/routing_algorithms/shortest_path_link_state_routing/one_iface_free_bw_allocation_only_over_isls.py

@@ -40,6 +40,7 @@ def algorithm_free_one_only_over_isls(
     gsl_attachment_strategy: GSLAttachmentStrategy,
     current_time: Time,
     list_gsl_interfaces_info: list,  # Info about bandwidth per node/interface
+    ground_station_satellites_in_range: list | None = None,
 ) -> dict:
     """
     Calculates bandwidth and forwarding state (shortest paths via ISLs only, no GS relaying)
@@ -71,7 +72,11 @@ def algorithm_free_one_only_over_isls(
         constellation_data, ground_stations, list_gsl_interfaces_info
     )
     fstate = _calculate_forwarding_state(
-        topology_with_isls, ground_stations, gsl_attachment_strategy, current_time
+        topology_with_isls,
+        ground_stations,
+        gsl_attachment_strategy,
+        current_time,
+        ground_station_satellites_in_range,
     )
 
     return {
@@ -121,6 +126,7 @@ def _calculate_forwarding_state(
     ground_stations: list[GroundStation],
     gsl_attachment_strategy: GSLAttachmentStrategy,
     current_time: Time,
+    ground_station_satellites_in_range: list | None = None,
 ) -> dict:
     """
     Returns the forwarding state object using shortest path calculation.
@@ -131,6 +137,7 @@ def _calculate_forwarding_state(
             ground_stations,
             gsl_attachment_strategy,
             current_time,
+            ground_station_satellites_in_range,
         )
         log.debug("Calculated forwarding state object.")
         return fstate

leopath/network_state/routing_algorithms/shortest_path_link_state_routing/shortest_path_link_state_routing.py

@@ -38,6 +38,10 @@ def compute_state(
         epoch = Time("2000-01-01 00:00:00", scale="tdb")
         current_time = epoch + time_since_epoch_ns * astro_units.ns
 
+        # Route toward the ground station, not toward one chosen satellite:
+        # any satellite currently above the destination's horizon is a valid
+        # egress, and the fstate calculation picks whichever minimises path
+        # length plus GSL length.
         return algorithm_free_one_only_over_isls(
             time_since_epoch_ns,
             constellation_data,
@@ -46,4 +50,5 @@ def compute_state(
             gsl_strategy,
             current_time,
             list_gsl_interfaces_info,
+            ground_station_satellites_in_range,
         )