hierarchy: vectorize _get_branch_stats per-label loops (ADR 0013)#242
Merged
Conversation
Replaces 4 per-label Python loops in `Branches._get_branch_stats` (`hierarchical.py`) with vectorized numpy + scipy.ndimage equivalents: - **base_lengths gather**: `label_lengths[unique_labels_int]` with a bounds-mask, replacing `for i, lbl in enumerate(unique_labels): ...`. - **Tip-radius adjustment** (lone + regular tip loops): `np.searchsorted(unique_labels, tip_labels)` + `np.add.at`. Multiple tips per label accumulate via `np.add.at`'s sequential semantics (bare indexed `+=` would only apply the last addend for duplicates). - **Median thickness**: `scipy.ndimage.median(thicknesses, labels, index=unique_labels)`. Bit-identical to per-group `np.median` on the yeast 2D + 3D fixtures and on a 5000-element synthetic case. Defensive `np.bincount`-driven NaN-fill for absent labels (no-op on the test fixtures since `unique_labels = np.unique(L[L>0])` guarantees every label has at least one matching voxel). - **Tortuosity**: stable `np.argsort(tip_labels)` + boundary scan to grab the first two tips per qualifying label, then a single vectorized distance + division. 2D/3D split falls out of `tip_coords.shape` directly (the now-removed `no_z` local). The float32 cast point shifts: legacy casts after every `+=` (one cast per tip), rewrite accumulates in float64 and casts once at the end. **More numerically accurate** in the IEEE sense (single rounding vs. accumulated round-after-every-step), but drifts by 1 float32 ULP on labels with multiple tips. Observed on yeast 3D fixture, label 2 at t=0 (the longest branch): `branch_length_raw` 16.6453686 → 16.6453667, a 1.9e-6 absolute / 2e-7 relative diff. Tortuosity inherits the same drift in the same row. Bar: `rtol=1e-5, atol=1e-5` (see ADR 0013); 2D fixture is bit-identical (no multi-tip labels); `branch_thickness` is bit-identical on both fixtures. New imports: `from scipy import ndimage as ndi_cpu` (scipy already a transitive dependency). Tests: 3 new regression tests pin the equivalence bar: `test_get_branch_stats_2d_post_rewrite_bit_identical`, `test_get_branch_stats_3d_post_rewrite_approx_equivalent`, and `test_get_branch_stats_run_is_deterministic` (in-band determinism across two independent runs on fresh fixtures). ADR 0013 documents the trade-off and the 4 considered options (rejected: per-tip Python loop with intermediate cast, pure-float32 add.at, single+multi-tip detection split, defer entirely). Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
…hetic multi-tip test Frangi runs at session-start (`frangi_3d_path` → `_run_filter_to_disk`) and is SIMD-sensitive across macOS / Linux / Windows. Different Frangi → different segmentation → different branch counts on the 3D yeast fixture (Linux/Windows produce 4 branches, macOS produces 6), so a hardcoded 3D baseline can't be cross-platform stable. Replace `test_get_branch_stats_3d_post_rewrite_approx_equivalent` with `test_get_branch_stats_multi_tip_label_synthetic_post_rewrite`: constructs a 9-voxel "+" skeleton (4 tips on one label) directly, bypassing all SIMD-sensitive upstream stages. Asserts the post-rewrite single-end-cast value equals the formula-derived expected value, and that the legacy per-tip-cast pattern stays within `rtol=1e-5, atol=1e-5` of the new value. Platform-stable by construction (all inputs are hand-set integers + irrational-at- float32 distances chosen to make the cast point matter). Updates ADR 0013 to reflect the test-bar shift and the cross- platform Frangi limitation. The 2D bit-identity + determinism + synthetic tests together cover the equivalence claim without depending on SIMD-stable Frangi output. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
2 tasks
aelefebv
added a commit
that referenced
this pull request
May 12, 2026
…ncy (#243) Two safe-and-free cleanups in `hierarchical.py`: - **`_get_ref_coords` redundant gather** — `vals_a` and `vals_b` were computed by the *same* `idxmin[branch_labels_clipped]` expression (same input array, same lookup). Compute once, share the NaN mask. ~6 lines deleted; fewer allocations on the hot motility path (called per-frame inside `_get_motility_stats`). - **`_save_adjacency_maps` voxel→node Python loop** — replaces the nested `for voxel_idx, nodes in enumerate(...): for n in nodes: edges_vn.append((voxel_idx, int(n)))` with `np.repeat` + `np.concatenate`. Only fires when `not skip_nodes` (production default in the napari pipeline). Removes per-edge Python attribute + tuple-allocation overhead — a 1MB-edge frame goes from millions of Python ops to two C-level numpy calls. Bit-identical: verified locally on yeast 2D + 3D fixtures via SHA over all 6 hierarchy outputs (`features_*` + `adjacency_maps`) — all 12 hashes match the post-PR-#242 baseline exactly. Co-authored-by: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
Merged
1 task
aelefebv
added a commit
that referenced
this pull request
May 12, 2026
…ation audit close-out (#244) - now.md: new "Recently shipped" entry covering PRs #241/#242/#243 (vectorized per-label group construction, vectorized _get_branch_stats loops + ADR 0013, _get_ref_coords + vox→node cleanups). Watch list updated to reflect hierarchical's pinning status and the deferred wins (regionprops solidity ~17%, flow_interpolation _get_vector_weights ~28%) that the next perf audit could target. - feature-extraction.md: new Performance section documenting the 4 shipped optimizations + the 3 deferred items, plus a Benchmarks subsection pointing at tests/test_hierarchical_perf.py. Co-authored-by: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Replaces the 4 per-label Python loops in `Branches._get_branch_stats`
(`hierarchical.py`) with vectorized numpy + scipy.ndimage equivalents:
Equivalence bar
Approx-equivalent (`rtol=1e-5, atol=1e-5`) per ADR 0013. Float32 cast point shifts: the legacy loop casts after every `+=` (one cast per tip), the rewrite accumulates in float64 and casts once at the end. More numerically accurate in the IEEE sense (single rounding vs. accumulated round-after-every-step), but drifts by 1 float32 ULP on multi-tip labels.
Observed drift on yeast 3D fixture, label 2 at t=0 (the longest branch):
Bit-identical for:
Why
cProfile on the yeast 3D fixture (`Hierarchy.run()`, 0.42 s total) flagged `_get_branch_stats` at 62 ms (15%) — second only to motility (which is in `flow_interpolation.py`, out of scope). The 4 per-label loops are O(N · L) (full-volume comparison per unique label); vectorized variants are O(N log N) or O(N + L). On the test fixtures the wall-clock win is modest (small label count); on production data with hundreds-to-thousands of branches per frame, the savings should be visible.
Test plan
🤖 Generated with Claude Code