perf: choose pivots before sorting, and fuse phase 3 into phase 1 - #9
Closed
BenjaminDEMAILLE wants to merge 41 commits into
Closed
perf: choose pivots before sorting, and fuse phase 3 into phase 1#9BenjaminDEMAILLE wants to merge 41 commits into
BenjaminDEMAILLE wants to merge 41 commits into
Conversation
Two prerequisites for the performance work that follows, neither of which changes behaviour. The crate had no test covering the LCP array. That is the riskiest possible gap for this algorithm: the public entry points discard the array, but the *next* merge level consumes it in the three-case decision, so a single wrong LCP entry silently reorders suffixes one level up and the SA comes out subtly wrong. Add four tests that check `lcp[0] == 0` and `lcp[i] == lcp(text[sa[i-1]..], text[sa[i]..])` against a naive oracle, over fixtures, random texts across four alphabet sizes, long runs and periodic text, and finite `max_context`. `bench/README.md` claims the published numbers were taken with fat LTO and one codegen unit, supplied by a parent workspace. That workspace is not in this repo, so every build made from it since the crate went standalone has used `lto = false, codegen-units = 16`. Pin the profile here. In a library crate `[profile.release]` applies only when this crate is the workspace root, so it affects this repo's own tests, examples and benches and is invisible to downstream consumers. Measured on Apple M4 Max, 12 threads, chr21 FASTA (47.5 MB): 27.8 s -> 24.0 s wall. Neutral on N-free DNA. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The CaPS-SA merge kernel stays the general path, but it is not the right
algorithm for the most common request: the standard lexicographic suffix
array of a byte text, unsegmented, with no context bound. Add
`src/radix.rs` and route that case through it.
Profiling separated two distinct costs, and the merge kernel pays both:
* Step count. `n log n` merge steps, most resolved by a symbol
comparison at a random text address. 80 MB of N-free DNA is ~2.1e9
steps at ~13 ns each.
* Scan length. Every leaf merge starts at `m = 0`, so two suffixes
sharing a long prefix cost a scan proportional to that prefix. Genome
FASTA carries megabyte-scale runs of `N` (period-61 once line wrapping
is included), where one comparison scans millions of bytes. That drives
the cost per merge step from 13 ns to 222 ns, a 16x penalty which is
entirely scan time. This is what made real chr21 20x slower than
N-free DNA of comparable size.
The new path removes both. It sorts by a packed fixed-depth key, then
resolves the remainder by prefix doubling on ranks. The packing picks the
narrowest field width in {1,2,4,8} bits that holds the alphabet, so DNA
over {0,1,2,3} resolves 32 symbols per key rather than the 8 a raw byte
key gives. After the seed no comparison reads the text again, so a
megabyte run of `N` costs exactly what random DNA costs.
The seed sorts by `(key, min(n - p, k))`. The second component is
required, not cosmetic: zero-padding makes a short suffix share a key
with any suffix continuing in zeros, and `0` is a real symbol in every
DNA encoding. Ordering by visible length puts the proper prefix first,
which is the crate's shorter-is-smaller convention. Without it `[0, 0]`
leaves two positions permanently tied and doubling cannot terminate.
Guards are soundness conditions, not heuristics, and all three default
to declining:
* `max_context` must be unbounded; a finite bound makes the merge's
comparator fall through to `boundary_order`, which compares lengths,
so it is not lexicographic.
* `LimitProvider::plain_lex_len` must report the full text. New method,
defaulting to `None`, overridden only by `PlainText`. An impl that
delegates `lim_at` to `PlainText` but overrides `boundary_order` for a
different convention (STAR's spacer-as-largest) inherits `None` and
stays on the merge kernel without changing a line.
* `S` must be exactly `u8`. Packing wider symbols into an order-
preserving key is endianness-dependent: for `u16` on a little-endian
host `0x0100 > 0x0001` as values but their byte views compare the
other way. The rest of the crate avoids this only because
`LcpDispatch` resolves equality over bytes and recovers ordering
through `S: Ord`.
Tests: exhaustive over every binary text to length 10 and every ternary
text to length 6, random texts across seven alphabet widths, texts where
a real `0` collides with padding, long runs, periodic text, and the
wrapped-FASTA `N`-block shape.
Measured on Apple M4 Max, 12 threads, against the previous kernel, with
byte-identical suffix arrays on both real inputs:
chr21 fwd+revcomp, N-free, 80 MB 6.08 s -> 1.14 s CPU 28.1 s -> 5.0 s
chr21 FASTA, 47.5 MB 27.8 s -> 1.16 s CPU 283.5 s -> 5.2 s
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Checking a new construction path against the old one only shows the two agree. Checking adjacent suffixes directly is O(n · lcp), which is unusable on exactly the repetitive inputs that need checking most: on a chr21 FASTA a single adjacent pair can share megabytes. Use the fixpoint characterisation instead. With `rank` the inverse of `sa` and `f(p) = (text[p], rank[p + 1])`, taking `rank[n]` as less than every real rank, a permutation of `0..n` is the suffix array of `text` if and only if `f` is strictly increasing along it. That is one pass to invert plus one pass to compare, independent of any construction algorithm and independent of LCP length. Exposed as `caps_sa::verify_sa` and wired to a `--verify` flag on the bench CLI, off by default so it never contaminates a timing run. Full-scale results on Apple M4 Max, 12 threads: chr21 fwd+revcomp, N-free, 80,177,238 entries verify OK in 1.14 s chr21 FASTA, 47,488,540 entries verify OK in 0.57 s Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Phase timings (`CAPS_SA_PROFILE=1`) showed the two rank-scatter passes
were fully sequential and had become the largest single cost after the
seed sort: 0.34 s of a 1.17 s build on 80 MB of DNA.
Both scatter through a permutation -- the target index is `sa[i]`, not
`i` -- so the writes are not expressible as disjoint sub-slices and
`split_at_mut` does not apply. They are nonetheless disjoint: `sa` is a
permutation and the ranges being processed partition its index space, so
every slot is written exactly once. Introduce a small `Scatter` wrapper
that encodes precisely that contract in its `unsafe fn set`, and drive
both passes with rayon.
Grouping now has each index decide for itself whether it starts a group;
the index that does owns the group, walks it to find the end, and writes
its members' ranks. Exactly one owner per group, and `collect` on an
indexed parallel iterator preserves order, so the group list still comes
out sorted, which is what `split_disjoint` relies on.
Also materialise the successor ranks before sorting each group.
`sort_unstable_by_key` re-evaluates its key function O(len log len)
times and every evaluation was a random probe into `rank`; paying once
per element makes the sort's memory traffic sequential.
Apple M4 Max, 12 threads, suffix arrays byte-identical to the previous
kernel and independently `--verify`-checked:
chr21 fwd+revcomp, N-free, 80 MB grouping 0.341 s -> 0.075 s
total 1.17 s -> 0.89 s
chr21 FASTA, 47.5 MB grouping 0.172 s -> 0.037 s
total 1.10 s -> 1.04 s
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The repository had no Rust CI at all: the only workflow deployed the docs site, while 73 tests sat in the tree with nothing running them. That is not a safe baseline for changing the sorting kernel. Covers both architectures that matter here, since the LCP kernel and the pooled external-memory bucket path are the parts that diverge per platform: macOS is aarch64/NEON, Ubuntu is x86_64/AVX2. Runs the tests in debug as well as release, because debug is what exercises the `debug_assert`s guarding the unchecked scatter in `radix.rs` and the buffer-length invariants in the merge kernel. Adds fmt, clippy with warnings denied, a check against the declared 1.89 MSRV, and rustdoc with broken links denied. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
`bench/run.sh` compares against upstream C++ and needs both binaries
prebuilt. Add a self-contained script for the case that prompted this
work: it fetches hg38 chr21 and prepares *both* inputs, which is the
distinction that made the original slowdown report hard to interpret.
chr21.0123 forward ++ revcomp, one byte per base, codes 0..=3,
ambiguous bases dropped. ~80 MB, alphabet 4, no long runs.
The input libsais is normally benchmarked on.
chr21.fa the raw FASTA, still carrying its ~6.6 Mb of `N`. Wrapped
at 60 columns, so the `N` blocks are a period-61 repeat
rather than a plain run.
Benchmarking one implementation on the first and another on the second
compares two different problems. The second is the realistic input and
the one that used to be pathological.
Builds with `-C target-cpu=native` and runs each case through
`--verify`, so the harness reports correctness alongside timing.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Adds to README: the in-memory fast path with measured numbers, a "Choosing a path" table giving the three soundness conditions and why each one exists, and `verify_sa` usage. Adds to bench/README: the chr21 section, including the per-merge-step arithmetic that separates the two inputs (13 ns/step on N-free DNA against 222 ns/step on the FASTA, same kernel and same machine), and the phase breakdown of the fast path. Corrects two claims that were misleading: The build paragraph attributed `lto = "fat"` and `codegen-units = 1` to a parent workspace. No such workspace is in the repository, so from the commit that made the crate standalone until the profile was added, every build made from this repo used `lto = false, codegen-units = 16`. Numbers taken in that window are not comparable with numbers taken now. The 97.54% `lcp_u8_avx2` profile was read as "LCP scanning is expensive", which motivated widening the scan through AVX2, AVX-512 and the hybrid. The LCP kernel is also where the two random text loads happen, so on short-LCP input those samples are load stalls and a wider vector cannot help. The existing AVX-512 ablation already showed this: the 64-byte-only variant was 16% slower on rand100m and only the long-LCP human slice gained. Both readings are now stated. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Public documentation linked to `FilteredSource`, which is private, so `cargo doc` fails under `RUSTDOCFLAGS=-D warnings`. Pre-existing, but it blocks the rustdoc job added in the previous commit. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The seed was a `par_sort_unstable` over a materialised `Vec<(u64, u32, I)>`, which was both the largest remaining cost and the peak-memory driver. Three things change. The source is never materialised: keys are recomputed from `text` in the histogram pass and again in the scatter, which trades a random read of an n-element key array for a sequential read of the text. Peak memory drops to the two destination buffers, 12 bytes per position at `I = u32` against the 16 a `(u64, u32, I)` record costs. And the top-level partition becomes a counting pass, which parallelises evenly, where a parallel comparison sort's first partitioning steps are close to serial. 11 bits (2048 buckets) keeps the write-combining state near 512 KB and inside a core's private cache. 16 bits would need 16 MB of open write lines and thrash the TLB instead. The visible-length tie-break no longer needs storing. Only the last `k - 1` positions can have a visible length below `k`, so it is a function of the position alone and is applied in the per-bucket sort and in the group scan. Apple M4 Max, 12 threads, chr21 fwd+revcomp 80 MB, output unchanged and `--verify` clean: peak RSS 2.83 GB -> 2.21 GB (-22%) seed 0.341 s -> 0.289 s total 0.892 s -> 0.844 s Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Closes two of the three gaps listed as limitations when the fast path landed. **Subsets.** Doubling cannot be restricted to a subset directly, because a round compares `rank[p + d]` and that successor is generally not in the subset, so ranks must be defined for every position in the text. Building the whole array and filtering it sidesteps that, and the filter is one O(n) pass since the full array is already ordered. Worth it when the subset is a real fraction of the text, which is what this API exists for: STAR-style indexing keeps every ACGT position and drops only spacers. Below one eighth of the text it declines, because O(n) to build and discard would dwarf the O(m log m) the merge kernel needs. That ratio is a performance heuristic; the guards it sits behind remain correctness conditions. It also declines on duplicate or out-of-range positions. The output is a permutation of the input *multiset*, which a membership filter cannot reproduce. **In-memory sample sort.** `build_in_memory_sample_sort` exists to sort in RAM, so where the doubling path applies it is strictly better: same output, no bucket machinery, and none of the scan cost on repeat-heavy text. `build_ext_mem` deliberately does *not* get this. Its purpose is to bound peak memory, and routing it through an in-memory algorithm would defeat exactly that. It stays on the merge kernel, and the remaining limitation is now stated as a deliberate choice rather than an omission. Apple M4 Max, 12 threads, chr21 fwd+revcomp 80 MB. `--in-mem-ss` output verified identical to the in-memory path's: --in-mem-ss 3.41 s -> 1.18 s wall, 34.5 s -> 7.2 s CPU Tests: duplicate positions keep their multiplicity, a tiny subset of a large text still matches brute force through the merge kernel, and out-of-range positions still panic rather than silently returning a wrong answer. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Records the counting-sort seed (peak RSS 2.83 GB -> 2.21 GB, total 0.89 s -> 0.84 s on the 80 MB input), adds the `--in-mem-ss` row, and replaces the "everything else uses the merge kernel" line with what is now actually true: subsets and in-memory sample sort are covered, and `build_ext_mem` stays on the merge kernel deliberately, because bounding peak memory is the whole point of that path. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
`build_ext_mem` was left on the merge kernel because bounding peak
memory is its purpose and prefix doubling needs a rank for every
position in the text. That left it paying the full scan cost on
repeat-heavy input. Profiling put 94% of an ext-mem run on chr21 FASTA
in phase 1:
chr21.0123 (N-free, 80 MB) phase1 1.15 s of 3.49 s total
chr21.fa (47.5 MB) phase1 22.85 s of 24.19 s total
Normalised that is 0.014 s/MB against 0.48 s/MB, the same 34x scan
penalty the in-memory path had, and for the same reason: two suffixes
inside a long repeat agree for as far as it continues, so one comparison
scans megabytes.
Fix it in the comparator rather than the algorithm, which keeps the
memory bound intact. If `text[s..e)` has period `q` and two suffixes
start at `a < b` inside it with `(b - a) % q == 0`, they agree until the
later one reaches `e`, so `lcp(a, b) >= e - b` is known in O(1) from the
run's bounds with no scanning. When the phase does not match, the two
must differ within `q` symbols and the ordinary scan is already short.
The scan is additionally bounded so it stops at a run's start rather
than traversing it.
Detecting only single-symbol runs would have missed the case that
actually occurs: in wrapped FASTA an `N` block is 60 `N`s then a
newline, which is period 61, not period 1. Periods up to 64 are
considered.
Detection is two-stage so texts without runs pay almost nothing. A
sampling pass looks for any periodic window and collects the periods
that occur; the full scan runs only for those. On N-free DNA the sample
finds nothing, the table is empty, and every query short-circuits on a
slice-empty check. The table itself is a few dozen entries, so the
memory bound is untouched.
`Cmp` bundles the SIMD dispatch with the run table and replaces the bare
`LcpDispatch` threaded through the merge kernel, so phase 1, the phase-3
pivot searches and the phase-4 cascade all benefit. It stays `Copy` and
still travels through the recursion in registers.
Apple M4 Max, 12 threads. Ext-mem output verified identical to the
in-memory suffix array on both inputs:
chr21.fa 24.19 s -> 2.48 s wall, 268 s -> 23.3 s CPU
phase 1 22.85 s -> 0.95 s
peak RSS 147 MB -> 151 MB
chr21.0123 3.49 s -> 3.55 s wall (unchanged; no runs to find)
Tests: the run-aware LCP is checked against a byte-at-a-time oracle over
sampled position pairs on homopolymers, wrapped-FASTA blocks and
multi-period texts, plus detection shape (sorted, disjoint, period claim
actually holds) and `max_context` behaviour inside a run.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Adds the ext-mem phase breakdown that located the cost (94% of the chr21 FASTA run in phase 1), the before/after table, and the detail worth keeping: a homopolymer detector would not have worked, because in 60-column wrapped FASTA the longest single-byte run is 60. Each line of `N`s ends in a newline, so the real structure is a period-61 repeat spanning 6.6 Mb. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Phase 4 had become the largest ext-mem cost, and its merge CPU (19.2 s
on 80 MB of DNA) exceeded the entire CPU of the in-memory path. The
merge is latency-bound: the tied branch dereferences the text at two
random addresses, and the address for step i+1 is not known until step i
retires, so there is no memory-level parallelism and the hardware
prefetcher cannot see the pattern.
The candidate positions themselves live in the two index arrays, which
are sequential and already in cache, so the addresses several steps
ahead are known even though the dependent loads are not. Issue them as
prefetches, offset by the current boundary LCP `m`, which estimates
where the next scans start.
This is not the prefetch recorded as a negative result in `lcp.rs`. That
one sat inside the strided scan loop, which the hardware prefetcher
already covers. This one targets the random access, which it cannot.
Apple M4 Max, 12 threads, ext-mem, output verified identical to the
in-memory suffix array on both inputs:
chr21.0123, 80 MB phase4 merge CPU 19.20 s -> 12.42 s
phase4 wall 2.17 s -> 1.46 s
total 3.55 s -> 2.90 s CPU 33.8 -> 27.0 s
chr21.fa, 47.5 MB phase4 merge CPU 7.99 s -> 6.76 s
total 2.56 s -> 1.99 s CPU 23.3 -> 18.0 s
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Sweeping `p` directly showed the previous default sitting on the wrong
side of a flat region. Total work is `n log n` either way, since a
smaller `p` moves levels out of phase 4's per-partition cascade and into
phase 1's merge sort, but the constants differ: phase 3 shrinks
quadratically in `p`, and phase 4's cascade does a full pass over its
partition per level.
Peak RSS is set by phase 4 holding `4 x threads` partitions of `n / p`
records at once, so it only starts growing once `p` is small enough for
that product to rival the text. Measured on chr21 forward ++ revcomp
(80 MB), 12 threads:
p total peak RSS
48 2.58 s 987 MB
96 2.49 s 538 MB
306 2.85 s 282 MB
612 2.80 s 202 MB <- 131072
1224 3.05 s 205 MB <- 65536 (previous default)
128Ki is the largest step that costs nothing in memory: same peak RSS,
~8% less wall. Going further trades real memory for speed, which is the
opposite of what this path is for, so it stays available through
`ExtMemOpts::subproblem_count` rather than becoming the default.
At genome scale this changes nothing. `PHASE1_MAX_PARTITIONS` already
binds for any n above ~1 GB, so GRCh38 still gets p = 8192.
Output verified identical to the in-memory suffix array on both inputs.
chr21.0123, 80 MB 3.01 s -> 2.77 s
chr21.fa, 47.5 MB 2.06 s -> 1.86 s
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Phase 1 had become the largest external-memory phase (1.20 s of 2.54 s on 80 MB of DNA) and was already parallelising at ~97%, so the way forward was less work rather than more threads. It was sorting every subarray from singletons, which is the case the merge kernel handles worst: each leaf merge starts at `m = 0` and orders two suffixes by scanning the text at two random addresses. Sorting by the packed key first resolves the leading `k` symbols with no text access at all (32 symbols for DNA at 2 bits each), and yields the LCP between adjacent runs for free from `(key_a ^ key_b).leading_zeros()`. Only suffixes agreeing through all `k` symbols reach the merge kernel, on the short slice they occupy. This is the bounded-memory counterpart to the prefix doubling in `build_in_memory`. Doubling itself is not available here: it needs a rank for every position in the text, which is exactly the memory this path refuses to spend. A fixed-depth key needs none. The cross-run LCP is capped by both suffixes' lengths. Padding can agree with a real `0` symbol past the end of the shorter suffix, so the raw `leading_zeros` count can overstate it, and a wrong LCP would silently corrupt the order at the next merge level. Gated on the same conditions as the other fast paths (`u8` symbols, plain lexicographic comparator, unbounded `max_context`) and falls back to `merge_sort` otherwise. The alphabet scan that picks the field width runs once per build, not once per subarray. Apple M4 Max, 12 threads, output verified identical to the in-memory suffix array on both inputs: chr21.0123, 80 MB phase1 1.201 s -> 0.316 s total 2.54 s -> 2.11 s chr21.fa, 47.5 MB phase1 0.898 s -> 0.415 s total 1.86 s -> 1.65 s Peak RSS 205 MB -> 233 MB: the key vector is 16 bytes per record over one subarray per worker. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
`CascadeWorkspace::merge_one_level` walked its run pairs sequentially, even though the pairs at a level are independent and write to disjoint destination ranges. The only thing forcing the order was the running `src_off` / `dst_off`, and both are prefix sums, so they can be computed up front and each pair handed its own sub-slices. This is what capped phase 4's parallel efficiency. Partition-level parallelism (`4 x threads` at once) hides it while there are many partitions in flight, but each partition's cascade ends in a single 2-way merge over the whole partition, and those tails serialise. Apple M4 Max, 12 threads, output verified identical to the in-memory suffix array on both inputs: chr21.0123, 80 MB 2.11 s -> 2.02 s chr21.fa, 47.5 MB 1.65 s -> 1.55 s A grain-size threshold was tried, on the theory that the small early levels would not pay for their rayon tasks. It measured worse on both inputs (2.22 s and 1.63 s), so the pairs are merged in parallel at every level. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Emitting is Theta(n) and single-threaded by construction: the caller's closure is `FnMut` and its ordering is the point. But it also did not overlap anything, so once per chunk every worker sat idle while the main thread drained the merged results. A scoped producer now merges chunk c+1, itself rayon-parallel, while the main thread emits chunk c. The channel bound of one keeps at most two chunks resident, so the transient cost is one extra chunk of merged positions rather than the unbounded queue an unsynchronised producer would build. This is not the existing `ordered_phase4_emit` path, which coordinates at *partition* granularity through an mpsc channel and a `BTreeMap` and measured slower than plain collect-then-emit. Here the producer hands over whole chunks that are already in order, so the consumer only drains them and no reordering structure is needed. That path is left untouched behind its opt-in flag. Apple M4 Max, 12 threads, output verified identical to the in-memory suffix array on both inputs: chr21.0123, 80 MB 2.02 s -> 1.82 s peak RSS 233 -> 240 MB chr21.fa, 47.5 MB 1.55 s -> 1.45 s peak RSS 182 -> 194 MB Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
A partition arrives as `p` sorted sub-subarrays, and the cascade merged
them pairwise in `log2(p)` levels, each a full LCP-enhanced pass. That
was the largest single cost left in the ext-mem build: 15.4 CPU-seconds
of a 15.1-second run on 80 MB of DNA.
When a packed key applies, discarding that sortedness and re-sorting the
partition outright is much cheaper. One key sort resolves the leading
`k` symbols with no text access, and only suffixes agreeing through all
of them reach the merge kernel, so `log2(p)` passes collapse into one.
Gated on the run table being empty, which is the point worth recording.
A long periodic run is exactly a stretch where a fixed-depth key
resolves nothing, since every suffix inside it shares the whole key, so
the re-sort would hand the merge kernel one enormous tied group and
throw away ordering phase 1 had already established. Measured
unconditionally it was a large regression on the `N`-heavy input:
cascade key re-sort
chr21.0123 (no runs) 1.95 s 1.65 s
chr21.fa (6.6 Mb N) 1.45 s 2.23 s
So the run table, which already exists to make scans cheap, doubles as
the predicate for whether a fixed-depth key can be expected to pay off.
Apple M4 Max, 12 threads, output verified identical to the in-memory
suffix array on both inputs:
chr21.0123, 80 MB 1.85 s -> 1.65 s phase4 merge CPU 15.4 -> 7.8 s
chr21.fa, 47.5 MB 1.45 s -> 1.47 s (unchanged; takes the cascade)
Peak RSS on the run-free input rises 240 MB -> 285 MB: the key vector
and the sort buffers are live for each of the `4 x threads` partitions
in flight. Still an order of magnitude under the in-memory path's
2.2 GB, which is the comparison that matters for choosing this path.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The field width was chosen from the largest byte value in the text
rather than from how many distinct symbols it uses, and for the most
natural input format those are very different numbers. A plain ACGT
sequence has four symbols, but its largest byte is `'T'` (84), so the
packer used 8-bit fields and fit 8 symbols per key instead of the 32
that four symbols allow.
Rank the bytes that actually occur onto a dense code range first. The
map is monotone by construction, since codes are assigned in ascending
byte order, so a packed key stays order-preserving and the zero-padding
argument carries over unchanged: code 0 is still the minimum.
Apple M4 Max, 12 threads, chr21 forward ++ revcomp written as ASCII
ACGT (80 MB, four symbols, largest byte 84):
before after
seed sort 0.480 s 0.281 s
doubling 0.867 s 0.488 s
total 1.404 s 0.827 s -41%
The two inputs benchmarked so far both happened to hide this. One is
pre-encoded to codes 0..3 and is already dense; the other is dominated
by its `N` runs, where key depth is irrelevant because every suffix
inside a run shares the whole key. A FASTA-derived ACGT text is the
common case that neither covered.
Cross-check: the suffix array of the ASCII text is identical to the
suffix array of the 0..3-coded text, as it must be, the two encodings
being order-isomorphic. Both the in-memory and ext-mem paths agree.
Tests: the field width follows alphabet size rather than byte value
across four alphabets, and the remap is asserted monotone.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Building a key was a chain of `k` dependent shift-or-lookup steps, 32 of
them for 2-bit DNA fields, and that chain set the cost of the seed sort
in every path that uses a packed key.
Two changes remove it. The alphabet map is applied to the whole text
once, up front, so the inner loop no longer carries a dependent table
load. And the fields are gathered by a binary-tree SWAR shuffle: each
step folds neighbouring fields together and halves the stride, so eight
symbols cost three shift-or-mask pairs rather than eight dependent
steps. The load is big-endian, which puts the text's first byte in the
result's most significant field, the order the key already needed.
The ranked copy is only materialised when the identity map does not
already rank the text. A `0..3` DNA encoding is already dense and pays
nothing; an ASCII text pays one byte per symbol, which is what buys the
narrower fields in the first place.
Note the group concatenation assigns the first group rather than
shifting it in. With 8-bit fields there is exactly one group and the
shift distance would be 64, which is not legal for `u64`: release builds
mask it to zero and happen to produce the right answer, debug builds
panic. This was caught by running the tests in debug, which is why the
CI workflow does.
Apple M4 Max, 12 threads, 80 MB chr21 forward ++ revcomp, output
verified identical on all three encodings:
seed sort total
ASCII ACGT 0.281 s -> 0.177 s 0.827 s -> 0.702 s
coded 0..3 0.278 s -> 0.168 s 0.802 s -> 0.719 s
Peak RSS is 2.21 GB for the coded input, unchanged, and 2.29 GB for the
ASCII one, the difference being the ranked copy.
Tests: a new case checks the SWAR gather against the obvious shift-or
loop for every field width, every alignment and every tail length, and
the monotonicity of the alphabet map is now asserted through `key_at`
rather than through the internal table.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Instrumenting the rounds (`CAPS_SA_PROFILE=1` now prints one line each)
showed where the time actually went, and it was not where I assumed.
On 80 MB of DNA the first round has 3.5 million tied groups averaging
**four elements**:
round depth=32: 14193102 tied in 3526079 groups (avg 4.0)
round depth=64: 9929726 tied in 2924078 groups (avg 3.4)
round depth=128: 7480152 tied in 2543160 groups (avg 2.9)
At that size neither the sort nor the rank probes dominate. The
bookkeeping around them does: one heap allocation per group for the key
vector, and a sequential `split_at_mut` chain over all 3.5 million
groups to hand each one its sub-slices.
Both go. `Scatter` already encodes "disjoint ranges, one owner each",
which is exactly the property the groups have, so a group takes its own
sub-slices directly and the sequential prepass disappears. Groups up to
32 elements build their key vector in a stack buffer, so nearly every
group avoids the allocator. `flat_map_iter` replaces collecting a
`Vec<Vec<_>>` of mostly-empty vectors.
Apple M4 Max, 12 threads, output verified identical on all three
encodings and both paths:
doubling total
chr21.0123 0.502 s -> 0.234 s 0.733 s -> 0.566 s
ASCII ACGT 0.764 s -> 0.611 s
The first round alone goes from 0.124 s to 0.053 s. The `N`-heavy FASTA
is unchanged: its groups are large, so it was never paying the
per-group overhead.
The per-round log stays in, behind `CAPS_SA_PROFILE`. It is what located
this, and the group-size distribution is the thing worth looking at
first for any further work here.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Addresses the blockers @rob-p raised in COMBINE-lab#7. **Signed symbols.** `seed_params` accepted any `Symbol` one byte wide and reinterpreted it as `u8`, but `Symbol` is implemented for `i8` too and a packed key orders its fields as unsigned: `-1` has byte `0xFF` and sorts above `1`. The in-memory doubling guard already required exact `u8`; the packed-key paths now match it. Rob's reproducer `[-1, 0, -1, 1, -2, 0, 1, -1, 0, -2, 1, 0]` is in the tests, alongside external-memory and sample-sort cases at several partition counts and an `i8::MIN`/`i8::MAX` case. Reverting the guard fails exactly the two end-to-end tests and leaves the in-memory one passing, matching the report that only the packed-key paths were affected. `runs.rs` deliberately keeps its one-byte-wide check: a run table is built from byte *equality*, which coincides with value equality for `i8`, and ordering is recovered by the caller through `S: Ord`. **Run skipping taxed the common path.** Consulting the table cost up to three binary searches before any comparison, and nearly every LCP call in real sequence mismatches within a few symbols and never reaches a run. Probe with the ordinary bounded scan first, and consult the table only once a match has survived 256 symbols. The probe is not wasted: whatever it matches counts towards the answer. Measured on a fixture matching the one in COMBINE-lab#7 (parsed chr21 forward ++ revcomp, `N` kept as a symbol, ACGT-start filter, 80,177,238 retained positions, 12 threads), with byte-identical output in all three cases: run table disabled 1.97 s table consulted eagerly 2.40 s (+22%) probe first 2.02 s (+2.5%, within noise) And it keeps the benefit where runs matter: the raw FASTA takes 43.8 s with the table disabled against ~1.5 s with it. **Packer built before its guards.** `Packer::new` may materialise a ranked copy of the whole text, and on a segmented or context-bounded build that copy can never be used, since every packed-key path declines those. It is now constructed only after the `plain_lex_len` and `max_context` guards pass. **Docs.** The claim that periods up to 64 also cover satellite arrays was too broad: the canonical alpha-satellite monomer is 171 bases and measures at parity with no detection. Both the module docs and the README now state the measured coverage and its limit. Also adds `--filter-acgt` to the bench CLI, so the STAR-shaped filtered build this all turns on can be benchmarked directly. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The packed-key paths declined every segmented build, which is the comparator a splice-junction index actually uses, so none of this work reached it. The obstacle was that zero-padding a short suffix encodes shorter-is-smaller, while STAR's `boundary_order` is `lim_b.cmp(&lim_a).then(p_a.cmp(&p_b))` — longer-is-smaller. That is expressible in the key once the provider says which side it wants. `LimitProvider::boundary_rank` returns `ShorterFirst` or `LongerFirst`, defaulting to `None` so existing implementations keep today's behaviour. A segmented key then packs `min(k, lim_at(p))` symbols — never reading into the next segment — and pads with a reserved sentinel below every real code under `ShorterFirst`, above every real code under `LongerFirst`. The position tie-break does not need to be in the key. Keys that tie defer to `boundary_order` itself, so the key only has to avoid contradicting the convention, never to reproduce it. The sentinel needs one code above the alphabet, so a segmented build sizes its field to hold `alphabet` rather than `alphabet - 1`, and `Packer::new` takes that as a parameter: a plain build must not pay for it, since widening the field halves the symbols a key carries. Where no spare code fits (256 distinct symbols in 8-bit fields) the path declines. Tests compare against the provider's own comparator over random segmented texts on the ruSTAR alphabet, under both conventions, with and without the ACGT-start filter, at three partition counts. Two notes on how those tests are written, both learned the hard way: - They assert the *property* (permutation, and no adjacent pair out of order) rather than equality with a canonical answer. `SegmentedText`'s default `boundary_order` returns `Equal` for suffixes that end together with equal content, so their relative order is genuinely free and a stable-sort oracle is not a valid reference. - Sensitivity was checked by breaking `key_at_bounded` deliberately. A *constant* wrong key still passes, because it collapses everything into one tied group that the fallback sorts correctly — so that check proves nothing. An order-*inverting* key does fail the test, which is what establishes the fast path is both taken and verified. Not benchmarked end to end on an annotation-shaped fixture: that needs rustar-aligner's junction pipeline, which I cannot run here. @rob-p, if the GENCODE harness from COMBINE-lab#7 is available this is the change that should finally move phase 1 on it. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
… branch Measures the case @rob-p correctly identified as the one none of this work reached: segmented text, STAR's spacer-as-largest `boundary_order`, ACGT-start filter, external-memory construction. `bench/gsj_fixture.py` builds the fixture from a FASTA and a GTF the way a STAR index is laid out — genome sequence, then one 2*overhang flank per distinct junction, then the reverse complement — and `gsj_bench` runs it through `build_ext_mem_for_positions_with` under a `LongerFirst` provider. Apple M4 Max, 12 threads, chr21 with its GENCODE annotation: 95,409,966 symbols, 9,952 segments, 82,167,238 retained ACGT-start positions. segmented keys declined (previous behaviour) 3.84 s segmented keys enabled 2.57 s -33% Output order verified directly against the provider's comparator. This is a *shape* reproduction, not a size one. A genome-wide annotation yields far more junctions than a single chromosome's: this fixture has 4,975 junctions and 9,952 segments against the 698,597 and 1,397,196 of the GENCODE v50 primary-assembly fixture in COMBINE-lab#7. The segmented-key path should matter more, not less, as the segment count grows, but that is a prediction and not something these numbers establish. Also runs CI on every branch rather than only `main`. A pull request opened from a fork by a first-time contributor does not run workflows until a maintainer approves them, which is why the API showed no check runs at the pinned tips. Building on push makes the fork produce evidence that can be linked from the PR. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
@rob-p asked, correctly, what happens to *memory* once a realistically sized junction library is included, and suspected an interaction with prefix doubling and radix sorting scaling in the number of distinct elements rather than in `n`. The answer is that there is one, and it cuts against making this a default. Fixture: chr21 plus a 698,597-junction library, 372,858,766 symbols, 1,397,196 segments, 320,856,244 retained ACGT-start positions. That is the same junction and segment count as the GENCODE v50 primary-assembly fixture in COMBINE-lab#7, and within 0.4% of its symbol count. Apple M4 Max, 12 threads: segmented keys off segmented keys on phase 1 9.544 s 3.313 s phase 4 11.616 s 11.421 s total 22.718 s 16.256 s -28% peak RSS 3.45 GB 5.33 GB +54% So the segmented key does reach the splice-junction case and does cut phase 1 by 65%. But it buys that with half again as much resident memory, from the ranked text copy and the per-subarray key vectors, both proportional to the input. For a constructor whose argument against libsais and against STAR is memory, that is not a trade to make unconditionally. This is the same shape of objection raised against #4's dense-subset heuristic, and it wants the same answer: a memory budget the caller sets, not a silent default. Not implemented here — recording the measurement first so the decision is made against a number. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…on target
The RSS number I reported was wrong, and it was wrong in the direction
that changes the recommendation, so correcting it first.
I compared one run against one run. Replicating both sides three times
shows peak memory is unchanged:
segmented keys off segmented keys on
total 23.5-24.0 s 15.8-16.2 s -33%
peak RSS 3.22-3.94 GB 3.21-3.23 GB
The earlier "3.45 GB -> 5.33 GB, +54%" was one outlier measured against
another. There is no memory cost to weigh against the speedup, so the
argument I made for keeping segmented keys opt-in on memory grounds does
not hold.
Also sweeps the partition count on the same segmented fixture, which is
the re-evaluation the review asked for: the 128 Ki target was tuned on a
plain workload and needed checking on a segmented one.
p = 2845 (128 Ki target, current default) 15.26 s 3.22 GB
p = 5690 (64 Ki target, previous default) 19.09 s 3.71 GB
p = 11380 26.97 s 5.19 GB
It holds up: the current default is both faster and smaller than the one
it replaced, and the trend continues in the same direction.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The review objected that at `subset >= text / 8` this silently constructed the complete suffix array plus an `n`-element membership structure and full rank arrays, on the caller's behalf and without saying so. The objection is right: those arrays are sized by the *text*, not by the subset, so the rule made a large allocation decision that only the caller can actually make. Replaced by `Opts::subset_full_sa_budget`, the peak extra bytes the subset path may spend. `None` is the default and never takes the trade; `Some(budget)` allows it when the estimated footprint fits, which the docs spell out as `n * (3 * size_of::<I>() + 9)`: three index-wide arrays for the suffix array, the ranks and the round scratch, eight bytes of key per position, and a one-byte membership flag. Tested at three text sizes with the budget set exactly to the estimate, one byte below it, and left at the default: all three produce the same answer, so the gate changes resource use and never correctness. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
First step towards the subsumed-LCP idea @rob-p described, and an honest one: the mechanism is right and the measured effect is small. The merge already exploits the LCP triangle inequality (`lcp(a,c) = min(lcp(a,b), lcp(b,c))` when the two differ) for *adjacent* pairs — that is exactly its three-case rule, which resolves two of three cases with no text access. What it did not exploit is a known LCP arriving from outside the merge: a tied group coming out of a packed-key sort agrees on at least `k` symbols by construction, and the fallback merge rescanned them from zero. `merge_from` takes that as a `base`. It costs one substitution in the existing invariant: the three-case rule is stated against the last-output element, initialised to the empty string, and is now initialised to the length-`base` prefix every element shares. `lcp(B, s) = base` for all `s`, `B` precedes every element, and `lcp_x[0] = base` is what the first iteration reads. The loop body is untouched, because every case in it argues about relative offsets and none mentions zero. `merge` stays as a `base = 0` shim, so the cascade and phase-1 call sites are unchanged. Measured, 12 threads, output byte-identical: chr21 coded, ext-mem 1.69 s -> 1.61 s ruSTAR-shaped filtered 2.02 s -> 2.05 s (no change) Small, and worth saying why: after a 32-symbol key sort the tied groups are short, and the merge's own three-case rule already starts subsequent scans at the right offset. Only the first comparison in each group saves the rescan. The gain would grow with the key depth relative to typical LCPs. The value here is the mechanism rather than the number. A general subsumed-LCP cache — inferring `lcp(a,c)` from cached `lcp(a,b)` and `lcp(b,c)` for non-adjacent pairs — needs exactly this: a merge that can be told a lower bound it did not derive itself. Run skipping is the other special case already in the tree, deriving its bound from periodic structure instead. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Closes the structural gap behind the subsumed-LCP discussion, and comes
with a finding that changes what is left to do.
**The finding.** In a 2-way merge the LCP triangle inequality is already
saturated. After computing `lcp(a_i, b_j) = L` and advancing A, the next
pair is `(a_{i+1}, b_j)`, and `lcp(a_i, a_{i+1})` is already in the
source LCP array; combining them is exactly what the three-case rule
does. There is no residual inference to harvest inside a merge. The
places a cached LCP can still pay are *outside* one: a bound arriving
from elsewhere, which is what `merge_from`'s `base` and run skipping
both supply.
**The gap.** Prefix doubling answers comparisons from ranks and never
computes an LCP, which is the structural reason the external-memory path
could not be routed through it — that path needs the array the merge
kernel yields as a byproduct. And the merge's array was internal, so no
caller could obtain one at all.
Kasai's algorithm closes both: one linear pass from the suffix array.
The bound is the point, and it is exactly what the scanning merge lacks
— `h` falls by at most one per position and rises only while matching,
so symbol comparisons total at most `2n` however repetitive the text is.
Tested against a naive per-pair scan on the inputs that make the naive
version expensive: long runs, period-3 and period-61 text, three
alphabet widths, plus the empty and single-symbol cases.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Wall times on the external-memory path drift 30% between sessions on the same machine and commit. Diagnosed after a phantom regression that bisected to a commit predating it: 32 GB of benchmark *output* had accumulated in the working directory, and each build streams several GB of bucket data through the page cache. Deleting the outputs restored the original numbers. with 32 GB resident 20.5-21.1 s after deleting them 15.8-17.5 s, 154-166 s CPU The temp files are not implicated and do not leak: `$TMPDIR` measured 454 MB and 801 entries both before and after a run, so the pooled anonymous buckets are released as intended. Writes down the four rules that follow: send output to /dev/null unless checking correctness, interleave A against B rather than batching each, report CPU alongside wall, and treat a single reading as a hypothesis. That last one is not abstract. Two figures published in this file came from unreplicated single runs and both were wrong: a claimed +54% RSS that replication showed to be flat, and an 11.4 s phase-4 time that was 15-16 s on repetition. The interleaved comparison of the segmented-key change gave the correct 0.67 ratio even on a session where neither side's absolute number was right. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Phase 3 asks for `p` upper bounds per subarray, once per pivot, and the pivots are sorted, so the answers are non-decreasing. Each search ignored that and bisected the whole array: at `p = 612` over 131 K-record subarrays that is ~17 suffix comparisons per pivot no matter how close the answer is to the previous one. Galloping forward from the previous split costs one comparison when the bracket is empty, which it usually is once `p` is large enough that consecutive pivots land within a few records of each other. Measured, 12 threads, output verified identical to the in-memory suffix array: phase 3 0.394 s -> 0.362 s on chr21. That is 8%, not the 8x the comparison count alone would predict, and the gap is the useful part: phase 3 is bound by data movement, not by searching. It writes every one of the n records into its partition bucket — 640 MB at this input size — and 0.36 s for that is roughly 3.5 GB/s, which is the real floor. An earlier attempt to accelerate the same searches with packed keys was reverted for measuring nothing, and this explains why: both were optimising the part that was not the cost. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…cket count Adds per-pass timing to the seed sort behind `CAPS_SA_PROFILE`, which locates its cost precisely on 80 MB of DNA: histogram 0.058 s prefix sum 0.001 s scatter 0.080 s bucket sort 0.127 s <- 48% The per-bucket comparison sort dominates, so the obvious lever is the bucket count: more buckets, smaller and more cache-resident sorts. The 11-bit choice had been reasoned about (write-combining state stays inside a core's private cache) but never measured. Measured now, and it makes no difference. A first sweep suggested 13 bits was 18% better, but that was single runs. Interleaving three pairs per the methodology this repo now documents: 11 bits 0.582 0.654 0.613 mean 0.616 s 13 bits 0.640 0.623 0.595 mean 0.619 s So the bucket sort is 48% of the seed sort and is not reachable by changing how many buckets there are. Left at 11 bits, and recorded here so the sweep is not repeated. This is the fifth time in this branch that a single reading pointed the wrong way. The interleaved comparison has been right every time. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Per-phase thread scaling, 1 to 12 threads on 80 MB of DNA:
threads phase 1 phase 3 phase 4 total
1 2.088 s 0.756 s 6.413 s 9.288 s
12 0.289 s 0.361 s 0.954 s 1.640 s
speedup 7.2x 2.1x 6.7x 5.7x
Phases 1 and 4 scale acceptably. Phase 3 saturates by six threads and
gains 14% from six to twelve. It is 22% of the wall here and would be a
larger share on a wider machine, which makes it the thing to fix rather
than the seed sort or the cascade.
Two hypotheses tested, both wrong, both recorded so they are not
retried. Physical file contention: raising `CAPS_SA_N_PHYS` from 12
through 192 does not help and mildly hurts (0.356 s to 0.420 s). Search
cost: galloping cut pivot comparisons roughly eightfold for 8%, and an
earlier packed-key attempt at the same searches measured nothing.
What remains is the write path. Phase 3 reads every record back and
writes every record out again, ~1.3 GB here, through the page cache,
which does not parallelise with threads. 3.5 GB/s is far below this
machine's DRAM, which points at the kernel path rather than memory.
So the fix is not to speed phase 3 up but to not run it: fuse it into
phase 1, which needs pivots before phase 1 rather than after. A cheap
pre-sampling pass over the raw positions supplies them, and correctness
does not depend on their quality — any splitters give a correct sample
sort, only the balance changes. That removes one full write-and-read
round trip of every record. Not attempted here: it restructures the
ext-mem driver and the filtered and positions variants with it, and
that wants its own change with its own verification.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Phase 3 was the worst-scaling part of the external-memory build: 2.1x
on twelve threads against 7.2x and 6.7x for phases 1 and 4, saturating
by six threads and gaining 14% from six to twelve. Two explanations were
tested and both were wrong — raising the physical file count from 12 to
192 does not help, and galloping cut its pivot comparisons roughly
eightfold for 8%. What was left is the write path: phase 3 read every
record back and wrote every record out again, ~1.3 GB at this input
size, through a page cache that does not parallelise with threads.
So do not make phase 3 faster; remove it.
It existed only because pivots were sampled from the already-sorted
subarrays, which forced the order: sort everything and spill it, choose
pivots, read all of it back to distribute. But **pivots do not have to
come from sorted data**. Any splitters produce a correct sample sort;
only the balance of the partitions changes. A cheap pre-pass over the
raw positions supplies them, and phase 1 can then sort and distribute in
one go, never materialising the subarrays at all.
Sampling is by strided blocks of 64 rather than strided singletons: a
`PositionSource` fills a contiguous run cheaply but pays per call, and
the `Filtered` variant especially so. Coverage of the position space is
the same and splitter quality is not sensitive to the difference.
Two spilled copies of every record disappear — one write and one read —
and with them an entire bucket file pool, since there are no subarray
buckets left to hold.
Apple M4 Max, chr21 forward ++ revcomp (80 MB), interleaved A/B, three
pairs, fused winning every pair:
before 1.975 1.733 1.790 s mean 1.833
after 1.671 1.527 1.609 s mean 1.602 -12.6%
Scaling improves as well, which is the point rather than the wall time:
1 thread 12 threads speedup
before 9.288 s 1.640 s 5.7x
after 9.563 s 1.352 s 7.1x
The phase 1 + phase 3 pair went from a combined ~3.9x to 6.3x for the
fused phase alone.
Verified: byte-identical output against the in-memory suffix array on
N-free DNA and on raw FASTA, against the previous ext-mem path on the
ACGT-filtered ruSTAR-shaped fixture, and `--verify` clean on the
segmented 698,597-junction splice-junction fixture.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Contributor
|
Thanks @BenjaminDEMAILLE — the core restructuring is sound, and it produced a clear benefit on the complete ruSTAR input. I integrated it onto I did not merge this PR head verbatim because it contains 40 prerequisite commits from #8. Doing so would have merged unreviewed work and subsumed #10, #11, and #13 before their requested independent reviews. Instead I ported only the phase-0 presampling and fused external-memory sort/distribution design onto the pre-review Full ruSTAR fixtureGENCODE human v50, GRCh38 primary-assembly FASTA plus comprehensive primary-assembly GTF, using ruSTAR's actual FASTA/GTF parsing, genome bin padding, junction preparation/deduplication,
One warm-up, then three interleaved measured pairs:
The preprocessing portion fell from a 107.426 s mean to 102.787 s. The deterministic presample also improved final partition balance: phase 4 fell from 280.993 to 258.097 s. Temporary filesystem writes fell from about 193 GB to 96.5 GB, system CPU time from 102.01 to 55.12 s, and mean peak RSS from about 10.47 to 10.15 GiB. Every full run emitted 6,176,694,310 positions with hash For completeness, the isolated change regressed the smaller focused fixture (14.43 s historical Closing this contributor PR as superseded by the clean |
This was referenced Aug 12, 2026
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
2 participants
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.
Depends on #8. Removes phase 3 from the external-memory build entirely.
Why
Phase 3 was the worst-scaling part of the build:
It saturates by six threads and gains 14% from six to twelve. Two explanations were tested and both are wrong, recorded so they are not retried:
CAPS_SA_N_PHYSfrom 12 through 192 does not help and mildly hurts (0.356 s to 0.420 s).What remains is the write path. Phase 3 read every record back and wrote every record out again — ~1.3 GB at this input size — through a page cache that does not parallelise with threads. 3.5 GB/s is far below this machine's DRAM, which points at the kernel rather than memory.
What
Phase 3 existed only because pivots were sampled from the already-sorted subarrays, which forced the order: sort everything and spill it, choose pivots, read all of it back to distribute.
Pivots do not have to come from sorted data. Any splitters produce a correct sample sort; only the balance of the partitions changes. That is what makes this restructuring far safer than it looks — correctness does not depend on the sample being good.
So a cheap pre-pass over the raw positions supplies them, and phase 1 sorts and distributes in one go. The subarrays are never materialised.
Sampling is by strided blocks of 64 rather than strided singletons: a
PositionSourcefills a contiguous run cheaply but pays per call, and theFilteredvariant especially so. Coverage is the same and splitter quality is not sensitive to the difference.Two spilled copies of every record disappear — one write and one read — and with them an entire bucket file pool, since there are no subarray buckets left to hold.
Measured
chr21 forward ++ revcomp (80 MB), 12 threads, interleaved A/B, three pairs, fused winning every pair:
Scaling is the real result:
The phase 1 + phase 3 pair went from a combined ~3.9x to 6.3x for the fused phase alone, so the gain grows with core count. @rob-p, this should show up more strongly at 32 or 64 cores than it does at 12.
Verified
--verifyclean on the segmented 698,597-junction splice-junction fixture under STAR's boundary order.Note
Single-thread time is slightly worse (9.288 s to 9.563 s): the pre-sample is extra work that the old flow got for free from data it was sorting anyway. It pays for itself by two threads and grows from there.
🤖 Generated with Claude Code