genomeGenerate: libsais suffix-array builder, selected by --limitGenomeGenerateRAM

[BenjaminDEMAILLE]

genomeGenerate: build the suffix array with libsais when it fits, selected by --limitGenomeGenerateRAM.

Corrections to the previous revision of this description

Three claims in the earlier body did not hold against the diff. They are fixed in the branch rather than in the prose.

--limitGenomeGenerateRAM did not select anything. The body said selection was driven by the flag and that this turned it "from warned-and-ignored into a real flag". It did not: build_libsais was called from nothing but its own tests, generate_streaming went to build_streaming unconditionally, and lib.rs still logged "accepted but not enforced". The builder was dead code and the flag was inert. It is wired now (0e34aab), and the warning is gone.

The streaming verification level did not exist. The body listed a check described as "build_libsais driven exactly as GenomeIndex::generate_streaming does, matching the caps-sa streaming byte stream (the path that ships, not a unit-test-only one)". The streaming test in the file compares caps-sa against caps-sa and says nothing about libsais. That test exists now (a9aaaa9), and it is only meaningful because the dispatch above makes it the shipping path.

The RAM estimate was half the truth. libsais_peak_bytes took n_genome and multiplied by 17, but the text libsais sorts is both strands. Any caller passing n_genome would have understated the requirement by a factor of two and handed a mammalian genome to the in-memory builder at the default limit. The argument is now n_entries, the call site doubles, and a test picks a limit that separates the two readings instead of one where they agree.

I also said the multi-genome benchmark was blocked on hardware I did not have. That was wrong; the machine has 128 GB and 16 cores. The full three-tier benchmark CONTRIBUTING asks for is below.

What changed

sa_build::build_libsais builds the suffix array with libsais instead of caps-sa, choosing the alphabet width automatically (u8/u16 for few segments, i32 large-alphabet for many-segment genomes).

GenomeIndex::generate_streaming picks between the two from --limitGenomeGenerateRAM: libsais when 17 × 2N fits the limit, caps-sa's external-memory path otherwise. --genomeSAsparseD above 1 always takes caps-sa, the only builder with a sparse arm. Both the chosen path and the two numbers are logged, so the decision is visible in the run rather than inferred:

Building suffix array with libsais (estimated peak 0.5 GB, --limitGenomeGenerateRAM 33.3 GB)...
Building suffix array with caps-sa external memory (libsais would need an estimated 106.7 GB, --limitGenomeGenerateRAM 33.3 GB)...

The new dependency

libsais 0.2.0, MIT OR Apache-2.0, wrapping IlyaGrebnov/libsais through feldroop/libsais-rs. The C is vendored in libsais-sys and compiled by cc: no system library to install, no bindgen, no libclang, no pkg-config. Declared default-features = false, which turns off the crate's default openmp feature, so the C builds single-threaded and needs nothing outside the vendored sources. Adds bytemuck and either, plus cc as a build dependency.

CONTRIBUTING.md asks for a non-Rust dependency to be raised in an issue before the PR. This predates that rule; the issue is #162, with the trade-offs and the alternatives laid out. That decision gates this PR, and if the answer is no it closes.

Why the output cannot change

A suffix array is unique for a given text. build_libsais sorts the same per-segment sentinel-transformed text the caps-sa sentinel arm uses, so its packed output is byte-for-byte identical.

Checked at four levels, the last two of which exercise the path that ships:

Level	Test
unit, three fixtures	build_libsais_matches_caps_sa_{single_chr,multi_chr_with_n,repetitive}
i32 large-alphabet arm forced	build_libsais_i32_matches_caps_sa_{single_chr,multi_chr_with_n}
streamed bytes, as generate_streaming writes them	libsais_streamed_matches_the_caps_sa_stream_byte_for_byte
real genomes	the benchmark below, SA and SAindex compared with cmp

Benchmark

Three genomes, as CONTRIBUTING asks: yeast, chromosome-scale, mammalian-scale. Same machine (M-series, 16 cores, 128 GB), --runThreadN 8, warm cache, release build, peak RSS from /usr/bin/time -l. Every row had its SA and SAindex compared with cmp between the two builders.

genome	size	libsais	caps-sa	ratio
yeast R64-1-1	12 Mb	1.0 s / 0.39 GB	1.7 s / 0.22 GB	1.6x
human chr1	249 Mb	27.4 s / 6.5 GB	47.2 s / 5.6 GB	1.7x
GRCh38 primary	3.1 Gb	641.6 s / 86.8 GB	1040.4 s / 14.4 GB	1.6x

The ratio holds at 1.6-1.7x across three orders of magnitude, so this is not a small-genome artefact. The reason it holds even though default-features = false leaves libsais single-threaded is visible in the CPU time: on GRCh38, libsais spends 726 s of user time for 641 s of wall, while caps-sa spends 6135 s of user time for 1040 s across 8 threads. caps-sa buys its small footprint with roughly six times the total CPU work for the same array.

The cost is memory: 86.8 GB against 14.4 GB, a factor of six. The 17 x 2N estimate came out 19% above the measured peak, which is the right direction for a guard.

The default declines the mammalian row. At --limitGenomeGenerateRAM 31000000000 a GRCh38 build takes the caps-sa path; the libsais number above needed an explicit --limitGenomeGenerateRAM 0. So the speedup is available to anyone up to roughly chromosome scale without touching a flag, and above that only to a user who states they have the RAM.

Reproduce from the branch, one genome per invocation:

test/bench_sa_builders.sh "$FASTA" out_bench 14 8

It drives both arms through --limitGenomeGenerateRAM, so it exercises the dispatch rather than bypassing it, reports wall clock and peak RSS for each, and compares SA and SAindex with cmp.

A caveat on the table above, which the script now guards against. These runs were taken before the script existed, on a shared workstation, and I did not record the load average at the time. That machine runs unrelated multi-core jobs, and one landing on a single arm of a minutes-long build would move that arm without leaving any trace in the numbers. The 1.6-1.7x ratio is consistent across three genomes spanning three orders of magnitude, which is not what intermittent contention usually produces, and the byte-identity checks are unaffected either way. But the timings should be re-taken on an idle machine before anyone leans on the exact figures, and the script refuses to run above a load average of 2 so the next set cannot have this problem.

Gate

568 lib + 21 integration tests, cargo clippy --all-targets -- -D warnings, cargo fmt --check, MSRV 1.89. CI green on all five platforms including windows-x86_64. On-disk index format unchanged.

This PR previously also carried --genomeType SuperTranscriptome and --genomeTransformOutput SAM; those are now #161 and #160, per the one-theme rule in CONTRIBUTING.md.