Skip to content

v0.15.0

Choose a tag to compare

@github-actions github-actions released this 24 Aug 00:56
· 35 commits to main since this release
bdb1d92

Build / Tooling

  • The POD5 compat job stops rebuilding its dependency graph every run.
    0.14.0 moved it off --release onto an optimised-but-not-LTO profile, which
    did not work: Swatinem/rust-cache derives its key from Cargo.lock and the
    toolchain, not from the cargo profile, and GitHub refuses to overwrite an
    existing cache key. The key therefore still held the old release artefacts,
    so every run restored artefacts it could not use, rebuilt everything, and
    then declined to save because the key already existed. Measured on main:
    363 s of a 400 s job was cargo build, the compat test itself was 1 s,
    and the cache post-step wrote nothing. Warm cost went 316 s → 400 s — the
    change cost more than it saved.

    Two fixes, both needed. The cache key now carries an explicit suffix that is
    bumped whenever the profile changes, so a profile switch can actually be
    saved. And the job builds the dev profile rather than an optimised one,
    because the suite round-trips a five-read fixture and the binary's
    throughput is irrelevant to it; third-party crates still compile at
    opt-level = 2 through [profile.dev.package."*"], so only escapepod's own
    crates drop to -O0, and those have to rebuild on any source change anyway.
    The now-unused ci-bin profile is removed.

Added

  • A cache of open, indexed readers, so the read-id index is built once per
    file instead of once per reader (#258).
    Reader caches its index in a
    OnceLock on the instance, so a consumer that opens a reader per batch
    throws the index away and rebuilds it on the next batch. That is not a small
    constant: on a 145 GB POD5 on a network filesystem it was minutes of
    uninterruptible sleep in folio_wait_bit_common per batch at ~0.6% of one
    core — the 10–80x data-preparation regression in rnabioco/leech#176. #251
    fixed the other half of it (the scan variants are gone and lookups index
    unconditionally), but "one reader per file per process" was left to every
    consumer, and each consumer that did not write it silently got the slow path.
    leech wrote it in Rust, and then wrote the same idea again, independently, in
    Python.

    Both shapes ship, because they answer different questions. cached_reader()
    is the process-global convenience, and it is the one that makes consumers
    actually stop hand-rolling this. ReaderCache is the owned type underneath
    it, for a library that needs the lifetime bounded or a process where one
    stage must not share readers with another; global_reader_cache() reaches
    the global's len() / clear().

    The value is in the ordering and the failure semantics rather than in the
    static, so those are the parts worth stating:

    • The file is opened outside the lock, which guards only the map. A slow
      open on one path never blocks a lookup on another, and the lock is never
      held across I/O, so this cannot deadlock. Two threads racing the same path
      cost one redundant open and both get the winner's Arc — publication goes
      through entry, not insert, so a race can never leave two live readers
      (and two indexes) for one file.
    • The index is warmed before the entry is published, so N workers hitting
      their first batch together find it built instead of piling up inside one
      lazy init. The warm-up respects autoindex_max(): above that read count it
      is skipped, because warming is a guess that random access is coming and a
      huge file that is only iterated should not pay for an index nobody asked
      for. Skipping only defers the build to the first lookup that demands one,
      and because the reader is now shared that build still happens once per file
      rather than once per batch — the cache keeps its whole value above the
      threshold, it just stops guessing.
    • A failed index build is logged, not propagated. Reader::open failing
      is an error, because there is no reader to hand back. An un-indexable
      POD5 is still a perfectly good reader for iteration, metadata, and signal
      access, and failing an open for a caller that may never do a lookup is
      worse than the slowdown; a caller that does demand a lookup sees the same
      error then, from the call that needs it. (One correction to the issue's
      framing: after #251 such a file is not "readable, just slowly" — the error
      surfaces from reads_by_ids rather than degrading to a scan. The reader
      stays usable; lookups by read id do not.)

    Keys are canonicalized, falling back to the path as given if that fails, so
    reads.pod5, ./reads.pod5, and a symlink to it are one entry rather than
    three readers with three indexes. The reader is opened on the canonical path
    too, so .p5s sidecar resolution does not depend on which spelling happened
    to arrive first. What stays resident is the index and not the file —
    ~24 bytes/read, so a few tens of MB even for a multi-million-read POD5 — and
    entries are never evicted, with clear() as the escape hatch for a process
    that walks an unbounded set of files.

  • Reader::read_index_if_built() — the non-committing half of
    read_index(): it never loads a sidecar and never scans, so it is the only
    way to ask whether a reader is warm without making it warm. Without it the
    warm-before-publish ordering above is unobservable, and a test that "checked"
    it by calling read_index() would only be asserting its own side effect.

  • escapepod_signal::mapping: the two Oxford Nanopore coordinate
    conventions that produce a resquiggle's input.
    refine_signal_map has
    always taken a sequence→signal map; nothing in the workspace produced
    one. So every consumer wrote its own eight lines off the mv/ns/ts tags
    and its own CIGAR walk — three copies in this repo alone (the charging
    classifier's anchoring, the resquiggle command, a test helper), plus the
    ones downstream. Each is a shifted map away from answering a different
    question than the caller thinks, with no error to show for it, which is the
    same argument that moved the k-mer level primitives here.

    • seq_to_signal_from_moves(moves, stride, trim_offset, num_samples)
      Remora's query_to_signal = np.nonzero(mv)[0] * stride, returned in
      trimmed-signal coordinates with num_samples - trim_offset as the
      closing boundary, because that is the frame the move table is in and the
      frame refine_signal_map is handed. A caller indexing the untrimmed POD5
      array adds trim_offset back; the charging anchoring now does that
      explicitly instead of folding + ts into the map's construction, where
      the frame was invisible.
    • ref_to_signal(query_to_signal, cigar) — reference→signal by the Remora
      knot convention: trailing non-match ops stripped, knots at the start and
      end - 1 of each match block (not end, which stretches every gap by a
      position), exact 1:1 integer lookup inside a block, and linear
      interpolation only across indel gaps.

    The CIGAR arrives as a local CigarOp { kind, len } rather than the
    (op, len) integer pair the convention is usually written with: the crate
    takes no alignment-library dependency for this, and a bare pair of integers
    is exactly what a caller transposes without the compiler noticing.

    ref_to_signal is integer arithmetic throughout except the one ratio each
    gap position needs — deliberately not the ref → float query → float signal
    chain that a pair of np.interp calls performs. Both interpolations there
    evaluate slope * (x - x0) + y0 with a pre-rounded slope, and the result is
    floored, so a one-ulp difference in the intermediate query coordinate
    becomes a one-sample difference in the answer: with the map [0, 7, 8] and
    a CIGAR of 1M 6D 1M, the float chain puts reference position 5 at sample 4
    instead of 5. It is rare — a 200 000-case sweep of realistic random CIGARs
    found no difference at all, and it takes a long deletion spanned by short
    dwells — which is precisely what makes it expensive to find once two
    consumers have each written their own version. It is pinned by a test here
    rather than rediscovered downstream.

Changed

  • features::span_stats gains a median, a range, a fill policy and an
    out-of-range policy, and takes a SpanConfig instead of a bare
    Normalization (#260).
    The reduction was already the right one — 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, and a consumer that
    disagreed with any of them could not use the function at all. leech therefore
    carried its own copy, then a second copy, and the two 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 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 to the crate that owns the reduction rather than
    re-derived in each caller.

    The three gaps, all now named fields on SpanConfig rather than assumptions.
    SpanStatsOut grows optional median and range buffers, built through
    SpanStatsOut::new(..).with_median(..).with_range(..) — 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), so a caller wanting only
    dwell/mean/sd does not pay for them. SpanFill { Nan, Zero, Value(f32) }
    chooses what an unresolved span gets: Nan stays the default and stays the
    honest answer — an unresolved base has no observation, and a substituted
    value is indistinguishable from a real one — but that argument does not
    survive contact with a neural network, where one NaN poisons the forward
    pass, so the alternatives exist for a caller feeding these arrays to a model.
    SpanBounds { Skip, Clamp } chooses what happens to a span hanging off the
    end: Skip (the default, and the old behaviour) treats an out-of-range
    coordinate as evidence the map is broken; Clamp intersects with
    [0, len) and summarises what survives, which is what a reference-anchored
    map needs once the aligned region is cropped and entries can legitimately go
    negative while the truncated span still carries real signal. Under Clamp,
    dwell is the clamped length, not the requested width — every other
    output is computed from exactly those samples, and pairing a sample count
    with a mean not taken over that many samples would be a contradiction a model
    could read.

    Both median conventions ship, rather than one being chosen silently.
    MedianConvention::SelectTotalCmp (the default) is
    stats::median_via_select, i.e. 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, for a consumer that cross-checks against a Python reference. Both
    average the two middle order statistics on an even-length span and take the
    middle one on an odd-length span; neither picks one middle and discards the
    other. Measured, the two are bit-identical over any span of finite values,
    including the even-length ulp-separated f32 spans where a ~1e-7 split was
    expected — total_cmp and partial_cmp induce the same order on non-NaN
    values, and numpy.median's float32 two-element mean is bit-for-bit
    (a + b) / 2.0. 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, SortPartialCmp propagates it. That 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. Both behaviours are pinned by tests, the numpy arm against goldens
    generated from numpy 2.5.1.

    The API break is behaviour-preserving by construction and proven so.
    SpanConfig::default() is the old behaviour exactly (Nan fill, Skip
    bounds, no normalisation), so span_stats(sig, spans, SpanConfig::new(norm), ..) is the old call. A test keeps the pre-SpanConfig implementation
    verbatim as an oracle and asserts the default path is bit-for-bit identical
    to it across three normalisations and six signal/span fixtures — including
    with the optional outputs requested, since the point of making them optional
    is that they cannot perturb the prefix-sum path. The guardrail was checked
    for teeth by flipping each default in turn and confirming it fails.

    The Python binding exposes the same knobs, keyword-only and 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, and bounds / median_convention take the policy by name. Callers
    that unpack three arrays are unaffected.

Fixed

  • One named refinement preset, with a per-read dwell target (#257). The
    settings block for refining a basecaller move table existed twice — once in
    escapepod's own Python binding (py_refine_signal_map) and once in a
    downstream Rust consumer — each carrying a comment asserting that it matched
    the other. The binding's docstring went further and promised the two paths
    matched "bit-for-bit". They did not: dwell_target had drifted, a fixed
    4.0 in the binding against the 0.0 sentinel that asks escapepod to
    resolve the target from the read's own move-table median dwell.

    That one field is not cosmetic. The dwell penalty is asymmetric — quadratic
    below target, logarithmic above — so a target set too low does not merely
    weaken the prior, it actively drags boundaries toward dwells the pore never
    produced. RNA004 at 130 bases/s and 4 kHz sits near 31 samples/base, so a
    target of 4.0 treated every base as roughly 8x too long. Measured across
    the two backends on the same reads with the same flags: max |signal delta|
    3.44 in normalized units, every dwell different, max |feature delta| 3.57
    .
    The two paths refined the same data to different boundaries for four
    releases, and the comment saying they agreed was there the whole time.

    RefineSettings::move_table_refinement(half_bandwidth, n_iters, seed) is now
    that configuration, as a value rather than a convention: fixed banding, a
    least-squares rough rescale over the 0.05–0.95 quantiles clipped 10 bases
    with use_base_center, a Theil-Sen inter-iteration rescale over at most 200
    points, level normalization off, and the asymmetric dwell penalty at weight
    0.5 with the per-read target. The sentinel gets a name —
    RefineAlgo::PER_READ_DWELL_TARGET — because 0.0 at a call site does not
    say what it means, and this is the field where that cost something.
    RefineAlgo's shape is unchanged, so nothing downstream has to move.

    Behaviour change for escapepod.refine_signal_map. dwell_target and
    dwell_weight become Optional[float], default None, meaning "use the
    preset"; passing a number still overrides it. The old default of 4.0 is
    gone rather than preserved. It is simply wrong for RNA004, it silently
    corrupted a production corpus, and a caller who wants it back can pass it
    explicitly — which is a better trade than making every future caller inherit
    a known-wrong number for bug-compatibility.

    The docstring stops promising bit-for-bit parity, since that promise is not
    enforceable from inside a docstring and was false when written; it now names
    the preset both paths construct, which is checkable. It also settles what
    (scale, shift, drift) are for. The return tuple is unchanged, and it had
    instructed callers to apply the rescale as (signal[i] - shift - drift*i) / scale while the downstream Rust path deliberately discarded those same
    values. Both readings were defensible because escapepod never said which it
    intended. It now does: the values are returned for inspection, applying
    them is the caller's decision, and the failure mode is documented — a
    per-read affine fit estimated over a near-constant stretch of signal (a 3'
    adapter, a homopolymer) is weakly identified, with observed scales ranging
    from 15 to 1084 and frequently negative.

    Three tests pin this. The preset's fields are asserted one by one, including
    the rescale filter constants and the quantile grid, so a future edit to any
    default cannot quietly redefine the preset. Refining an RNA004-like synthetic
    read under the preset must reproduce refining it under an explicitly named
    target equal to the input map's median dwell — "per-read" stated as something
    observable rather than as prose. And the same read must refine differently
    under a fixed 4.0; restoring the old default fails that test, which was
    confirmed by restoring it.