Skip to content

3/4 perf: external memory, 24.2 s -> 1.56 s on chr21 FASTA - #5

Open
BenjaminDEMAILLE wants to merge 22 commits into
COMBINE-lab:mainfrom
BenjaminDEMAILLE:stack/3-external-memory
Open

3/4 perf: external memory, 24.2 s -> 1.56 s on chr21 FASTA#5
BenjaminDEMAILLE wants to merge 22 commits into
COMBINE-lab:mainfrom
BenjaminDEMAILLE:stack/3-external-memory

Conversation

@BenjaminDEMAILLE

Copy link
Copy Markdown
Contributor

Part 3 of a four-PR stack. Merge in order: 1, 2, 3, 4.

I do not have write access to this repository, so the parts cannot use each
other as base branches; each targets main directly. Consequence: this diff still shows parts 1-2's commits.
Every diff collapses to just its own work as soon as the part before it merges.

Stack: 1/4 in-memory fast path · 2/4 seed + coverage · 3/4 (this) · 4/4 key packing

build_ext_mem stays on the merge kernel by design — prefix doubling needs a rank for every text position, which would defeat the memory bound the path exists to provide. So it was still paying the full scan cost. Profiling put 94% of a chr21 FASTA run in phase 1:

                     phase1     phase2   phase3   phase4    total
chr21.0123  before    1.15 s    0.005 s  0.215 s  2.080 s   3.49 s
chr21.fa    before   22.85 s    0.058 s  0.317 s  0.948 s  24.19 s

Six changes, each measured separately:

change chr21.0123 (80 MB) chr21 FASTA (47.5 MB)
baseline 3.49 s 24.19 s
skip periodic runs 3.55 s 2.48 s
prefetch the next candidates' text 2.90 s 1.99 s
subarray target 64Ki → 128Ki records 2.54 s 1.86 s
seed phase-1 subarrays with the packed key 2.11 s 1.65 s
merge cascade run pairs in parallel 2.02 s 1.55 s
pipeline the emit against the next merge 1.85 s 1.45 s
re-sort run-free partitions by key 1.65 s 1.47 s

The interesting ones

Run skipping. 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 reaches e, so lcp(a, b) >= e - b is known in O(1) with nothing scanned. A homopolymer detector would not have worked: in 60-column wrapped FASTA the longest single-byte run is 60, because each N line ends in a newline. The real structure is a period-61 repeat spanning 6.6 Mb. Detection is two-stage, so run-free texts pay one sampling sweep and get an empty table.

Prefetch. The merge is latency-bound: the tied branch dereferences the text at two random addresses and the next address is unknown until the current step retires. The candidate positions live in the index arrays, which are sequential and cached, so addresses several steps ahead are known. This is not the prefetch recorded as a negative result in lcp.rs — that one sat inside the strided scan, which hardware already covers.

Subarray target. Sweeping p showed the old default on the wrong side of a flat region. 128Ki is the largest step that costs nothing in memory; going further trades real memory for speed, so it stays caller-controlled.

Partition re-sort, gated on the run table. Unconditionally it was a large regression on the N-heavy input (1.45 s → 2.23 s): a long run is exactly where a fixed-depth key resolves nothing, so the re-sort hands the merge kernel one enormous tied group and discards ordering phase 1 had established. The run table, which exists to make scans cheap, doubles as the predicate for whether a fixed-depth key will pay off.

Negative result, recorded

A grain-size threshold on the parallel cascade levels measured worse on both inputs (2.22 s and 1.63 s), so every level goes parallel.

Peak RSS 214 MB → 285 MB, still an order of magnitude under the in-memory path. Ext-mem output verified identical to the in-memory suffix array on both inputs at every step.

🤖 Generated with Claude Code

BenjaminDEMAILLE and others added 22 commits August 11, 2026 21:33
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>
@BenjaminDEMAILLE

Copy link
Copy Markdown
Contributor Author

Part of the stack tracked in #7.

@BenjaminDEMAILLE

Copy link
Copy Markdown
Contributor Author

Status against 0.7.0: this does not rebase as it stands, and most of it should
not. Of its eleven commits, three are already settled: the merge-loop prefetch
is on main via #12, run skipping was closed as #10, and the parallel cascade
was closed as #13. The rest patch the phase 1 and phase 3 that 0.7.0 replaced
with the fused sort-and-distribute, and main already gallops from the
previous split in the pivot search.

What survives that is the phase-1 packed-key seed, re-derived against 0.7.0 in
#15 and measured on the annotated human genome there: 454.7 s to 342.3 s, minus
24.7%, peak RSS flat.

Leaving this open but stale rather than force-pushing something unreviewable.
Happy to close it if you would rather the re-derived piece stand alone.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant