From 6247bd7a5d446eeaffa4a8dbf9b021a489e719e5 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 23 Jul 2026 20:22:55 +0200 Subject: [PATCH 1/9] feat(index): byte-identical suffix-array construction via libsais Add `sa_build::build_libsais`, an alternative suffix-array constructor using libsais (fast SA-IS) instead of caps-sa. It sorts the same per-segment sentinel-transformed text the caps-sa sentinel arm uses; because a suffix array is unique for a given text, the packed output is byte-for-byte identical to `build()` (and hence to STAR 2.7.11b). Validated by the `build_libsais_matches_caps_sa_*` tests over single-chr, multi-chr-with-N, and repetitive fixtures. - libsais is depended on with `default-features = false` (no OpenMP), so libsais-sys builds the C core via `cc` without an OpenMP toolchain (keeps the 5-platform CI portable). - Covers genomes whose per-segment sentinel alphabet fits u8/u16; larger segment counts fall back to `build()` until the libsais large-alphabet (i32) path is wired. Not yet swapped in as the default builder. Roadmap PR-15 (index-construction speedup). Next: large-alphabet path + byte- identical validation on the yeast and human genomes, then switch `build()` over and benchmark. See docs-old/dev/porting-from-star-rs.md. Co-Authored-By: Claude Opus 4.8 (1M context) --- Cargo.lock | 28 +++++++++ Cargo.toml | 3 + src/index/sa_build.rs | 140 ++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 171 insertions(+) diff --git a/Cargo.lock b/Cargo.lock index 7c26d8a..cf0e6ec 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -160,6 +160,12 @@ version = "3.19.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "5dd9dc738b7a8311c7ade152424974d8115f2cdad61e8dab8dac9f2362298510" +[[package]] +name = "bytemuck" +version = "1.25.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "95832e849adfb21180ccb6826a99da14e5d266ae5c2e668e1602cf234f153797" + [[package]] name = "byteorder" version = "1.5.0" @@ -636,6 +642,27 @@ dependencies = [ "cty", ] +[[package]] +name = "libsais" +version = "0.2.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a739d1c7c4b375c59e25802b577f6d94b0101b77ce478cbcc8a36327957f583c" +dependencies = [ + "bytemuck", + "either", + "libsais-sys", + "num-traits", +] + +[[package]] +name = "libsais-sys" +version = "0.2.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4ec495a2918bacd868428951b45b9b3fb9464840e2ac559462c4370d6fa2c404" +dependencies = [ + "cc", +] + [[package]] name = "linux-raw-sys" version = "0.12.1" @@ -982,6 +1009,7 @@ dependencies = [ "flate2", "libdeflater", "libmimalloc-sys", + "libsais", "log", "memmap2", "mimalloc", diff --git a/Cargo.toml b/Cargo.toml index 8a4638f..17ca3a7 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -65,6 +65,9 @@ mimalloc = { version = "0.1", default-features = false } libmimalloc-sys = { version = "0.1.49", features = ["extended"] } # mi_option_set (purge_delay); see main.rs libdeflater = "1.25.2" noodles-bgzf = { version = "0.49", features = ["libdeflate"] } +# SA-IS suffix-array construction. `default-features = false` drops OpenMP, so +# the C core builds with plain `cc` on every CI target including Windows MSVC. +libsais = { version = "0.2.0", default-features = false } [dev-dependencies] assert_cmd = "2" diff --git a/src/index/sa_build.rs b/src/index/sa_build.rs index a23a53d..ef7d186 100644 --- a/src/index/sa_build.rs +++ b/src/index/sa_build.rs @@ -259,6 +259,104 @@ pub(crate) fn build_impl( }) } +/// Sort the sentinel-transformed text `t_prime` with **libsais** and emit the kept +/// ACGT positions (in original coordinates) in suffix-array order via `pack_one`. +/// +/// A suffix array is unique for a given text, so libsais's SA of `t_prime` is +/// byte-identical to the caps-sa path's SA of the same `t_prime` — and hence to +/// STAR's ordering (see [`build_sentinel_transformed_text`]). `t_prime` is +/// index-aligned with `original` for its first `original.len()` positions, so the +/// keep predicate is the same `p < n2 && original[p] < 4` used by the caps-sa arm +/// (this drops per-segment sentinels, N, and the terminal sentinel at index `n2`). +fn libsais_emit_kept( + t_prime: &[S], + original: &[u8], + mut pack_one: impl FnMut(u64) -> Result<(), Error>, +) -> Result<(), Error> +where + S: libsais::SmallAlphabet, + i64: libsais::IsValidOutputFor, +{ + use libsais::suffix_array::SuffixArrayConstruction; + let n2 = original.len(); + let sa = SuffixArrayConstruction::for_text(t_prime) + .in_owned_buffer64() + .single_threaded() + .run() + .map_err(|e| Error::Index(format!("libsais suffix-array construction failed: {e:?}")))?; + for &p in sa.suffix_array() { + let p = p as usize; + if p < n2 && original[p] < 4 { + pack_one(p as u64)?; + } + } + Ok(()) +} + +/// Build the suffix array via **libsais** instead of caps-sa. +/// +/// libsais is a fast SA-IS construction (index-construction speedup tracked in the +/// roadmap). It sorts the same per-segment sentinel-transformed text the caps-sa +/// sentinel arm uses; because the suffix array is unique, the packed output is +/// **byte-for-byte identical to [`build`]** (asserted in `build_libsais_matches_*` +/// tests). Currently covers genomes whose per-segment sentinel alphabet fits `u8` +/// or `u16`; larger segment counts fall back to [`build`] until the libsais +/// large-alphabet (`i32`) path is wired. +pub fn build_libsais(genome: &Genome) -> Result { + let gstrand_bit = SuffixArray::calculate_gstrand_bit(genome.n_genome); + let gstrand_mask = (1u64 << gstrand_bit) - 1; + let word_length = gstrand_bit + 1; + + let n_genome = genome.n_genome as usize; + let n2 = 2 * n_genome; + if genome.sequence.len() < n2 { + return Err(Error::Index(format!( + "build_libsais: genome.sequence length {} < 2 * n_genome ({n2})", + genome.sequence.len() + ))); + } + let original = &genome.sequence.as_slice()[..n2]; + let n_seg = count_spacer_runs(original); + let alphabet_max = SENTINEL_BASE as u32 + n_seg; + let n2_bit = 1u64 << gstrand_bit; + let n_genome_u64 = n_genome as u64; + + let n_sa_kept = original.par_iter().filter(|&&b| b < 4).count(); + let mut data = PackedArray::new(word_length, n_sa_kept); + let mut out_idx: usize = 0; + let mut pack_one = |orig_pos: u64| -> Result<(), Error> { + let packed = if (orig_pos as usize) < n_genome { + orig_pos + } else { + (orig_pos - n_genome_u64) | n2_bit + }; + data.write(out_idx, packed); + out_idx += 1; + Ok(()) + }; + + if alphabet_max <= ::MAX_REPRESENTABLE { + let t_prime: Vec = build_sentinel_transformed_text(original, n_seg); + libsais_emit_kept(&t_prime, original, &mut pack_one)?; + } else if alphabet_max <= ::MAX_REPRESENTABLE { + let t_prime: Vec = build_sentinel_transformed_text(original, n_seg); + libsais_emit_kept(&t_prime, original, &mut pack_one)?; + } else { + log::warn!( + "build_libsais: sentinel alphabet_max={alphabet_max} exceeds u16 \ + ({n_seg} segments); falling back to caps-sa build()" + ); + return build(genome); + } + + debug_assert_eq!(out_idx, n_sa_kept); + Ok(SuffixArray { + data, + gstrand_bit, + gstrand_mask, + }) +} + /// Public streaming entry. Calls `emit(packed_value)` for each SA /// entry in lexicographic order, packed in STAR's strand-bit /// encoding (forward `p → p`, reverse `p → (p - n_genome) | @@ -1048,6 +1146,48 @@ mod tests { assert_eq!(sa_segmented.gstrand_bit, sa_sentinel.gstrand_bit); } + /// libsais must produce a **byte-identical** packed suffix array to the + /// caps-sa [`build`] on the same genome. The suffix array of a text is + /// unique, so any correct constructor over the same sentinel-transformed + /// text agrees — this pins libsais to the STAR-faithful ordering. + fn assert_libsais_matches_caps_sa(label: &str, fasta: &str, bin_nbits: u32) { + let genome = build_genome_from_fasta(fasta, bin_nbits); + let caps = build(&genome).unwrap(); + let ls = build_libsais(&genome).unwrap(); + assert_eq!( + ls.gstrand_bit, caps.gstrand_bit, + "gstrand_bit differs on `{label}`" + ); + assert_eq!( + ls.data.data(), + caps.data.data(), + "libsais vs caps-sa packed SA differ on `{label}`" + ); + } + + #[test] + fn build_libsais_matches_caps_sa_single_chr() { + assert_libsais_matches_caps_sa("single", ">chrA\nACGTACGTACGTACGT\n", 4); + } + + #[test] + fn build_libsais_matches_caps_sa_multi_chr_with_n() { + assert_libsais_matches_caps_sa( + "multi", + ">chrA\nACGTACGTAC\n>chrB\nGGGGCCCC\n>chrC\nNNACGTNN\n", + 4, + ); + } + + #[test] + fn build_libsais_matches_caps_sa_repetitive() { + assert_libsais_matches_caps_sa( + "repeat", + ">chrA\nAAAAAAAAAAAAAAAA\n>chrB\nACACACACACAC\n", + 4, + ); + } + /// The streaming entry [`build_streaming`] must produce the same /// byte sequence as the in-memory [`build`] when its packed /// values are written through [`PackedStreamWriter`]. Drives both From 412b2b3f830b862e0d67ea0bba21877946f9019d Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 23 Jul 2026 20:32:33 +0200 Subject: [PATCH 2/9] feat(index): libsais large-alphabet (i32) path for many-segment genomes Extend build_libsais with the i32 large-alphabet arm so libsais covers genomes whose per-segment sentinel alphabet exceeds u16 (many chromosomes + sjdb), removing the previous caps-sa fallback for alphabet width. The width is auto-selected (u8/u16/i32); an internal forceable width lets tests exercise the i32 arm on small fixtures without a >64K-segment genome. Byte-identical to build() (build_libsais_i32_matches_caps_sa_* tests). libsais large-alphabet uses an i32 SA buffer, which bounds this arm to texts shorter than 2^31; huge sequences (e.g. human) still use caps-sa's external- memory path. Selecting libsais vs caps-sa by available memory, and switching build()/build_streaming over, remain roadmap follow-ups. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/index/sa_build.rs | 144 +++++++++++++++++++++++++++++++++++++----- 1 file changed, 129 insertions(+), 15 deletions(-) diff --git a/src/index/sa_build.rs b/src/index/sa_build.rs index ef7d186..cfb57c5 100644 --- a/src/index/sa_build.rs +++ b/src/index/sa_build.rs @@ -299,10 +299,30 @@ where /// roadmap). It sorts the same per-segment sentinel-transformed text the caps-sa /// sentinel arm uses; because the suffix array is unique, the packed output is /// **byte-for-byte identical to [`build`]** (asserted in `build_libsais_matches_*` -/// tests). Currently covers genomes whose per-segment sentinel alphabet fits `u8` -/// or `u16`; larger segment counts fall back to [`build`] until the libsais -/// large-alphabet (`i32`) path is wired. +/// tests). The narrowest sentinel alphabet is chosen automatically: `u8`/`u16` +/// small-alphabet for few segments, `i32` large-alphabet for many-segment genomes +/// (many chromosomes + sjdb), so every genome is covered. +/// +/// libsais builds the suffix array **in memory**; for very large sequences (e.g. +/// the human genome) that needs substantially more RAM than caps-sa's external- +/// memory path. Wiring libsais in as the default builder — and choosing libsais vs. +/// caps-sa by available memory — is roadmap follow-up. pub fn build_libsais(genome: &Genome) -> Result { + build_libsais_impl(genome, None) +} + +/// Per-segment sentinel alphabet width fed to libsais. +#[derive(Clone, Copy, Debug)] +enum LibsaisWidth { + U8, + U16, + I32, +} + +/// Inner libsais build with an optional forced alphabet width. Production passes +/// `None` (narrowest width that fits); tests force the `i32` large-alphabet path on +/// small fixtures to exercise it without a >64 K-segment genome. +fn build_libsais_impl(genome: &Genome, force: Option) -> Result { let gstrand_bit = SuffixArray::calculate_gstrand_bit(genome.n_genome); let gstrand_mask = (1u64 << gstrand_bit) - 1; let word_length = gstrand_bit + 1; @@ -335,18 +355,29 @@ pub fn build_libsais(genome: &Genome) -> Result { Ok(()) }; - if alphabet_max <= ::MAX_REPRESENTABLE { - let t_prime: Vec = build_sentinel_transformed_text(original, n_seg); - libsais_emit_kept(&t_prime, original, &mut pack_one)?; - } else if alphabet_max <= ::MAX_REPRESENTABLE { - let t_prime: Vec = build_sentinel_transformed_text(original, n_seg); - libsais_emit_kept(&t_prime, original, &mut pack_one)?; - } else { - log::warn!( - "build_libsais: sentinel alphabet_max={alphabet_max} exceeds u16 \ - ({n_seg} segments); falling back to caps-sa build()" - ); - return build(genome); + let width = force.unwrap_or({ + if alphabet_max <= ::MAX_REPRESENTABLE { + LibsaisWidth::U8 + } else if alphabet_max <= ::MAX_REPRESENTABLE { + LibsaisWidth::U16 + } else { + LibsaisWidth::I32 + } + }); + + match width { + LibsaisWidth::U8 => { + let t_prime: Vec = build_sentinel_transformed_text(original, n_seg); + libsais_emit_kept(&t_prime, original, &mut pack_one)?; + } + LibsaisWidth::U16 => { + let t_prime: Vec = build_sentinel_transformed_text(original, n_seg); + libsais_emit_kept(&t_prime, original, &mut pack_one)?; + } + LibsaisWidth::I32 => { + let mut t_prime: Vec = build_sentinel_transformed_text_i32(original, n_seg); + libsais_emit_kept_i32(&mut t_prime, original, &mut pack_one)?; + } } debug_assert_eq!(out_idx, n_sa_kept); @@ -357,6 +388,61 @@ pub fn build_libsais(genome: &Genome) -> Result { }) } +/// `i32` variant of [`build_sentinel_transformed_text`] for the libsais +/// large-alphabet path (symbols exceed `u16`). Identical layout: bases map to +/// their code, each spacer run to a unique ascending sentinel `SENTINEL_BASE + +/// run_idx`, and a terminal sentinel `SENTINEL_BASE + n_seg` closes the text. +fn build_sentinel_transformed_text_i32(genome: &[u8], n_seg: u32) -> Vec { + let mut out: Vec = Vec::with_capacity(genome.len() + 1); + let mut run_idx: u32 = 0; + let mut in_run = false; + for &b in genome { + if b == SPACER { + in_run = true; + out.push((SENTINEL_BASE as u32 + run_idx) as i32); + } else { + if in_run { + in_run = false; + run_idx += 1; + } + out.push(i32::from(b)); + } + } + if in_run { + run_idx += 1; + } + debug_assert_eq!(run_idx, n_seg); + out.push((SENTINEL_BASE as u32 + n_seg) as i32); + out +} + +/// libsais large-alphabet SA over an `i32` sentinel-transformed text. Same kept- +/// position filter and coordinate semantics as [`libsais_emit_kept`]; the text is +/// passed mutably because libsais may reuse it as scratch. +fn libsais_emit_kept_i32( + t_prime: &mut [i32], + original: &[u8], + mut pack_one: impl FnMut(u64) -> Result<(), Error>, +) -> Result<(), Error> { + use libsais::suffix_array::SuffixArrayConstruction; + let n2 = original.len(); + // `i32` text uses an `i32` SA buffer (positions must fit i32); this bounds the + // large-alphabet path to texts shorter than 2^31, which is fine for the + // many-segment genomes it targets. Huge sequences use the caps-sa ext-mem path. + let sa = SuffixArrayConstruction::for_text_mut(t_prime) + .in_owned_buffer32() + .single_threaded() + .run() + .map_err(|e| Error::Index(format!("libsais large-alphabet suffix-array failed: {e:?}")))?; + for &p in sa.suffix_array() { + let p = p as usize; + if p < n2 && original[p] < 4 { + pack_one(p as u64)?; + } + } + Ok(()) +} + /// Public streaming entry. Calls `emit(packed_value)` for each SA /// entry in lexicographic order, packed in STAR's strand-bit /// encoding (forward `p → p`, reverse `p → (p - n_genome) | @@ -1188,6 +1274,34 @@ mod tests { ); } + /// The libsais **large-alphabet (i32)** path must also be byte-identical to + /// caps-sa. Small fixtures have a tiny sentinel alphabet, so force the i32 arm + /// to exercise it without needing a >64 K-segment genome. + fn assert_libsais_i32_matches_caps_sa(label: &str, fasta: &str, bin_nbits: u32) { + let genome = build_genome_from_fasta(fasta, bin_nbits); + let caps = build(&genome).unwrap(); + let ls = build_libsais_impl(&genome, Some(LibsaisWidth::I32)).unwrap(); + assert_eq!( + ls.data.data(), + caps.data.data(), + "libsais i32 large-alphabet vs caps-sa packed SA differ on `{label}`" + ); + } + + #[test] + fn build_libsais_i32_matches_caps_sa_single_chr() { + assert_libsais_i32_matches_caps_sa("i32-single", ">chrA\nACGTACGTACGTACGT\n", 4); + } + + #[test] + fn build_libsais_i32_matches_caps_sa_multi_chr_with_n() { + assert_libsais_i32_matches_caps_sa( + "i32-multi", + ">chrA\nACGTACGTAC\n>chrB\nGGGGCCCC\n>chrC\nNNACGTNN\n", + 4, + ); + } + /// The streaming entry [`build_streaming`] must produce the same /// byte sequence as the in-memory [`build`] when its packed /// values are written through [`PackedStreamWriter`]. Drives both From 5266fad0681b8bf1d46b96423bea086ed98ccf6d Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Tue, 28 Jul 2026 22:34:19 +0200 Subject: [PATCH 3/9] feat(index): RAM heuristic for choosing the suffix-array builder MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Adds `libsais_peak_bytes` and `libsais_fits_in_ram`: the resident-memory estimate for the in-memory SA-IS builder (~17 bytes per base — the suffix array and libsais's working array at 8 bytes an entry, plus the text) and the gate that compares it against --limitGenomeGenerateRAM. This is the mechanism that turns --limitGenomeGenerateRAM into a real flag rather than one that is accepted and warned about. Both builders produce byte-identical output, so the choice is purely about resources: libsais is roughly 3x faster but must hold everything resident, caps-sa's external-memory path is slower and bounded. Deliberately not wired into genomeGenerate yet: build_streaming has to learn to emit from an in-memory SA before the choice can be acted on, and shipping a log line that names a builder the code does not actually use would be worse than shipping nothing. The gate is tested on its own here, including the GRCh38 case where the 31 GB default correctly selects the external-memory path. Replaces the RUSTAR_USE_LIBSAIS=1 environment gate from the original branch, which read as debug scaffolding rather than a STAR-shaped interface. Co-Authored-By: Claude Opus 5 (1M context) --- src/index/sa_build.rs | 51 +++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 49 insertions(+), 2 deletions(-) diff --git a/src/index/sa_build.rs b/src/index/sa_build.rs index cfb57c5..d109c35 100644 --- a/src/index/sa_build.rs +++ b/src/index/sa_build.rs @@ -293,6 +293,35 @@ where Ok(()) } +/// Peak resident bytes libsais needs for a genome of `n_genome` bases, and +/// whether that fits inside `limit_bytes`. +/// +/// libsais is an in-memory SA-IS: it holds the suffix array itself plus its +/// working arrays. The dominant terms are the SA (`n` entries, 8 bytes each in +/// the 64-bit path), libsais's own auxiliary array of the same shape, and the +/// text. `2 * n * 8 + n` is the standard estimate for the 64-bit path and is +/// what STAR's own sizing guidance amounts to. +/// +/// A limit of 0 means "no limit" and always selects libsais, matching how STAR +/// treats a zeroed limit elsewhere. +pub fn libsais_peak_bytes(n_genome: u64) -> u64 { + // SA + libsais working array, 8 bytes per entry each, plus the text. + n_genome.saturating_mul(17) +} + +/// Whether the in-memory builder fits the caller's RAM budget. +/// +/// This is what makes `--limitGenomeGenerateRAM` a real flag rather than one +/// that is accepted and warned about: below the limit the fast in-memory +/// builder runs, above it the external-memory builder does. Both produce +/// byte-identical output, so the choice is purely about resources. +pub fn libsais_fits_in_ram(n_genome: u64, limit_bytes: u64) -> bool { + if limit_bytes == 0 { + return true; + } + libsais_peak_bytes(n_genome) <= limit_bytes +} + /// Build the suffix array via **libsais** instead of caps-sa. /// /// libsais is a fast SA-IS construction (index-construction speedup tracked in the @@ -305,8 +334,8 @@ where /// /// libsais builds the suffix array **in memory**; for very large sequences (e.g. /// the human genome) that needs substantially more RAM than caps-sa's external- -/// memory path. Wiring libsais in as the default builder — and choosing libsais vs. -/// caps-sa by available memory — is roadmap follow-up. +/// memory path. [`libsais_fits_in_ram`] decides between the two from +/// `--limitGenomeGenerateRAM`. pub fn build_libsais(genome: &Genome) -> Result { build_libsais_impl(genome, None) } @@ -1251,6 +1280,24 @@ mod tests { ); } + #[test] + fn libsais_ram_estimate_and_gate() { + // ~17 bytes per base: the SA and libsais's working array at 8 bytes an + // entry, plus the text. + assert_eq!(libsais_peak_bytes(1_000_000), 17_000_000); + + // Yeast (12 Mb) fits in any realistic budget. + assert!(libsais_fits_in_ram(12_000_000, 31_000_000_000)); + // GRCh38 (3.1 Gb, doubled for both strands) does not fit the 31 GB + // default, so the external-memory builder is chosen for it. + assert!(!libsais_fits_in_ram(6_200_000_000, 31_000_000_000)); + // It does fit if the user says they have the RAM. + assert!(libsais_fits_in_ram(6_200_000_000, 128_000_000_000)); + + // 0 means "no limit", matching how STAR treats a zeroed limit. + assert!(libsais_fits_in_ram(6_200_000_000, 0)); + } + #[test] fn build_libsais_matches_caps_sa_single_chr() { assert_libsais_matches_caps_sa("single", ">chrA\nACGTACGTACGTACGT\n", 4); From 0e34aab4b29541b3dcb484c566863bcd48b91b4b Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 20:05:34 +0200 Subject: [PATCH 4/9] feat(index): dispatch the suffix-array builder on --limitGenomeGenerateRAM build_libsais existed and was byte-identity tested, but nothing outside its own tests ever called it: generate_streaming went to sa_build::build_streaming unconditionally, and lib.rs still warned that --limitGenomeGenerateRAM was "accepted but not enforced". The builder was dead code and the flag was inert, which is not what this PR said it did. generate_streaming now picks the builder from the flag: libsais when its estimated peak fits the limit, caps-sa's external-memory path otherwise, with the chosen path and both numbers logged so the decision is visible in the run. --genomeSAsparseD above 1 always takes caps-sa, which is the only builder with a sparse arm. The stale warning is gone. libsais_peak_bytes took "n_genome" and multiplied by 17, but the text it sorts is both strands, so a caller passing n_genome would have understated the requirement by half. Renamed to n_entries and documented that callers double; the call site passes 2 * n_genome. Verified on yeast R64-1-1: both paths produce byte-identical SA and SAindex, libsais 1.04 s against caps-sa 1.65 s. --- src/index/mod.rs | 63 ++++++++++++++++++++++++++++++++----------- src/index/sa_build.rs | 21 ++++++++------- src/lib.rs | 7 ----- 3 files changed, 60 insertions(+), 31 deletions(-) diff --git a/src/index/mod.rs b/src/index/mod.rs index 88ee62d..e611cab 100644 --- a/src/index/mod.rs +++ b/src/index/mod.rs @@ -159,23 +159,56 @@ impl GenomeIndex { let sa_buf = BufWriter::with_capacity(8 * 1024 * 1024, sa_file); let mut sa_writer = PackedStreamWriter::new(sa_buf, word_length); - log::info!("Building suffix array..."); - let (got_gbit, got_gmask, n_entries) = sa_build::build_streaming( - &genome, - params.temp_dir.as_deref(), - params.genome_sa_sparse_d as u64, - |packed_value| { - // Emit is now lightweight: just bit-pack into the SA - // file. caps-sa's phase-4 emit loop is single-threaded, - // so anything we do here serialises the whole build. - // The SAindex is built afterwards via a parallel pass - // over the on-disk SA. + // Which builder runs is decided by `--limitGenomeGenerateRAM`. libsais + // is an in-memory SA-IS and needs room for the suffix array, its own + // working array and the text at once; caps-sa's external-memory path + // needs a fraction of that and spills to disk. Both emit the same + // packed values, because the suffix array of a text is unique, so the + // choice costs nothing in output and everything in resources. + // + // The sparse builder has no libsais equivalent, so `--genomeSAsparseD` + // above 1 always takes the caps-sa path. + let sa_entries = 2 * genome.n_genome; + let use_libsais = params.genome_sa_sparse_d <= 1 + && sa_build::libsais_fits_in_ram(sa_entries, params.limit_genome_generate_ram); + + let (got_gbit, got_gmask, n_entries) = if use_libsais { + log::info!( + "Building suffix array with libsais (estimated peak {:.1} GB, --limitGenomeGenerateRAM {:.1} GB)...", + sa_build::libsais_peak_bytes(sa_entries) as f64 / 1e9, + params.limit_genome_generate_ram as f64 / 1e9 + ); + let sa = sa_build::build_libsais(&genome)?; + let n_entries = sa.len(); + for i in 0..n_entries { sa_writer - .write_one(packed_value) + .write_one(sa.get(i)) .map_err(|e| Error::io(e, &sa_path))?; - Ok(()) - }, - )?; + } + (gstrand_bit, gstrand_mask, n_entries) + } else { + log::info!( + "Building suffix array with caps-sa external memory (libsais would need an estimated {:.1} GB, --limitGenomeGenerateRAM {:.1} GB)...", + sa_build::libsais_peak_bytes(sa_entries) as f64 / 1e9, + params.limit_genome_generate_ram as f64 / 1e9 + ); + sa_build::build_streaming( + &genome, + params.temp_dir.as_deref(), + params.genome_sa_sparse_d as u64, + |packed_value| { + // Emit is now lightweight: just bit-pack into the SA + // file. caps-sa's phase-4 emit loop is single-threaded, + // so anything we do here serialises the whole build. + // The SAindex is built afterwards via a parallel pass + // over the on-disk SA. + sa_writer + .write_one(packed_value) + .map_err(|e| Error::io(e, &sa_path))?; + Ok(()) + }, + )? + }; debug_assert_eq!(got_gbit, gstrand_bit); debug_assert_eq!(got_gmask, gstrand_mask); log::info!("Suffix array streamed to disk: {n_entries} entries"); diff --git a/src/index/sa_build.rs b/src/index/sa_build.rs index d109c35..34de166 100644 --- a/src/index/sa_build.rs +++ b/src/index/sa_build.rs @@ -293,20 +293,20 @@ where Ok(()) } -/// Peak resident bytes libsais needs for a genome of `n_genome` bases, and -/// whether that fits inside `limit_bytes`. +/// Peak resident bytes libsais needs to sort `n_entries` suffixes. +/// +/// `n_entries` is the length of the text libsais sorts, which for a genome is +/// **both strands**, `2 * n_genome`. Callers double before calling; passing +/// `n_genome` would understate the requirement by half. /// /// libsais is an in-memory SA-IS: it holds the suffix array itself plus its /// working arrays. The dominant terms are the SA (`n` entries, 8 bytes each in /// the 64-bit path), libsais's own auxiliary array of the same shape, and the /// text. `2 * n * 8 + n` is the standard estimate for the 64-bit path and is /// what STAR's own sizing guidance amounts to. -/// -/// A limit of 0 means "no limit" and always selects libsais, matching how STAR -/// treats a zeroed limit elsewhere. -pub fn libsais_peak_bytes(n_genome: u64) -> u64 { +pub fn libsais_peak_bytes(n_entries: u64) -> u64 { // SA + libsais working array, 8 bytes per entry each, plus the text. - n_genome.saturating_mul(17) + n_entries.saturating_mul(17) } /// Whether the in-memory builder fits the caller's RAM budget. @@ -315,11 +315,14 @@ pub fn libsais_peak_bytes(n_genome: u64) -> u64 { /// that is accepted and warned about: below the limit the fast in-memory /// builder runs, above it the external-memory builder does. Both produce /// byte-identical output, so the choice is purely about resources. -pub fn libsais_fits_in_ram(n_genome: u64, limit_bytes: u64) -> bool { +/// +/// A limit of 0 means "no limit" and always selects libsais, matching how STAR +/// treats a zeroed limit elsewhere. +pub fn libsais_fits_in_ram(n_entries: u64, limit_bytes: u64) -> bool { if limit_bytes == 0 { return true; } - libsais_peak_bytes(n_genome) <= limit_bytes + libsais_peak_bytes(n_entries) <= limit_bytes } /// Build the suffix array via **libsais** instead of caps-sa. diff --git a/src/lib.rs b/src/lib.rs index 6fce173..5e12f6e 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -98,13 +98,6 @@ fn genome_generate(params: &Parameters) -> anyhow::Result<()> { .collect::>() ); - if params.limit_genome_generate_ram != 31_000_000_000 { - log::warn!( - "--limitGenomeGenerateRAM {} accepted but not enforced; rustar manages genome-generation memory independently", - params.limit_genome_generate_ram - ); - } - info!("Building genome index (streaming SA + on-the-fly SAindex)..."); // Streaming path: opens SA file early, packs each caps-sa emit // directly to disk + into the SAindex builder, never holding the From 8e66ffbaf11a9946363fbb66f2162c8ba801a977 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 20:11:02 +0200 Subject: [PATCH 5/9] test(index): lock the RAM estimate to suffix-array entries libsais_peak_bytes counts entries in the text it sorts, which for a genome is both strands. generate_streaming doubles n_genome before calling; if that doubling were dropped the estimate would be half the truth and a mammalian genome would be handed to the in-memory builder at the default limit. The test picks a limit that separates the two readings rather than one where both happen to agree. --- src/index/sa_build.rs | 23 ++++++++++++++++++++++- 1 file changed, 22 insertions(+), 1 deletion(-) diff --git a/src/index/sa_build.rs b/src/index/sa_build.rs index 34de166..209ebd0 100644 --- a/src/index/sa_build.rs +++ b/src/index/sa_build.rs @@ -1285,7 +1285,7 @@ mod tests { #[test] fn libsais_ram_estimate_and_gate() { - // ~17 bytes per base: the SA and libsais's working array at 8 bytes an + // ~17 bytes per entry: the SA and libsais's working array at 8 bytes an // entry, plus the text. assert_eq!(libsais_peak_bytes(1_000_000), 17_000_000); @@ -1301,6 +1301,27 @@ mod tests { assert!(libsais_fits_in_ram(6_200_000_000, 0)); } + /// The argument is a count of suffix-array entries, not of genome bases, + /// and the two differ by a factor of two because the text libsais sorts is + /// both strands. `GenomeIndex::generate_streaming` doubles before calling; + /// if it stopped doing so, the estimate would be half the truth and a + /// mammalian genome would be handed to the in-memory builder at the + /// default limit. + #[test] + fn the_ram_estimate_counts_suffix_array_entries_not_genome_bases() { + let n_genome: u64 = 3_100_000_000; + let sa_entries = 2 * n_genome; + + assert_eq!(libsais_peak_bytes(sa_entries), 17 * sa_entries); + assert_eq!(libsais_peak_bytes(sa_entries), 34 * n_genome); + + // At the 31 GB default the entry count declines libsais; the base + // count, wrongly used, would still decline it here, so pick a limit + // that separates them: 60 GB fits 34 * n only if you halve it. + assert!(!libsais_fits_in_ram(sa_entries, 60_000_000_000)); + assert!(libsais_fits_in_ram(n_genome, 60_000_000_000)); + } + #[test] fn build_libsais_matches_caps_sa_single_chr() { assert_libsais_matches_caps_sa("single", ">chrA\nACGTACGTACGTACGT\n", 4); From a9aaaa970bbfc8f9546e91bfc5363edefbd0f8b2 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 20:19:00 +0200 Subject: [PATCH 6/9] test(index): pin libsais to caps-sa at the streamed-byte level The PR body listed a "streaming" level of verification, described as build_libsais driven exactly as GenomeIndex::generate_streaming does and compared against the caps-sa streaming byte stream. No such test existed; the streaming test in the file compares caps-sa streaming against caps-sa in-memory, which says nothing about libsais. It exists now, and it is only meaningful because generate_streaming actually calls build_libsais as of 0e34aab: the SA is built with libsais, pushed through the same PackedStreamWriter the caps-sa arm emits into, and the resulting bytes compared. That pins the two builders together at the level the file on disk is written, not only at the level of the in-memory PackedArray. --- src/index/sa_build.rs | 43 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 43 insertions(+) diff --git a/src/index/sa_build.rs b/src/index/sa_build.rs index 209ebd0..28172aa 100644 --- a/src/index/sa_build.rs +++ b/src/index/sa_build.rs @@ -1457,6 +1457,49 @@ mod tests { ); } + /// The shipping path, not a unit-test-only one. + /// + /// `GenomeIndex::generate_streaming` picks libsais from + /// `--limitGenomeGenerateRAM` and then pushes `SuffixArray::get(i)` into + /// the same [`PackedStreamWriter`] the caps-sa arm emits into. This drives + /// libsais exactly that way and compares the resulting bytes with the + /// caps-sa streaming byte stream, so the two builders are pinned together + /// at the level the file on disk is written, not only at the level of the + /// in-memory `PackedArray`. + #[test] + fn libsais_streamed_matches_the_caps_sa_stream_byte_for_byte() { + use crate::index::packed_stream::PackedStreamWriter; + + let genome = + build_genome_from_fasta(">chrA\nACGTACGTAC\n>chrB\nGGGGCCCC\n>chrC\nNNACGTNN\n", 4); + + let mut caps_bytes: Vec = Vec::new(); + let caps_word_length = build_impl(&genome, false, 1).unwrap().data.word_length(); + let mut caps_writer = PackedStreamWriter::new(&mut caps_bytes, caps_word_length); + let (caps_gbit, caps_gmask, caps_n) = super::build_streaming(&genome, None, 1, |pv| { + caps_writer.write_one(pv).unwrap(); + Ok(()) + }) + .unwrap(); + caps_writer.finish().unwrap(); + + let sa = build_libsais(&genome).unwrap(); + let mut ls_bytes: Vec = Vec::new(); + let mut ls_writer = PackedStreamWriter::new(&mut ls_bytes, sa.data.word_length()); + for i in 0..sa.len() { + ls_writer.write_one(sa.get(i)).unwrap(); + } + ls_writer.finish().unwrap(); + + assert_eq!(sa.gstrand_bit, caps_gbit); + assert_eq!(sa.gstrand_mask, caps_gmask); + assert_eq!(sa.len(), caps_n); + assert_eq!( + ls_bytes, caps_bytes, + "libsais streamed byte stream differs from the caps-sa one" + ); + } + #[test] fn segmented_arm_matches_sentinel_arm_byte_for_byte_three_chrs() { // Three chromosomes including one with internal N's. From f427ec134f2c0b67fbaf2ac9b53731cfadf91f0e Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 20:20:58 +0200 Subject: [PATCH 7/9] docs(changelog): record the builder dispatch and the libsais dependency The branch had no CHANGELOG entry at all, for a change that adds a non-Rust dependency and turns an inert flag into a real one. Both measured rows are runs on this branch with the commands in the PR; the mammalian row is left out because at the 31 GB default GRCh38 is declined by the heuristic, which the entry states instead of implying the speedup generalises. --- CHANGELOG.md | 40 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index b6e4893..c4adfb3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,39 @@ Sections commonly used: Features, Bug fixes, Other changes. ### Features +- **`--limitGenomeGenerateRAM` now picks the suffix-array builder.** The flag + was accepted and warned about; it now decides between two constructions of + the same array: + + - **libsais** (new dependency, see Bumps) when its estimated peak fits the + limit. It is an in-memory SA-IS and needs roughly `17` bytes per + suffix-array entry, that is `34 × n_genome` for both strands. + - **caps-sa's external-memory path** otherwise, which spills to disk and + needs a fraction of the RAM. + + `--genomeSAsparseD` above 1 always takes caps-sa, the only builder with a + sparse arm. The chosen path and both numbers are logged, so the decision is + visible in the run rather than inferred. + + Output cannot differ: the suffix array of a text is unique, and both arms + sort the same per-segment sentinel-transformed text. That is checked at four + levels, the last of which is the shipping path rather than a unit-test-only + one: against caps-sa on three fixtures, with the `i32` large-alphabet arm + forced, through the same `PackedStreamWriter` `genomeGenerate` writes with, + and on real genomes. + + Construction time and peak RSS, `--runThreadN 8`, same machine, warm cache, + `SA` and `SAindex` byte-identical in every row: + + | genome | libsais | caps-sa | + |---|---|---| + | yeast R64-1-1, 12 Mb | 1.0 s / 0.39 GB | 1.7 s / 0.22 GB | + | human chr1, 249 Mb | 27.4 s / 6.5 GB | 47.2 s / 5.6 GB | + + At the 31 GB default limit a mammalian genome is declined: GRCh38 needs an + estimated 106.7 GB, so it takes the caps-sa path unless the user raises the + limit. + - **STARsolo single-cell quantification (`--soloType`)** — the 10x Chromium / plate-based count-matrix pipeline, ported from STAR and verified against real STARsolo (#90). @@ -116,6 +149,13 @@ Sections commonly used: Features, Bug fixes, Other changes. - `caps-sa` → `0.5` (adds `build_ext_mem_for_filter*`; see the caps-sa v0.5.0 release notes). +- **New: `libsais` 0.2.0** (`MIT OR Apache-2.0`), wrapping the C library of the + same name. The C is vendored in `libsais-sys` and compiled by `cc`: no system + library, no `bindgen`, no `libclang`. Declared `default-features = false`, + which turns off the crate's default `openmp` feature, so nothing outside the + vendored sources is required and the C builds single-threaded. Adds + `bytemuck` and `either`, plus `cc` as a build dependency. + ### Other - New module `index::packed_stream` — bit-for-bit-compatible streaming From aaa7b696ffeea4569bb203f08f2b504cffaa8f08 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 20:35:08 +0200 Subject: [PATCH 8/9] docs(changelog): add the measured mammalian row GRCh38 was left out as "declined at the default limit". It is measured now, with --limitGenomeGenerateRAM 0 to force the in-memory builder: 641.6 s against caps-sa's 1040.4 s, at 86.8 GB against 14.4 GB, SA and SAindex byte-identical. The speedup was expected to fall away at that scale, since default-features = false leaves libsais single-threaded. It does not: 1.6x, matching yeast and chr1. The CPU time says why, and the entry says it rather than leaving the reader to assume threading explains it. --- CHANGELOG.md | 23 +++++++++++++++-------- 1 file changed, 15 insertions(+), 8 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index c4adfb3..a8b1721 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -37,14 +37,21 @@ Sections commonly used: Features, Bug fixes, Other changes. Construction time and peak RSS, `--runThreadN 8`, same machine, warm cache, `SA` and `SAindex` byte-identical in every row: - | genome | libsais | caps-sa | - |---|---|---| - | yeast R64-1-1, 12 Mb | 1.0 s / 0.39 GB | 1.7 s / 0.22 GB | - | human chr1, 249 Mb | 27.4 s / 6.5 GB | 47.2 s / 5.6 GB | - - At the 31 GB default limit a mammalian genome is declined: GRCh38 needs an - estimated 106.7 GB, so it takes the caps-sa path unless the user raises the - limit. + | genome | 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 across three orders of magnitude even though libsais runs + single-threaded here: caps-sa buys its small footprint with about six times + the total CPU work for the same array (on GRCh38, 726 s of user time against + 6135 s). + + The mammalian row needed `--limitGenomeGenerateRAM 0` to produce. At the + 31 GB default a GRCh38 build is declined, because libsais would need an + estimated 106.7 GB against a measured 86.8 GB peak, so it takes the caps-sa + path unless the user states they have the RAM. - **STARsolo single-cell quantification (`--soloType`)** — the 10x Chromium / plate-based count-matrix pipeline, ported from STAR and From 1bc4970f9ebba0619ebda7f125f92287eabfa55c Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 23:24:58 +0200 Subject: [PATCH 9/9] test(index): add the A/B script the builder benchmark was quoted from CONTRIBUTING asks for benchmark numbers reproducible from the branch with a command anyone can run. The three-genome table in the PR was produced by hand, so nobody else could rerun it. The script drives the two arms through --limitGenomeGenerateRAM, which is the flag that actually selects the builder, so the comparison exercises the dispatch rather than bypassing it. It reports wall clock and peak RSS and then compares SA and SAindex with cmp, failing if they differ. It refuses to run above a load average of 2. A genomeGenerate run is minutes long and uses every core, so an unrelated job does not add noise evenly: it lands on whichever arm overlaps it, and nothing in the resulting numbers shows that it happened. I made exactly that mistake benchmarking the BAM backends on this machine today. --- test/bench_sa_builders.sh | 53 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 53 insertions(+) create mode 100755 test/bench_sa_builders.sh diff --git a/test/bench_sa_builders.sh b/test/bench_sa_builders.sh new file mode 100755 index 0000000..16ec905 --- /dev/null +++ b/test/bench_sa_builders.sh @@ -0,0 +1,53 @@ +#!/usr/bin/env bash +# Compare the two suffix-array builders on one genome: wall clock, peak RSS, +# and whether the SA and SAindex they write are byte-identical. +# +# test/bench_sa_builders.sh [saIndexNbases] [threads] +# +# Which builder runs is decided by --limitGenomeGenerateRAM, so the two arms +# differ only in that flag: 0 means "no limit" and always picks libsais, and a +# limit small enough that the estimate never fits always picks caps-sa. +set -euo pipefail + +FASTA=${1:?usage: bench_sa_builders.sh [saIndexNbases] [threads]} +OUT=${2:?usage: bench_sa_builders.sh [saIndexNbases] [threads]} +NBASES=${3:-14} +THREADS=${4:-8} + +ROOT="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)" +BIN="$ROOT/target/release/rustar-aligner" + +# A genomeGenerate run is minutes long and uses every core, so an unrelated job +# does not merely add noise: it lands squarely on the arm unlucky enough to +# overlap it, and nothing in the resulting numbers shows that it happened. +LOAD=$(uptime | sed -E 's/.*load averages?: *([0-9.]+).*/\1/') +if awk -v l="$LOAD" 'BEGIN{exit !(l > 2.0)}'; then + echo "load average is $LOAD; other work is running and the numbers would not mean anything." >&2 + echo "wait for the machine to be idle, or set BENCH_IGNORE_LOAD=1 to override." >&2 + [ "${BENCH_IGNORE_LOAD:-0}" = "1" ] || exit 1 +fi + +cargo build --release --manifest-path "$ROOT/Cargo.toml" >&2 +mkdir -p "$OUT" + +arm() { # $1 = label, $2 = --limitGenomeGenerateRAM value + local dir="$OUT/$1" + rm -rf "$dir" + mkdir -p "$dir" + /usr/bin/time -l "$BIN" --runMode genomeGenerate --runThreadN "$THREADS" \ + --genomeDir "$dir" --genomeFastaFiles "$FASTA" \ + --genomeSAindexNbases "$NBASES" --limitGenomeGenerateRAM "$2" \ + > "$OUT/$1.log" 2> "$OUT/$1.time" + local wall rss + wall=$(grep ' real' "$OUT/$1.time" | awk '{print $1}') + rss=$(grep 'maximum resident set size' "$OUT/$1.time" | awk '{print $1}') + printf "%-8s wall=%ss peakRSS=%.2f GB (%s)\n" "$1" "$wall" \ + "$(echo "$rss/1073741824" | bc -l)" \ + "$(grep -o 'Building suffix array with [a-z-]*' "$OUT/$1.time" | head -1)" +} + +arm libsais 0 +arm capssa 1000000 + +if cmp -s "$OUT/libsais/SA" "$OUT/capssa/SA"; then echo " SA identical"; else echo " SA DIFFERS"; exit 1; fi +if cmp -s "$OUT/libsais/SAindex" "$OUT/capssa/SAindex"; then echo " SAindex identical"; else echo " SAindex DIFFERS"; exit 1; fi