Address review findings

Author: constantinpape

MIGRATION_GUIDE.md

@@ -2382,11 +2382,11 @@ sources point to their own flat C-order index and background/unreachable voxels
 contain `-1`. Paths are `int64` NumPy-axis-order coordinates from the source to
 the selected target. Target ties are deterministic and use flat C-order indices.
 `number_of_threads=1` retains the specialized sequential heaps; `0` uses
-hardware concurrency. Threaded solves use exact FP64 delta stepping: distance
-fields are bitwise-identical to the sequential result, while an equal-cost
-predecessor or path may differ from the sequential heap but remains deterministic.
-Small workloads and underperforming cost modes automatically retain the
-sequential implementation.
+hardware concurrency. Broad multi-source physical fields can use exact FP64
+delta stepping: distance fields are bitwise-identical to the sequential result,
+while an equal-cost predecessor may differ from the sequential heap but remains
+deterministic. Paths, weighted modes, small workloads, and narrow source sets
+retain the faster specialized heaps.
 
 This Dijkstra field is deliberately distinct from
 `geodesic_distance_field` below. Dijkstra is exact for the chosen 4/8- or
@@ -2401,7 +2401,8 @@ are close on well-sampled, uniform domains but are not numerically equivalent;
 TEASAR therefore uses Dijkstra where it must reconstruct and join paths.
 
 See `development/distance/benchmark_dijkstra.py` for a standalone benchmark of
-the Dijkstra field, predecessor field, weighted field, and early-stopping path.
+the Dijkstra field, predecessor field, weighted field, and early-stopping path;
+pass `--broad-multisource` to exercise the delta-stepping dispatch.
 The implementation uses reusable fixed neighbor metadata and specialized heap
 strategies internally, but these do not change dtype, tie-breaking, or public
 results.
@@ -2557,19 +2558,19 @@ Important differences and current scope:
   change branch positions and vertex counts, so output is not expected to be
   vertex-for-vertex identical to kimimaro.
 - The C++ core remains dependency-free. `number_of_threads=1` is the default;
-  `0` uses hardware concurrency. One thread budget is shared by the exact
-  distance transform and the internal FP64 compact-foreground Dijkstra backend.
-  Dependent root sweeps and rails remain ordered, while sufficiently large
-  shortest-path wavefronts use deterministic delta stepping. Compact IDs avoid
-  full-volume shortest-path fields. The public Dijkstra functions remain dense
-  and exact FP64.
+  `0` uses hardware concurrency. The thread budget controls the exact distance
+  transform. Compact-foreground root sweeps and rails remain ordered on the
+  optimized heap: threshold benchmarks showed that staged delta stepping nearly
+  doubled runtime for large solid objects. Compact IDs avoid full-volume
+  shortest-path fields. The public Dijkstra functions remain dense and exact
+  FP64, with delta stepping reserved for broad physical multi-source fields.
 
 Correctness tests are under `tests/skeleton/test_teasar.py`. The independent
 Dijkstra benchmark is `development/distance/benchmark_dijkstra.py`; the
 end-to-end synthetic branching-tube benchmark is
 `development/skeleton/benchmark_teasar.py` and can optionally add kimimaro with
 `--kimimaro` when it is installed. Pass `--sequential-backends` to reproduce
-the dense/compact/precision design matrix; extended density, spacing, and PDRF
+the dense/compact FP64 design matrix; extended density, spacing, and PDRF
 regimes are in `development/skeleton/benchmark_teasar_sequential.py`.
 
 ## I/O and Build Dependencies

development/distance/DIJKSTRA_PERFORMANCE_NOTES.md

@@ -543,13 +543,16 @@ gate: matching contracted branch/leaf topology, bidirectional 95th-percentile
 distance at most one normalized voxel, Hausdorff distance at most two voxels,
 and total physical length within 2% of sequential compact FP64.
 
-## Implemented parallel follow-up
+## Initial parallel follow-up (dispatch later narrowed)
 
 Implemented after reviewing the draft above. The implementation scope was
 expanded to the dense public Dijkstra fields as well as the compact TEASAR
 backend. Both use one deterministic FP64 delta-stepping primitive in
 `distance/detail/delta_stepping.hxx`; the optimized heaps remain unchanged.
 
