2/4 perf: MSD counting-sort seed, and cover subsets + in-memory sample sort - #4
Open
BenjaminDEMAILLE wants to merge 8 commits into
Open
Conversation
Open
4 tasks
Contributor
Author
|
Part of the stack tracked in #7. |
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>
`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>
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>
BenjaminDEMAILLE
force-pushed
the
stack/2-faster-seed-and-coverage
branch
from
August 13, 2026 17:38
31c1bb4 to
14f3107
Compare
Contributor
Author
|
Rebased onto 0.7.0 and re-measured against it. Details in the updated Two things worth flagging on the other two parts of the stack, #6 and #8: |
Contributor
Author
|
One correction to the numbers below, which I would rather state than have you
Prefix doubling needs a rank for every position, so this is inherent to the That is a trade for you to accept or refuse explicitly, not one to make |
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.
Two follow-ups to the fast path: a cheaper seed, and coverage of the two entry points stack 1 left behind.
Seed sort: comparison sort → MSD counting sort
The seed was
par_sort_unstableover a materialisedVec<(u64, u32, I)>, which was both the largest remaining cost and the peak-memory driver.I = u32against the 16 a record costs.11 bits (2048 buckets) keeps the write-combining state near 512 KB, inside a core's private cache; 16 bits would need 16 MB of open write lines.
The visible-length tie-break no longer needs storing — only the last
k - 1positions can have a visible length belowk, so it is a function of position alone.Subsets and in-memory sample sort
Doubling cannot be restricted to a subset directly: a round compares
rank[p + d]and that successor is generally not in the subset, so ranks must be defined for every text position. Building the whole array and filtering it is oneO(n)pass instead. Below one eighth of the text it declines and the merge kernel runs — that ratio is a performance heuristic, unlike the guards it sits behind. It also declines on duplicate or out-of-range positions, since the output is a permutation of the input multiset.build_in_memory_sample_sortexists to sort in RAM, so where doubling applies it is strictly better.build_ext_memdeliberately does not get this: bounding peak memory is its purpose, and routing it through an in-memory algorithm would defeat exactly that. Stack 3 deals with it properly.Numbers (Apple M4 Max, 12 threads)
--in-mem-ssTests: 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 returning a wrong answer.
🤖 Generated with Claude Code