refactor(rust): one per-base stats impl, one POD5 batch-read helper - #200
Merged
jayhesselberth merged 1 commit intoAug 23, 2026
Merged
Conversation
signal_stats::compute_signal_stats and inference_pipeline::features::compute_per_base_stats were near-identical copies of the same four statistics, reached by different callers: the first by the Python fast path in features.py, the second by the Rust extraction pipeline. They disagreed on negative map entries -- one cast i64 to usize raw (wrapping, so the base was skipped and left zero), the other clamps to 0 and computes over the truncated span. Reference-anchored maps can carry negative entries, so the same read could get different level features depending on which one ran. The pyfunction is a wrapper now. The POD5 lookup -- parse ids as UUIDs, reads_by_ids, get_signal_bulk -- was written out four times: twice in pod5_io and once each in the training and inference entry points. Four copies is four places to forget cached_reader, and forgetting it is the ~10-80x regression in #176. One helper now, with the GIL-release requirement stated where it can be read. Also warn when SigMapRefiner.algo or sd_params are set. Neither reaches the DP since refine started delegating to escapepod, which builds its own settings and takes an asymmetric dwell penalty in place of leech's short-dwell table. They still feed the reference implementation the parity tests use, so they stay -- but a caller setting them expects refinement to change, and it will not.
jayhesselberth
added a commit
to rnabioco/escapepod-rs
that referenced
this pull request
Aug 24, 2026
`features::span_stats` had the right reduction -- one pass over the covered region with f64 prefix sums, O(1) per span, spans supplied by the caller -- but three of its choices were baked in, so a consumer that disagreed with any of them had to re-implement the whole thing. leech did, twice, and the two copies disagreed on exactly one of those choices (rnabioco/leech#200): the Python fast path skipped a span with a negative start and left zeros, the Rust pipeline computed over the truncated span, for the same read. The three choices are now named fields on a new `SpanConfig`, which replaces the bare `Normalization` parameter: - `SpanStatsOut` grows optional `median` and `range` buffers, built via `SpanStatsOut::new(..).with_median(..).with_range(..)`. Optional because neither can come from the prefix sums. - `SpanFill { Nan, Zero, Value(f32) }` -- `Nan` stays the default and stays the honest answer, but a consumer feeding a neural network needs a sentinel that does not poison the forward pass. - `SpanBounds { Skip, Clamp }` -- `Skip` is the default and the old behaviour; `Clamp` intersects with `[0, len)` and summarises what survives, which is what a reference-anchored map needs. `dwell` under `Clamp` is the clamped length. - `MedianConvention { SelectTotalCmp, SortPartialCmp }` -- both ship rather than one being chosen silently. Measured, they are bit-identical over finite spans (including ulp-separated f32, where the issue expected a ~1e-7 split) and diverge only on a span containing NaN, which SortPartialCmp propagates the way numpy.median does. `SpanConfig::default()` is the old behaviour exactly. A test keeps the pre-`SpanConfig` implementation verbatim as an oracle and asserts the default path is bit-for-bit identical across three normalisations and six fixtures, with and without the optional outputs; flipping either default makes it fail. The Python binding exposes the same knobs keyword-only, all defaulting to the historical behaviour, so callers unpacking three arrays are unaffected. Closes #260
jayhesselberth
added a commit
to rnabioco/escapepod-rs
that referenced
this pull request
Aug 24, 2026
`features::span_stats` had the right reduction -- one pass over the covered region with f64 prefix sums, O(1) per span, spans supplied by the caller -- but three of its choices were baked in, so a consumer that disagreed with any of them had to re-implement the whole thing. leech did, twice, and the two copies disagreed on exactly one of those choices (rnabioco/leech#200): the Python fast path skipped a span with a negative start and left zeros, the Rust pipeline computed over the truncated span, for the same read. The three choices are now named fields on a new `SpanConfig`, which replaces the bare `Normalization` parameter: - `SpanStatsOut` grows optional `median` and `range` buffers, built via `SpanStatsOut::new(..).with_median(..).with_range(..)`. Optional because neither can come from the prefix sums. - `SpanFill { Nan, Zero, Value(f32) }` -- `Nan` stays the default and stays the honest answer, but a consumer feeding a neural network needs a sentinel that does not poison the forward pass. - `SpanBounds { Skip, Clamp }` -- `Skip` is the default and the old behaviour; `Clamp` intersects with `[0, len)` and summarises what survives, which is what a reference-anchored map needs. `dwell` under `Clamp` is the clamped length. - `MedianConvention { SelectTotalCmp, SortPartialCmp }` -- both ship rather than one being chosen silently. Measured, they are bit-identical over finite spans (including ulp-separated f32, where the issue expected a ~1e-7 split) and diverge only on a span containing NaN, which SortPartialCmp propagates the way numpy.median does. `SpanConfig::default()` is the old behaviour exactly. A test keeps the pre-`SpanConfig` implementation verbatim as an oracle and asserts the default path is bit-for-bit identical across three normalisations and six fixtures, with and without the optional outputs; flipping either default makes it fail. The Python binding exposes the same knobs keyword-only, all defaulting to the historical behaviour, so callers unpacking three arrays are unaffected. Closes #260
jayhesselberth
added a commit
to rnabioco/escapepod-rs
that referenced
this pull request
Aug 24, 2026
`features::span_stats` had the right reduction -- one pass over the covered region with f64 prefix sums, O(1) per span, spans supplied by the caller -- but three of its choices were baked in, so a consumer that disagreed with any of them had to re-implement the whole thing. leech did, twice, and the two copies disagreed on exactly one of those choices (rnabioco/leech#200): the Python fast path skipped a span with a negative start and left zeros, the Rust pipeline computed over the truncated span, for the same read. The three choices are now named fields on a new `SpanConfig`, which replaces the bare `Normalization` parameter: - `SpanStatsOut` grows optional `median` and `range` buffers, built via `SpanStatsOut::new(..).with_median(..).with_range(..)`. Optional because neither can come from the prefix sums. - `SpanFill { Nan, Zero, Value(f32) }` -- `Nan` stays the default and stays the honest answer, but a consumer feeding a neural network needs a sentinel that does not poison the forward pass. - `SpanBounds { Skip, Clamp }` -- `Skip` is the default and the old behaviour; `Clamp` intersects with `[0, len)` and summarises what survives, which is what a reference-anchored map needs. `dwell` under `Clamp` is the clamped length. - `MedianConvention { SelectTotalCmp, SortPartialCmp }` -- both ship rather than one being chosen silently. Measured, they are bit-identical over finite spans (including ulp-separated f32, where the issue expected a ~1e-7 split) and diverge only on a span containing NaN, which SortPartialCmp propagates the way numpy.median does. `SpanConfig::default()` is the old behaviour exactly. A test keeps the pre-`SpanConfig` implementation verbatim as an oracle and asserts the default path is bit-for-bit identical across three normalisations and six fixtures, with and without the optional outputs; flipping either default makes it fail. The Python binding exposes the same knobs keyword-only, all defaulting to the historical behaviour, so callers unpacking three arrays are unaffected. Closes #260
jayhesselberth
added a commit
to rnabioco/escapepod-rs
that referenced
this pull request
Aug 24, 2026
`features::span_stats` had the right reduction -- one pass over the covered region with f64 prefix sums, O(1) per span, spans supplied by the caller -- but three of its choices were baked in, so a consumer that disagreed with any of them had to re-implement the whole thing. leech did, twice, and the two copies disagreed on exactly one of those choices (rnabioco/leech#200): the Python fast path skipped a span with a negative start and left zeros, the Rust pipeline computed over the truncated span, for the same read. The three choices are now named fields on a new `SpanConfig`, which replaces the bare `Normalization` parameter: - `SpanStatsOut` grows optional `median` and `range` buffers, built via `SpanStatsOut::new(..).with_median(..).with_range(..)`. Optional because neither can come from the prefix sums. - `SpanFill { Nan, Zero, Value(f32) }` -- `Nan` stays the default and stays the honest answer, but a consumer feeding a neural network needs a sentinel that does not poison the forward pass. - `SpanBounds { Skip, Clamp }` -- `Skip` is the default and the old behaviour; `Clamp` intersects with `[0, len)` and summarises what survives, which is what a reference-anchored map needs. `dwell` under `Clamp` is the clamped length. - `MedianConvention { SelectTotalCmp, SortPartialCmp }` -- both ship rather than one being chosen silently. Measured, they are bit-identical over finite spans (including ulp-separated f32, where the issue expected a ~1e-7 split) and diverge only on a span containing NaN, which SortPartialCmp propagates the way numpy.median does. `SpanConfig::default()` is the old behaviour exactly. A test keeps the pre-`SpanConfig` implementation verbatim as an oracle and asserts the default path is bit-for-bit identical across three normalisations and six fixtures, with and without the optional outputs; flipping either default makes it fail. The Python binding exposes the same knobs keyword-only, all defaulting to the historical behaviour, so callers unpacking three arrays are unaffected. Closes #260
jayhesselberth
added a commit
to rnabioco/escapepod-rs
that referenced
this pull request
Aug 24, 2026
#263) `features::span_stats` was already the right reduction — one pass over the covered region with `f64` prefix sums, O(1) per span, spans supplied by the caller — but three of its choices were baked in, so a consumer that disagreed with any of them had to re-implement the whole thing. leech did, twice, and the two copies disagreed on exactly one of those choices (rnabioco/leech#200): the Python fast path skipped a span with a negative start and left zeros, the Rust pipeline computed over the truncated span. Same read, different features, depending on which path reached it. The payoff here is not line count; it is that the numbers stop depending on which code ran. Precedent: #204, where the rule that decides what a model sees was moved into the crate that owns the reduction rather than re-derived in each caller. ## The three gaps, and how each is closed The bare `norm: Normalization` parameter becomes a `SpanConfig` carrying every knob, so a call site names what it asked for: ```rust pub struct SpanConfig { pub norm: Normalization, pub fill: SpanFill, // Nan (default) | Zero | Value(f32) pub bounds: SpanBounds, // Skip (default) | Clamp pub median: MedianConvention, // SelectTotalCmp (default) | SortPartialCmp } ``` **1. No median or range.** `SpanStatsOut` grows `median: Option<&mut [f32]>` and `range: Option<&mut [f32]>`, built through `SpanStatsOut::new(dwell, mean, sd).with_median(..).with_range(..)` so call sites do not churn on struct literals. They are optional because neither can come from the prefix sums: each needs its own pass over the span, and the median a select or a sort on top of that, which a consumer wanting only dwell/mean/sd should not pay for. `range` is `max - min` after normalisation, propagating `NaN` like `np.ptp`. When either is requested, one gather per resolved span serves both. **2. `NaN` was the only fill.** `SpanFill { Nan, Zero, Value(f32) }`, applied to *every* output array for an unresolved span. `Nan` stays the default and the docstring keeps its argument for why — an unresolved base has no observation, and substituting a value makes it indistinguishable from a real one — with the addition that the argument does not survive contact with a neural network, where a single `NaN` poisons the forward pass. That is why the alternatives exist. **3. Negative span starts were always skipped.** `SpanBounds { Skip, Clamp }`. `Skip` is the default and the old behaviour: an out-of-range coordinate is evidence the map is broken. `Clamp` intersects the span with `[0, signal.len())` and summarises what survives, treating it as unresolved only when the intersection is empty — the right answer under a reference-anchored map, where an entry can legitimately go negative once the aligned region is cropped and the truncated span still carries real signal. Documented explicitly: under `Clamp`, **`dwell` is the clamped length, not the requested width**, because every other output is computed from exactly those samples and pairing a sample count with a mean not taken over that many samples is a contradiction a model could read. ## Why both median conventions, rather than one chosen - `MedianConvention::SelectTotalCmp` (default) is `stats::median_via_select` — `select_nth_unstable` with `total_cmp`, the convention every other median in escapepod-signal already uses. - `MedianConvention::SortPartialCmp` is a full sort with `partial_cmp` plus numpy's own `NaN` check, reproducing `numpy.median` over a `float32` array exactly. Both **average the two middle order statistics** on an even-length span (`(lo + hi) / 2.0`, evaluated in `f32`) and take the middle one on an odd-length span; neither picks one middle and discards the other. That rule is stated in the docs rather than left implied. **A finding that contradicts the issue, reported rather than worked around.** The issue expects the two to split by ~1e-7 on even-length spans of near-equal `float32`. Measured, they do not: over any span of finite values the two are **bit-identical**, including ulp-separated spans. `total_cmp` and `partial_cmp` induce the same order on non-`NaN` values, so both end up averaging the same two elements, and `numpy.median`'s `float32` two-element mean is bit-for-bit `(a + b) / 2.0` in `f32` (checked against numpy 2.5.1 over 400k random pairs). Where they genuinely diverge is a span containing `NaN`: `SelectTotalCmp` sorts it to the high end and returns a finite median from the values below it, while `SortPartialCmp` propagates it. That case is not exotic — a caller padding a window with `NaN` hits it on every padded base — and the propagating answer is the one consistent with `mean`, which is already `NaN` there. So both still ship, and a caller that needs numpy parity can name it; the docs now say precisely where the choice matters instead of gesturing at a difference that is not there. ## The bit-exactness guardrail, and how it is proven `SpanConfig::default()` is the old behaviour exactly. `the_default_config_is_bit_identical_to_the_pre_config_implementation` keeps the **pre-`SpanConfig` function body verbatim** as a test-only `legacy_span_stats` oracle and compares `to_bits()` — not `==`, since most of these are `NaN` — across three normalisations (`None`, `MedianMad{1e-3}`, `MedianMad{1e9}`) and six signal/span fixtures: a 20k pseudo-random read, a short read, a flat read, an empty signal, an empty span list, mixed spans and an edge-case span set covering every flavour of out-of-range. It asserts identity both **without** the optional outputs and **with** them requested, since the entire point of making them optional is that they cannot perturb the prefix-sum path. The guardrail was checked for teeth rather than assumed: flipping the `SpanBounds` default to `Clamp` fails it (and two other tests), and flipping the `SpanFill` default to `Zero` fails it. Both mutations reverted. ## What the tests pin - Bit-exactness of the default path vs. the verbatim old implementation (above). - `sort_partial_cmp_reproduces_numpy_median_and_ptp` — goldens generated in the pixi `python-test` env with **numpy 2.5.1**. The fixture signal is pinned as raw `f32` bit patterns so no decimal literal has to round-trip, and the expected values are `np.median(sig[a:b])` / `np.ptp(sig[a:b])` pinned the same way; the generator snippet is reproduced in a comment above the constants. Spans cover eight ulp-separated values around 1.0 and around 100.7 (shuffled), an odd-length span, a single-element span, a `NaN`-bearing span, and a 33-element random span. - `the_median_conventions_agree_on_every_finite_span` — brute force over ulp-separated spans at five magnitudes and lengths 1–17, plus ordinary-magnitude spans up to length 33, asserting bit equality. This test *is* the documentation for the finding above. - `the_median_conventions_disagree_on_a_nan_span` — the real divergence, with the `mean` already being `NaN` asserted alongside as the argument for the propagating convention. - `median_even_averages_odd_picks_and_a_single_sample_is_itself`, under both conventions. - `range_is_max_minus_min` (including a constant span → `0.0` and a single-element span → `0.0`), and `range_and_median_are_normalised_like_the_mean`. - `the_fill_lands_in_every_output_including_the_optional_ones` — `Nan` / `Zero` / `Value(v)` in all five arrays, resolved spans untouched, and `Zero` proven bit-identical to `Value(0.0)`. - `clamp_summarises_a_truncated_span_where_skip_abstains` — a negative start and a past-the-end span computed correctly under `Clamp` where `Skip` abstains, with `dwell` asserted as the clamped length; a fully in-range span bit-identical under both policies. Plus `clamp_takes_the_fill_for_a_span_that_survives_nothing` and `clamp_widens_the_prefix_sum_region_to_the_clamped_spans`. - `scratch_reuse_is_bit_identical`, extended to cover the new per-span scratch buffer. - Python side: `test_span_statistics_median_range_fill_and_bounds` covers the appended outputs, the fill, both bounds policies, both median conventions by name, the `ValueError` on a bad policy string, the `NaN` divergence cross-checked against `np.median` itself, and the batch path carrying the same knobs. ## Callers The only in-repo caller is the Python binding — `escapepod-classify` has its own unrelated `features` module and does not use `span_stats`, and no bench does either. The binding exposes the same knobs keyword-only, all defaulting to the historical behaviour: `median=True` / `range=True` append a fourth and fifth array to the returned tuple, `fill=<float>` replaces the `NaN` sentinel (`None` keeps `NaN`), and `bounds` / `median_convention` take the policy by name. Callers unpacking three arrays are unaffected. The `.pyi` stub is updated to match. **The `escapepod-classify` charging goldens pass unchanged** — 54/54 including `charging_parity::charging_chain_matches_reference`. Nothing in `tests/fixtures/gen_charging_golden.py`'s output moved. ## Validation All under `srun -p rna`: `cargo fmt --all`; `cargo clippy --workspace --all-targets` clean; `cargo nextest run --workspace` 657/657; `cargo test --doc -p escapepod-signal` 20/20 (the new `span_stats` doctest included); `cargo build -p escapepod-python`; `pixi run -e python-test test-python` 83 passed / 44 skipped; ruff and ty clean. Re-verified after rebasing onto #262: `cargo nextest run -p escapepod-signal -p escapepod-classify -p escapepod-cli` 356/356 on the combined state. --- **Stacked PR.** Based on #262 (#259), not on `main` — review the top commit only. #257 follows on top of this one. (#261, the fourth of the #257–#260 audit, is independent and targets `main` separately.) Closes #260
This was referenced Aug 24, 2026
jayhesselberth
added a commit
that referenced
this pull request
Aug 24, 2026
…sible (#204) leech_core sat at 0.3.0 from v0.3.1 to v0.6.4 -- ten releases, spanning #176, #185, #187, #188, #192, #195, #200 and #202 -- while the Rust changed underneath it. The string is not decoration: uv keys its archive cache on it, so `uv sync` could restore a compiled extension built from any earlier revision that shared the version, over a current build. Caught it doing exactly that: 43 tests failing with pre-#188 behaviour (chunk_signal_kmer_inputs no longer snapping map[0] = 0) against an up-to-date working tree. rust/Cargo.toml is now the single source and tracks leech's version; rust/pyproject.toml takes it through `dynamic = ["version"]` rather than carrying a third copy to keep in sync. leech_core also exports __version__ now, via env!(CARGO_PKG_VERSION), and check_rust() compares it against leech.__version__. The two are separate distributions built from one repository, so a cross-revision pairing does not raise -- it produces different numbers, which is how #176 stayed hidden. check_rust() printed a bare "leech_core" with no version; it now names it and says which half to rebuild, including the case where the stale half is leech's own editable metadata. The release skill bumps both and re-verifies check_rust() afterwards, so this cannot drift again by omission.
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.
Stacked on #199. Last of the audit stack — the duplication that produced the divergences, rather than the divergences themselves. Protected by the parity suite added in #199.
1. Two implementations of the same four statistics
signal_stats::compute_signal_statsandinference_pipeline::features::compute_per_base_statsboth compute per-base mean/median/sd/range. Different callers reach each: the first is the fast path behindfeatures.py::compute_signal_features(the Python backend), the second is the Rust extraction pipeline.They disagreed on negative map entries:
Reference-anchored maps can carry negative entries (
reader.py:101does not clampsig_startwhereprocessing.rs:64does), so the same read could get different level features depending on which backend ran. The pyfunction is a wrapper now.2. Four copies of the POD5 lookup
Parse ids as UUIDs →
reads_by_ids→get_signal_bulkappeared inpod5_io.rstwice and once each intraining.rsandinference.rs. Four copies is four places to forgetcached_reader— and forgetting it is precisely the ~10-80x regression of #176, where a per-callReader::openthrew away the read-id index every batch.One helper (
pod5_cache::read_signals_by_ids, plus a map-returning variant for the three callers that do not need calibration), with the GIL-release requirement documented where someone adding a fifth caller will read it.Net:
pod5_io.rsloses ~60 lines, both pipeline entry points lose ~20 each.3.
SigMapRefiner.algo/sd_paramsare inertNeither has reached the DP since
refinestarted delegating to escapepod, which builds its ownRefineSettingsand uses an asymmetric dwell penalty rather than leech's short-dwell table. They still feed the reference implementation the parity tests exercise, so they stay — but a caller settingalgo="Viterbi"expects refinement to change and it will not, so it now says so.Not done here
The banded-DP scaffolding in
signal_refine.py(refine_signal_mapping,compute_sig_band,convert_to_seq_band,adjust_seq_band,_theil_sen_rescale) is reachable only from tests since refinement moved to escapepod. It reads as dead code but is not:test_rust_refinement_parity.pyuses it as the reference implementation escapepod's DP is checked against. Deleting it would remove a real golden test, so it stays.Beyond that, several of these belong upstream rather than deduplicated locally, and are now filed:
RefineSettingspreset that caused fix(refine): carry #168 across to the Python refiner #192/v0.6.2: Rust dwell/feature window still offset by kmer_context, plus ±1 dwell values #193pod5_cacheitselfbuild_seq_to_sig_map/compute_ref_to_signalspan_statsgaining median/range and a fill policy, which would retire the local implementation this PR consolidatesTests
144 passed across the parity, feature and refinement suites; Rust
fmt+clippy -D warningsclean.