+This section records the initial implementation. The review follow-up below
+supersedes its weighted-field and compact-TEASAR dispatch decisions.
+
 The concrete design closes several gaps in the draft:
 
 - frontier nodes and proposals are processed in fixed global batches independent
@@ -586,3 +589,59 @@ through the measured tiers while sharing the requested thread budget with EDT;
 the delta backend remains available for much larger compact domains. The
 development benchmarks now accept thread matrices so these thresholds can be
 revisited with evidence on other machines.
+
+## Review follow-up: retained delta-stepping scope
+
+The review follow-up on 2026-07-15 narrowed delta stepping to the case for which
+the benchmark shows a consistent benefit: broad multi-source **physical**
+fields. Paths, one/few-source fields, both weighted modes, small problems, and
+all compact TEASAR solves use their specialized heaps. The compact adapter and
+its duplicated sampled-cost-median heuristic were removed completely.
+
+`development/distance/benchmark_dijkstra.py --broad-multisource` adds a full
+source plane and can reproduce the dispatch. Three-repeat final medians were:
+
+| shape / mode | 1 thread | 2 threads | 4 threads | 1-to-4 scaling |
+| --- | ---: | ---: | ---: | ---: |
+| `256 x 256`, physical | 13.57 ms | 5.89 ms | 5.85 ms | 2.32x |
+| `256 x 256`, node-times-physical | 9.80 ms | 9.62 ms | 9.49 ms | 1.03x |
+| `32 x 96 x 96`, physical | 120.55 ms | 60.46 ms | 55.22 ms | 2.18x |
+| `32 x 96 x 96`, node-times-physical | 92.31 ms | 89.31 ms | 92.27 ms | 1.00x |
+
+The weighted rows now show heap noise rather than parallel scaling. Before the
+dispatch was narrowed, their 2/4-thread medians were 149.35/148.64 ms in 3D,
+versus the final 89.31/92.27 ms. Retaining delta stepping for that mode would
+have knowingly shipped a roughly 1.6x regression. The physical 3D row retains
+the useful broad-wavefront speedup (the pre-change medians were 122.79, 60.81,
+and 51.06 ms; the four-thread final minimum was 50.79 ms).
+
+The bucket-index quotient now uses `double`, matching the actual distance and
+delta types on every wheel target. Conversion uses an exclusive, exactly
+representable `2^64` upper bound to avoid undefined out-of-range conversion.
+A decimal-spacing case above the parallel threshold (`0.1, 0.5, 1.25`) asserts
+array-exact distances against the heap and deterministic predecessors between
+two and four threads.
+
+Two proposed changes remain intentionally unimplemented:
+
+- Equal-candidate predecessor handling still uses the existing strict-update
+  behavior. Distances are exact, two/four-thread predecessor fields are
+  deterministic, and refusing equal-distance rewrites prevents cycles in
+  zero-weight components. Harmonizing the small-frontier and merged-proposal
+  tie rules would change an observable predecessor/path choice, so it needs a
+  separately chosen public tie contract rather than being bundled into a
+  portability fix.
+- `std::map` buckets were not replaced speculatively. A flat vector gives O(1)
+  lookup for dense bucket ranges but can allocate excessively when FP64 weights
+  create large sparse index gaps; a ring or segmented design adds rebasing and
+  overflow complexity. Profile bucket lookup/allocation separately on the
+  retained broad physical workload before changing the data structure.
+
+The dispatch threshold (`foreground >= 32768` and at least approximately one
+source per 256 foreground voxels) remains benchmark-derived and intentionally
+conservative. The standalone benchmark records both heap and dispatched cases
+so it can be retuned on materially different wheel hardware.
+
+Final validation used a normal editable build and reported `1145 passed`. The
+optional array-view lifetime fix was also suite-tested in a separate build with
+`-fno-elide-constructors`, specifically removing any dependence on NRVO.

development/distance/benchmark_dijkstra.py

