diff --git a/CHANGELOG.md b/CHANGELOG.md index b6e4893..a8b1721 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,46 @@ 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 | 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 verified against real STARsolo (#90). @@ -116,6 +156,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 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/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 a23a53d..28172aa 100644 --- a/src/index/sa_build.rs +++ b/src/index/sa_build.rs @@ -259,6 +259,222 @@ 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(()) +} + +/// 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. +pub fn libsais_peak_bytes(n_entries: u64) -> u64 { + // SA + libsais working array, 8 bytes per entry each, plus the text. + n_entries.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. +/// +/// 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_entries) <= 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 +/// 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). 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. [`libsais_fits_in_ram`] decides between the two from +/// `--limitGenomeGenerateRAM`. +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; + + 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(()) + }; + + 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); + Ok(SuffixArray { + data, + gstrand_bit, + gstrand_mask, + }) +} + +/// `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) | @@ -1048,6 +1264,115 @@ 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 libsais_ram_estimate_and_gate() { + // ~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); + + // 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)); + } + + /// 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); + } + + #[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 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 @@ -1132,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. 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 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