1/4 perf: radix-seeded prefix doubling for the in-memory suffix array - #3
Open
BenjaminDEMAILLE wants to merge 5 commits into
Open
1/4 perf: radix-seeded prefix doubling for the in-memory suffix array#3BenjaminDEMAILLE wants to merge 5 commits into
BenjaminDEMAILLE wants to merge 5 commits into
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>
BenjaminDEMAILLE
force-pushed
the
stack/1-in-memory-fast-path
branch
from
August 13, 2026 17:38
72db165 to
cda07cb
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.
Introduces the in-memory fast path, plus the test and CI scaffolding the rest of the stack depends on.
Why
The report that started this was "caps-sa is ~18x slower than libsais on chr21". Reproducing it turned up something more useful: the two numbers were measured on different inputs, and the gap is one specific pathology.
Dividing CPU time by
n log2 nmerge steps, same kernel and same machine:16x apart, and it is entirely scan length. Every leaf merge starts with
m = 0, so a tied comparison scans the whole shared prefix. Genome FASTA carries ~6.6 Mb ofN(period-61 once 60-column wrapping is included), so one comparison can scan millions of bytes.What
Sort by a packed fixed-depth key, then resolve the rest by prefix doubling on ranks. After the seed no comparison reads the text again, so a megabyte-long
Nrun costs what random DNA costs.The seed orders by
(key, min(n - p, k)). The second component is required: zero-padding makes a short suffix share a key with any suffix continuing in zeros, and0is a real symbol in every DNA encoding. Without it[0, 0]leaves two positions permanently tied and doubling cannot terminate.Guards (all soundness conditions, all defaulting to declining): unbounded
max_context,LimitProvider::plain_lex_lenreporting the full text, andSexactlyu8. The existingStarConventiontest impl inheritsNonefrom the new trait method and stays on the merge kernel untouched.Also here
[profile.release]pinned.bench/README.mdcredited fat LTO to a parent workspace that is not in this repo, so every build since the crate went standalone usedlto = false, codegen-units = 16.verify_sa, anO(n)independent check that does not depend on LCP length.Numbers (Apple M4 Max, 12 threads)
N-free, 80 MBOutput byte-identical to the merge kernel on both inputs (127.7 M entries).
🤖 Generated with Claude Code