@@ -10,6 +10,8 @@
 --------
 python development/distance/benchmark_dijkstra.py --small --repeats 3
 python development/distance/benchmark_dijkstra.py --repeats 10 --include-geodesic
+python development/distance/benchmark_dijkstra.py --small --threads 1 2 4 \
+    --broad-multisource
 python development/distance/benchmark_dijkstra.py --large --repeats 3 \
     --include-geodesic --json /tmp/dijkstra.json
 """
@@ -64,7 +66,20 @@ def time_call(function, repeats: int, warmup: int) -> list[float]:
     return samples
 
 
-def build_cases(shape: tuple[int, ...], include_geodesic: bool, number_of_threads: int):
+def source_plane(mask: np.ndarray) -> np.ndarray:
+    """Return all foreground coordinates on the first axis-zero plane."""
+    plane_shape = mask.shape[1:]
+    tail = np.indices(plane_shape, dtype=np.int64).reshape(len(plane_shape), -1).T
+    sources = np.column_stack((np.zeros(len(tail), dtype=np.int64), tail))
+    return np.ascontiguousarray(sources[mask[tuple(sources.T)] != 0])
+
+
+def build_cases(
+    shape: tuple[int, ...],
+    include_geodesic: bool,
+    number_of_threads: int,
+    broad_multisource: bool,
+):
     mask = make_mask(shape)
     source = np.full(len(shape), 2, dtype=np.int64)
     target = np.asarray(shape, dtype=np.int64) - 3
@@ -124,6 +139,32 @@ def build_cases(shape: tuple[int, ...], include_geodesic: bool, number_of_thread
             ),
         ),
     ]
+    if broad_multisource:
+        sources = source_plane(mask)
+        cases.extend(
+            [
+                (
+                    "dijkstra/broad-multisource-physical-field",
+                    lambda: bic.distance.dijkstra_distance_field(
+                        mask,
+                        sources,
+                        spacing=spacing,
+                        number_of_threads=number_of_threads,
+                    ),
+                ),
+                (
+                    "dijkstra/broad-multisource-node-times-physical-field",
+                    lambda: bic.distance.dijkstra_distance_field(
+                        mask,
+                        sources,
+                        costs=costs,
+                        cost_mode="node_times_physical",
+                        spacing=spacing,
+                        number_of_threads=number_of_threads,
+                    ),
+                ),
+            ]
+        )
     if include_geodesic:
         cases.append(
             (
@@ -147,6 +188,11 @@ def main() -> int:
     parser.add_argument("--repeats", type=int, default=7)
     parser.add_argument("--warmup", type=int, default=1)
     parser.add_argument("--include-geodesic", action="store_true")
+    parser.add_argument(
+        "--broad-multisource",
+        action="store_true",
+        help="include source-plane fields that exercise parallel delta stepping",
+    )
     parser.add_argument(
         "--threads", type=int, nargs="+", default=[1],
         help="thread counts to benchmark (0 uses hardware concurrency)",
@@ -170,7 +216,10 @@ def main() -> int:
     for shape in shapes:
         for number_of_threads in args.threads:
             mask, cases = build_cases(
-                shape, args.include_geodesic, number_of_threads
+                shape,
+                args.include_geodesic,
+                number_of_threads,
+                args.broad_multisource,
             )
             foreground = int(np.count_nonzero(mask))
             for name, function in cases:

development/skeleton/PERFORMANCE_NOTES.md

@@ -302,7 +302,7 @@ remains FP64. Parallel shortest paths remain a separate follow-up.
 
 Final verification after selecting compact on-the-fly FP64: `1129 passed`.
 
-## Implemented parallel follow-up
+## Initial parallel follow-up (compact dispatch later removed)
 
 The shared thread budget for the distance transform and compact Dijkstra solves
 was implemented and measured on 2026-07-15. The reference matrix now passes
@@ -349,6 +349,9 @@ pool has one useful tracing task and additional workers mostly add process and
 shared-memory overhead. The only measured kimimaro improvement was about 3%
 at four workers on `256^3`.
 
+This section records the initial implementation. The threshold-focused review
+follow-up below supersedes its compact-Dijkstra dispatch decision.
+
 Worker count did not change either implementation's graph size:
 
 | volume | bioimage-cpp vertices / edges | kimimaro vertices / edges |
@@ -447,3 +450,66 @@ and preallocated outputs.
 
 Final verification: `1141 passed` with third-party pytest plugin autoload
 disabled. The normal editable build was restored after profiling.
+
+## Review follow-up: compact dispatch and deferred optimizations
+
+The post-optimization review was resolved on 2026-07-15. The compact
+delta-stepping adapter was removed rather than retained behind its
+`1 << 20`-foreground threshold. It was unreachable in the established sparse
+tube benchmarks and, when exercised by new threshold-focused cases, it caused
+a large regression. Compact root and rail solves now always use their optimized
+heap; `number_of_threads` still controls the exact distance transform.
+
+Reproduce the dispatch matrix with:
+
+```bash
+python development/skeleton/benchmark_teasar_dispatch.py \
+    --threads 1 2 4 --repeats 5 --warmup 1
+```
+
+The benchmark checks array-exact vertices, edges, and radii across worker
+counts. Five-repeat medians before and after removing the adapter were:
+
+| case | workers | before | after | change |
+| --- | ---: | ---: | ---: | ---: |
+| solid cube, 1,061,208 foreground | 1 | 1037.89 ms | 1059.10 ms | +2.0% |
+| solid cube, 1,061,208 foreground | 2 | 1970.56 ms | 1044.70 ms | -47.0% |
+| solid cube, 1,061,208 foreground | 4 | 1872.36 ms | 1028.86 ms | -45.1% |
+| solid cuboid, 1,048,576 foreground | 1 | 1066.11 ms | 1115.94 ms | +4.7% |
+| solid cuboid, 1,048,576 foreground | 2 | 1985.48 ms | 1130.47 ms | -43.1% |
+| solid cuboid, 1,048,576 foreground | 4 | 1943.07 ms | 1093.50 ms | -43.7% |
+| branching tube, 132,619 foreground | 1 | 305.56 ms | 313.56 ms | +2.6% |
+| branching tube, 132,619 foreground | 2 | 228.84 ms | 239.68 ms | +4.7% |
+| branching tube, 132,619 foreground | 4 | 198.84 ms | 199.86 ms | +0.5% |
+
+The small single-thread differences are within the predeclared 5% noise gate;
+the exact-threshold cases show that the removed dispatch nearly doubled
+runtime. The branching tube still scales through the threaded EDT without
+parallelizing its narrow shortest-path wavefronts.
+
+Other review-driven cleanup retained exact behavior: dense TEASAR releases its
+first root field after selecting the final root; dense and compact TEASAR share
+one invalidation-bounds calculation; invalid compact targets are checked before
+workspace mutation; and the failed experimental FP32 backend is no longer
+compiled or exposed by the development selector. Dense FP64 remains an
+independent correctness oracle rather than being folded into a template-heavy
+common rail driver.
+
+The following proposed optimizations were deliberately deferred:
+
+- A pre-sorted target list would replace an O(P x V) scan, but the measured
+  target-selection phase is only 0.8 ms on the `256^3` fixture. Its O(V log V)
+  setup, extra memory, and invalidation bookkeeping cannot pay back materially
+  at present.
+- Touched-node or generation-based compact state reset could avoid a V-byte
+  clear per rail, but all rail-path Dijkstra work is only 26.6 ms in the same
+  profile. Adding a branch/write to every discovery without first isolating
+  reset time risks moving cost into the hotter relaxation loop. Revisit only
+  after profiling reset as its own phase on a many-rail case.
+- A full dense/compact rail-loop abstraction was rejected for now. The indexing,
+  storage, crop handling, and invalidation strategies are intentionally
+  different, and keeping the dense oracle structurally independent helps it
+  catch compact-specific regressions. Only the identical, policy-free bounds
+  calculation was shared.
+- Candidate-target and state-reset changes are therefore benchmark follow-ups,
+  not pending correctness fixes.