From be615aa2ad693605733e05be162fb685bd5e2b7b Mon Sep 17 00:00:00 2001 From: alejandrogzi Date: Fri, 11 Sep 2026 10:21:29 +0200 Subject: [PATCH 1/3] =?UTF-8?q?feat:=20sync=20experimental=20rounds=2092?= =?UTF-8?q?=20and=2094=20=E2=80=94=20hspz=20hits-estimate=20and=20--query-?= =?UTF-8?q?list=20batch=20(experimental=2029e9e436fef97a39aa67909368d062bb?= =?UTF-8?q?679bd1f9)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/cli.rs | 108 +++- src/compare.rs | 3 +- src/hits_estimate.rs | 660 +++++++++++++++++++ src/main.rs | 2 + src/plan.rs | 29 +- src/run.rs | 1456 ++++++++++++++++++++++++++++++++++++++---- src/seed.rs | 310 +++++++++ 7 files changed, 2434 insertions(+), 134 deletions(-) create mode 100644 src/hits_estimate.rs diff --git a/src/cli.rs b/src/cli.rs index 7741eae..dd87419 100644 --- a/src/cli.rs +++ b/src/cli.rs @@ -36,14 +36,27 @@ pub(crate) enum Command { /// Run the C++ CUDA reference and this implementation back to back and /// compare runtime and HSP output. Compare(CompareArgs), + /// Exact per-unit seed-hit estimates before any GPU work (CPU-only). + HitsEstimate(HitsEstimateArgs), } #[derive(Args, Clone)] pub(crate) struct RunArgs { #[arg(short, long)] pub(crate) reference: PathBuf, - #[arg(short, long)] - pub(crate) query: PathBuf, + #[arg( + short, + long, + required_unless_present = "query_list", + conflicts_with = "query_list" + )] + pub(crate) query: Option, + /// One query FASTA path per line (blank lines and `#` comments ignored). + /// Bin-major batch: one reference × many queries, W=1 only. All-or-nothing: + /// a failing job aborts the batch, so partial per-job outputs (and a + /// truncated `-Z` archive) may remain on disk. + #[arg(long, conflicts_with = "query")] + pub(crate) query_list: Option, /// Directory to write `tmp.block.r.{plus,minus}.segments` into. #[arg(short, long, default_value = ".")] pub(crate) output: PathBuf, @@ -189,6 +202,9 @@ pub(crate) struct RunArgs { pub(crate) diagonal_partition: bool, /// Write output into one `.tar.gz` instead of a directory. Without a path, /// `.tar.gz` is used. Archived directly from the formatted bytes. + /// In batch mode (`--query-list`) one archive per job is written + /// (`OUT/000001.tar.gz` …) and a custom path value is ignored (a note is + /// printed). #[arg(short = 'Z', long, num_args = 0..=1, default_missing_value = "-", require_equals = false)] pub(crate) tarball: Option, } @@ -252,6 +268,53 @@ pub(crate) struct CompareArgs { pub(crate) tuning: Tuning, } +/// CPU-only exact seed-hit estimator (R92 PR1): same planning inputs as +/// `run` (bins/blocks identical to what `run` would freeze), no CUDA. +#[derive(Args, Clone)] +pub(crate) struct HitsEstimateArgs { + #[arg(short, long)] + pub(crate) reference: PathBuf, + #[arg(short, long)] + pub(crate) query: PathBuf, + /// Target bin size in bases, as in `run` (`-B 0` resolves to the default). + #[arg(short = 'B', long, default_value_t = 500_000_000)] + pub(crate) seq_block_size: u32, + /// Bin target for the query side; defaults to `--seq-block-size`. + #[arg(long)] + pub(crate) query_block_size: Option, + /// Bin records the way KegAlign does (sequential fill, not LPT). + #[arg(long)] + pub(crate) kegalign_bins: bool, + /// 12of19, 14of22, or an arbitrary pattern of 1s, 0s and Ts (k<=12). + #[arg(long, default_value = "12of19")] + pub(crate) seed: String, + #[arg(short = 'e', long, default_value_t = 1)] + pub(crate) step: u32, + /// Don't allow one transition in a seed hit. + #[arg(long)] + pub(crate) notransition: bool, + /// plus / minus / both + #[arg(short = 'S', long, default_value = "both")] + pub(crate) strand: String, + /// Worker threads for histogram passes; 0 uses available parallelism. + #[arg(short, long, default_value_t = 0)] + pub(crate) threads: usize, + /// LASTZ interval size, as in `run` (`-I`): interval boundaries decide + /// which shared window starts the pipeline seeds once and which twice. + #[arg(short = 'I', long, default_value_t = 10_000_000)] + pub(crate) lastz_interval_size: u32, + /// Query seed chunk size, as in `run` (`-C`): the interval schedule + /// (`seed::chunks`) is part of the count, so this must match `run`. + #[arg(short = 'C', long, default_value_t = 250_000)] + pub(crate) wga_chunk_size: u32, + /// Count every S-th query window start and scale by S (approximation). + #[arg(long, default_value_t = 1)] + pub(crate) stride: u32, + /// Write the plan's bin membership to this file and continue. + #[arg(long)] + pub(crate) dump_plan: Option, +} + #[derive(Args, Clone)] pub(crate) struct Tuning { #[arg(short, long, default_value = "12of19")] @@ -279,7 +342,7 @@ mod tests { { Command::Run(args) => { assert_eq!(args.reference, PathBuf::from("r.fa")); - assert_eq!(args.query, PathBuf::from("q.fa")); + assert_eq!(args.query, Some(PathBuf::from("q.fa"))); assert_eq!(args.output, PathBuf::from("out")); } _ => panic!("wrong subcommand"), @@ -312,5 +375,44 @@ mod tests { Command::Run(args) => assert_eq!(args.tarball, Some(PathBuf::from("o.tgz"))), _ => panic!("wrong subcommand"), } + // `--query-list` is mutually exclusive with `-q/--query`, and `-q` + // becomes optional so one of the two must be present. + assert!( + Cli::try_parse_from(["hspz", "run", "-r", "r.fa", "-o", "out"]).is_err(), + "one of -q/--query-list is required" + ); + match Cli::try_parse_from([ + "hspz", + "run", + "-r", + "r.fa", + "--query-list", + "q.txt", + "-o", + "out", + ]) + .unwrap() + .command + { + Command::Run(args) => { + assert_eq!(args.query, None); + assert_eq!(args.query_list, Some(PathBuf::from("q.txt"))); + } + _ => panic!("wrong subcommand"), + } + assert!( + Cli::try_parse_from([ + "hspz", + "run", + "-r", + "r.fa", + "-q", + "q.fa", + "--query-list", + "q.txt", + ]) + .is_err(), + "-q and --query-list must conflict" + ); } } diff --git a/src/compare.rs b/src/compare.rs index a0e78ac..d00cf99 100644 --- a/src/compare.rs +++ b/src/compare.rs @@ -30,7 +30,8 @@ pub(crate) fn compare(args: &CompareArgs) -> Fallible<()> { eprintln!("running Rust implementation ..."); let run_args = RunArgs { reference: args.reference.clone(), - query: args.query.clone(), + query: Some(args.query.clone()), + query_list: None, output: rust_dir.clone(), strand: "both".into(), seed: args.tuning.seed.clone(), diff --git a/src/hits_estimate.rs b/src/hits_estimate.rs new file mode 100644 index 0000000..e547f59 --- /dev/null +++ b/src/hits_estimate.rs @@ -0,0 +1,660 @@ +// Copyright (c) 2026 The Hiller Lab at the Senckenberg Gesellschaft für Naturforschung +// Distributed under the terms of the GNU General Public License, Version 3.0. + +// Author : Alejandro Gonzales-Irribarren +// Github : alejandrogzi +// Email : alejandrxgzi@gmail.com + +//! The `hits-estimate` command (R92 PR1): exact per-unit seed-hit counts, +//! computed on the host before any GPU work. +//! +//! CPU-only: no CUDA calls anywhere on this path, so it runs on a node +//! without a GPU. Planning reuses `plan::plan_with` with the same flags and +//! defaults as `run`, so bins/blocks are identical to what `run` would freeze +//! whenever `run` does not shrink the plan against a device budget (all PR1 +//! tests and the default whole-genome layout on a fitting GPU; see report.md +//! for the documented deviations). + +use crate::Fallible; +use crate::cli::HitsEstimateArgs; +use crate::plan::{self, PackedBin}; +use crate::seed::{self, Shape}; +use crate::sequence; +use crate::timing; +use std::io::Write; +use std::time::Instant; + +/// What one `estimate` call produced, for the stderr line. +pub(crate) struct EstimateReport { + pub(crate) r: usize, + pub(crate) q: usize, + pub(crate) threads: usize, + pub(crate) stride: u32, + pub(crate) load_ms: f64, + pub(crate) ref_hist_ms: f64, + pub(crate) qry_hist_ms: f64, + pub(crate) fold_dot_ms: f64, + pub(crate) total_ms: f64, + pub(crate) peak_rss_mib: u64, +} + +/// Loads the inputs, plans like `run`, histograms every bin and dots to +/// row-major `r*Q+q` predictions (`--stride` scaled). +pub(crate) fn estimate( + args: &HitsEstimateArgs, +) -> Fallible<(plan::Plan, Vec, EstimateReport)> { + let total = Instant::now(); + if args.stride == 0 { + return Err("--stride must be >= 1".into()); + } + if args.lastz_interval_size == 0 { + return Err("--lastz-interval-size must be >= 1".into()); + } + if args.wga_chunk_size == 0 { + return Err("--wga-chunk-size must be >= 1".into()); + } + let plus = args.strand == "plus" || args.strand == "both"; + let minus = args.strand == "minus" || args.strand == "both"; + if !plus && !minus { + return Err(format!("--strand must be plus, minus or both, got {}", args.strand).into()); + } + let shape = Shape::parse(&args.seed)?; + let transitions = !args.notransition; + let threads = crate::run::resolve_threads(args.threads); + + let t = Instant::now(); + let (_, ref_records, _) = sequence::read_records(&args.reference)?; + let (_, qry_records, _) = sequence::read_records(&args.query)?; + let load_ms = t.elapsed().as_secs_f64() * 1000.0; + + if args.seq_block_size == 0 && args.kegalign_bins { + return Err(plan::AUTO_KEGALIGN_ERROR.into()); + } + // Same resolution as `run` at one worker: `-B 0` means the default target. + let ref_target = if args.seq_block_size == 0 { + plan::DEFAULT_BLOCK_TARGET + } else { + u64::from(args.seq_block_size) + }; + let qry_target = args.query_block_size.map(u64::from).unwrap_or(ref_target); + let plan = plan::plan_with( + &crate::run::record_meta(&ref_records), + &crate::run::record_meta(&qry_records), + ref_target.max(1), + qry_target.max(1), + args.kegalign_bins, + ); + + if let Some(path) = &args.dump_plan { + let mut f = std::io::BufWriter::new(std::fs::File::create(path)?); + for (side, bins, recs) in [ + ("reference", &plan.reference_bins, &ref_records), + ("query", &plan.query_bins, &qry_records), + ] { + for b in bins { + for &id in &b.record_ids { + let (name, seq) = &recs[id as usize]; + writeln!(f, "{side}\t{}\t{name}\t{}", b.id, seq.len())?; + } + } + } + f.flush()?; + } + + // Histograms over exactly what the executor consumes: packed bins from + // the shared `PackedBin::build` (byte-identical to `run`). + let t = Instant::now(); + let mut ref_hists = Vec::with_capacity(plan.reference_bins.len()); + for b in &plan.reference_bins { + let packed = PackedBin::build( + b.record_ids.iter().map(|&id| { + let (n, s) = &ref_records[id as usize]; + (n.as_str(), s.as_slice()) + }), + "", + false, + ); + ref_hists.push(seed::count_ref_block_parallel( + &packed.buf[..packed.block_len], + &shape, + args.step, + threads, + )?); + } + let ref_hist_ms = t.elapsed().as_secs_f64() * 1000.0; + + let t = Instant::now(); + let mut qry_hists = Vec::with_capacity(plan.query_bins.len()); + for b in &plan.query_bins { + let packed = PackedBin::build( + b.record_ids.iter().map(|&id| { + let (n, s) = &qry_records[id as usize]; + (n.as_str(), s.as_slice()) + }), + "", + true, + ); + qry_hists.push(seed::count_query_block( + &packed.buf[..packed.block_len], + &packed.rc, + packed.block_len, + &shape, + plus, + minus, + args.lastz_interval_size, + args.wga_chunk_size, + args.stride, + threads, + )?); + } + let qry_hist_ms = t.elapsed().as_secs_f64() * 1000.0; + + let t = Instant::now(); + let hits = seed::unit_hits(&ref_hists, &qry_hists, &shape, transitions); + let fold_dot_ms = t.elapsed().as_secs_f64() * 1000.0; + + // `--stride` S counts every S-th query window start: scale back up. + let scale = u64::from(args.stride); + let mut pred = Vec::with_capacity(hits.len()); + for &h in &hits { + pred.push( + h.checked_mul(scale) + .ok_or("hits-estimate: hit count overflow")?, + ); + } + let report = EstimateReport { + r: plan.reference_bins.len(), + q: plan.query_bins.len(), + threads, + stride: args.stride, + load_ms, + ref_hist_ms, + qry_hist_ms, + fold_dot_ms, + total_ms: total.elapsed().as_secs_f64() * 1000.0, + peak_rss_mib: timing::peak_rss_kib() / 1024, + }; + Ok((plan, pred, report)) +} + +/// TSV rows to stdout, one stderr summary line (sufficient for whole-genome +/// timing on the 24-thread node). +pub(crate) fn run(args: &HitsEstimateArgs) -> Fallible<()> { + let (plan, pred, rep) = estimate(args)?; + let q = rep.q; + let mut out = std::io::BufWriter::new(std::io::stdout().lock()); + for (r, rb) in plan.reference_bins.iter().enumerate() { + for (qq, qb) in plan.query_bins.iter().enumerate() { + writeln!( + out, + "{r}\t{qq}\t{}\t{}\t{}", + rb.total_bp, + qb.total_bp, + pred[r * q + qq] + )?; + } + } + out.flush()?; + eprintln!( + "hits-estimate: R={} Q={} threads={} stride={} load_ms={:.0} ref_hist_ms={:.0} \ + qry_hist_ms={:.0} fold_dot_ms={:.0} total_ms={:.0} peak_rss_mib={}", + rep.r, + rep.q, + rep.threads, + rep.stride, + rep.load_ms, + rep.ref_hist_ms, + rep.qry_hist_ms, + rep.fold_dot_ms, + rep.total_ms, + rep.peak_rss_mib, + ); + Ok(()) +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::cli::{Cli, Command}; + use crate::seed::{SeedTable, chunk_seeds}; + use clap::Parser; + + /// Independent oracle for one packed (ref, query) pair: the existing + /// `cpu_stats` logic — `SeedTable::hit_count` summed over `chunk_seeds` + /// on both strands — which the estimator must reproduce exactly. The + /// window starts come from [`seed::chunks`], the same batch walk the + /// device seeder and `cpu_stats` both execute. + fn oracle_hits( + ref_buf: &[u8], + fwd: &[u8], + rc: &[u8], + block_len: usize, + shape: &Shape, + step: u32, + transitions: bool, + plus: bool, + minus: bool, + lastz: u32, + wga_chunk: u32, + ) -> u64 { + let table = SeedTable::build(ref_buf, shape, step); + let mut hits = 0u64; + let qbl = (block_len - shape.size) as u32; + for &(s, e) in &sequence::intervals(block_len, shape.size, lastz) { + if plus { + for (lo, hi) in seed::chunks(s, e, wga_chunk) { + for seed in chunk_seeds(fwd, shape, transitions, (lo, hi)) { + hits += table.hit_count((seed >> 32) as u32) as u64; + } + } + } + if minus { + for (lo, hi) in seed::chunks(qbl - e, qbl - s, wga_chunk) { + for seed in chunk_seeds(rc, shape, transitions, (lo, hi)) { + hits += table.hit_count((seed >> 32) as u32) as u64; + } + } + } + } + hits + } + + /// Hand-counted toy: short ACGT records with one N run and one lowercase + /// run per side, both strands, transitions on and off, steps 1-3, against + /// the `SeedTable` + `chunk_seeds` oracle. Runs in the normal suite. + #[test] + fn toy_estimator_matches_pipeline_oracle() { + let shape = Shape::parse("TTT0T").unwrap(); // size 5, k 4, all-transition + let lastz = 7; // force several intervals on ~20 bp blocks + let ref_recs = vec![ + ("r1", b"ACGTACGTACGTACGT".as_slice()), + ("r2", b"AAAANNNNAAAACCCC".as_slice()), + ("r3", b"acgtACGTacgtACGT".as_slice()), + ]; + let qry_recs = vec![ + ("q1", b"TTTTGGGGAAAACCCCGGGG".as_slice()), + ("q2", b"ACnnACGTACGTACGTAC".as_slice()), + ]; + let (ref_buf, _, ref_len) = sequence::pack(ref_recs.iter().map(|(n, s)| (*n, *s)), ""); + let (fwd, fwd_chrs, block_len) = sequence::pack(qry_recs.iter().map(|(n, s)| (*n, *s)), ""); + let (rc, _) = sequence::reverse_complement(&fwd, &fwd_chrs, block_len); + assert!(ref_buf.iter().any(|&c| c == b'N')); + assert!(ref_buf.iter().any(|&c| c == b'a')); + + for transitions in [true, false] { + for (plus, minus) in [(true, true), (true, false), (false, true)] { + for step in [1u32, 2, 3] { + let want = oracle_hits( + &ref_buf[..ref_len], + &fwd[..block_len], + &rc, + block_len, + &shape, + step, + transitions, + plus, + minus, + lastz, + 250_000, + ); + for threads in [1usize, 3] { + let rh = seed::count_ref_block_parallel( + &ref_buf[..ref_len], + &shape, + step, + threads, + ) + .unwrap(); + // Serial and parallel reference counts agree. + assert_eq!( + rh, + seed::count_ref_block(&ref_buf[..ref_len], &shape, step).unwrap() + ); + let qh = seed::count_query_block( + &fwd[..block_len], + &rc, + block_len, + &shape, + plus, + minus, + lastz, + 250_000, + 1, + threads, + ) + .unwrap(); + let got = seed::unit_hits(&[rh], &[qh], &shape, transitions); + assert_eq!( + got, + vec![want], + "transitions={transitions} plus={plus} minus={minus} step={step} threads={threads}" + ); + } + } + } + } + } + + /// What `seed::chunks` does at a shared endpoint is the entire boundary + /// question. Consecutive intervals share their endpoint (`intervals()` + /// emits `(0,7),(7,11),...`); the left interval seeds it only when its + /// chunk walk does not land exactly on it — i.e. when its length is not a + /// multiple of `wga_chunk` — otherwise the right interval seeds it, once. + /// With a chunk larger than the interval (the test's 250,000 against 7) + /// every endpoint lands inside the single chunk and is seeded twice. + /// Packed-bin counts therefore do NOT equal summed per-record counts in + /// general — the boundaries differ — which is why the estimator histograms + /// packed bins, not records. + #[test] + fn interval_boundary_windows_are_counted_twice_like_the_pipeline() { + let shape = Shape::parse("TTT0T").unwrap(); + // All-valid 16 bp record, lastz 7 -> intervals (0,7),(7,11), default + // chunk: start 7 is seeded in both intervals, 13 seeds from 12 windows. + let seq = b"ACGTACGTACGTACGT".as_slice(); + let bin = PackedBin::build([("r1", seq)], "", true); + let h = seed::count_query_block( + &bin.buf[..bin.block_len], + &bin.rc, + bin.block_len, + &shape, + true, + false, + 7, + 250_000, + 1, + 1, + ) + .unwrap(); + let total: u32 = h.iter().sum(); + assert_eq!(total, 13, "boundary start 7 must be double-counted"); + // And the oracle agrees (it walks the same chunks). + let want = oracle_hits( + &bin.buf[..bin.block_len], + &bin.buf[..bin.block_len], + &bin.rc, + bin.block_len, + &shape, + 1, + true, + true, + false, + 7, + 250_000, + ); + let got = seed::unit_hits( + &[seed::count_ref_block(&bin.buf[..bin.block_len], &shape, 1).unwrap()], + &[h], + &shape, + true, + ); + assert_eq!(got, vec![want]); + } + + /// Multi-interval exactness on a 200 bp record, CPU-only. With `-I 40` + /// and `-C 8` (`40 % 8 == 0`) the shared endpoints are seeded once and the + /// estimator must NOT double them; with the default `-C 250000` the single + /// chunk spans each interval and they are seeded twice. `-I 45` is not a + /// multiple of 8, so both chunk sizes seed every endpoint twice. The oracle + /// is the `cpu_stats`/device batch walk (`seed::chunks` + `chunk_seeds` + + /// `hit_count`); the chunk sizes are also shown to disagree at `-I 40`. + #[test] + fn multi_interval_endpoints_follow_wga_chunking() { + let shape = Shape::parse("TTT0T").unwrap(); + let seq: Vec = (0..200u32) + .map(|i| b"ACGT"[(i.wrapping_mul(7) % 4) as usize]) + .collect(); + let bin = PackedBin::build([("q1", seq.as_slice())], "", true); + let mut totals = Vec::new(); + for &(lastz, chunk) in &[(40u32, 8u32), (40, 250_000), (45, 8), (45, 250_000)] { + let want = oracle_hits( + &bin.buf[..bin.block_len], + &bin.buf[..bin.block_len], + &bin.rc, + bin.block_len, + &shape, + 1, + true, + true, + true, + lastz, + chunk, + ); + totals.push(want); + for threads in [1usize, 3] { + let rh = + seed::count_ref_block_parallel(&bin.buf[..bin.block_len], &shape, 1, threads) + .unwrap(); + let qh = seed::count_query_block( + &bin.buf[..bin.block_len], + &bin.rc, + bin.block_len, + &shape, + true, + true, + lastz, + chunk, + 1, + threads, + ) + .unwrap(); + let got = seed::unit_hits(&[rh], &[qh], &shape, true); + assert_eq!( + got, + vec![want], + "lastz={lastz} chunk={chunk} threads={threads}" + ); + } + } + assert!( + totals[0] < totals[1], + "C=8 must drop the duplicated endpoints: {totals:?}" + ); + assert_eq!( + totals[2], totals[3], + "I=45 endpoint rule is chunk-independent: {totals:?}" + ); + } + + #[test] + fn wide_seeds_error_cleanly() { + let shape = Shape::parse("14of22").unwrap(); + assert!(seed::count_ref_block(b"ACGTACGT", &shape, 1).is_err()); + assert!( + seed::count_query_block( + b"ACGTACGT", + b"ACGTACGT", + 8, + &shape, + true, + true, + 7, + 250_000, + 1, + 1 + ) + .is_err() + ); + } + + #[test] + fn hits_estimate_flags_parse() { + match Cli::try_parse_from(["hspz", "hits-estimate", "-r", "r.fa", "-q", "q.fa"]) + .unwrap() + .command + { + Command::HitsEstimate(a) => { + assert_eq!(a.seq_block_size, 500_000_000); + assert_eq!(a.query_block_size, None); + assert!(!a.kegalign_bins); + assert_eq!(a.seed, "12of19"); + assert_eq!(a.step, 1); + assert!(!a.notransition); + assert_eq!(a.strand, "both"); + assert_eq!(a.threads, 0); + assert_eq!(a.lastz_interval_size, 10_000_000); + assert_eq!(a.wga_chunk_size, 250_000); + assert_eq!(a.stride, 1); + assert_eq!(a.dump_plan, None); + } + _ => panic!("wrong subcommand"), + } + // `-I`/`-C` must be real flags with the same names as `run`. + match Cli::try_parse_from([ + "hspz", + "hits-estimate", + "-r", + "r.fa", + "-q", + "q.fa", + "-I", + "1000000", + "-C", + "8", + ]) + .unwrap() + .command + { + Command::HitsEstimate(a) => { + assert_eq!(a.lastz_interval_size, 1_000_000); + assert_eq!(a.wga_chunk_size, 8); + } + _ => panic!("wrong subcommand"), + } + } + + /// A zero interval size or chunk size would loop forever in + /// `sequence::intervals` / `seed::chunks`; the estimator must reject both + /// before touching any input. + #[test] + fn zero_interval_or_chunk_errors() { + let mut args = chr20_args(); + args.lastz_interval_size = 0; + assert!(estimate(&args).is_err()); + args.lastz_interval_size = 10_000_000; + args.wga_chunk_size = 0; + assert!(estimate(&args).is_err()); + } + + fn chr20_args() -> HitsEstimateArgs { + HitsEstimateArgs { + reference: "/tmp/hspz-cycle3-chr20/ref.fa".into(), + query: "/tmp/hspz-cycle3-chr20/qry.fa".into(), + seq_block_size: 500_000_000, + query_block_size: None, + kegalign_bins: false, + seed: "12of19".into(), + step: 1, + notransition: false, + strand: "both".into(), + threads: 0, + lastz_interval_size: 10_000_000, + wga_chunk_size: 250_000, + stride: 1, + dump_plan: None, + } + } + + /// chr20 pair, default flags: the single predicted unit equals `#seed + /// hits` = 192,899,566 (round-90 ledger; re-derived with + /// `hspz run --time --gpus 1` — see report.md). Needs the input files; + /// the GPU re-derivation itself is manual. + #[test] + #[ignore] + fn chr20_single_unit_matches_seed_hits() { + let (plan, pred, _) = estimate(&chr20_args()).unwrap(); + assert_eq!((plan.reference_bins.len(), plan.query_bins.len()), (1, 1)); + assert_eq!(pred, vec![192_899_566]); + } + + /// chr20, multi-interval: with `-I 1000000` (4 intervals) the pipeline + /// reports `#seed hits` = 192,899,602 in all three paths — device seeder, + /// CPU seeder and `--cpu-only` (report2.md) — which the estimator with the + /// same `-I`/`-C` must reproduce exactly. The shared endpoints are seeded + /// once because 1,000,000 is a multiple of the default chunk 250,000. + #[test] + #[ignore] + fn chr20_interval_1000000_matches_seed_hits() { + let args = HitsEstimateArgs { + lastz_interval_size: 1_000_000, + ..chr20_args() + }; + let (plan, pred, _) = estimate(&args).unwrap(); + assert_eq!((plan.reference_bins.len(), plan.query_bins.len()), (1, 1)); + assert_eq!(pred, vec![192_899_602]); + } + + /// chr20, `-I 300000`: the interval length is not a multiple of the + /// default chunk, so both seeders double-seed every shared endpoint and + /// `#seed hits` = 192,900,034 (report2.md). The estimator must agree. + #[test] + #[ignore] + fn chr20_interval_300000_matches_seed_hits() { + let args = HitsEstimateArgs { + lastz_interval_size: 300_000, + ..chr20_args() + }; + let (plan, pred, _) = estimate(&args).unwrap(); + assert_eq!((plan.reference_bins.len(), plan.query_bins.len()), (1, 1)); + assert_eq!(pred, vec![192_900_034]); + } + + /// The 5x3 synthetic plan (`-B 10000000 --query-block-size 1000000`): + /// all 15 predictions equal the ledger's 15 `hits` values (from + /// `hspz run --time --gpus 2 --max-hits 300000` on ZLUDA — see + /// report.md). Row-major `r*3+q`. + #[test] + #[ignore] + fn synthetic_5x3_plan_matches_ledger_hits() { + let args = HitsEstimateArgs { + reference: "/tmp/opencode/r90/ref5.fa".into(), + query: "/tmp/opencode/r90/qry3.fa".into(), + seq_block_size: 10_000_000, + query_block_size: Some(1_000_000), + kegalign_bins: false, + seed: "12of19".into(), + step: 1, + notransition: false, + strand: "both".into(), + threads: 0, + lastz_interval_size: 10_000_000, + wga_chunk_size: 250_000, + stride: 1, + dump_plan: None, + }; + let (plan, pred, _) = estimate(&args).unwrap(); + assert_eq!((plan.reference_bins.len(), plan.query_bins.len()), (5, 3)); + assert_eq!( + pred, + vec![ + 14_652_192, 14_903_352, 14_303_585, // + 15_314_304, 15_665_330, 15_130_480, // + 8_010_334, 8_393_551, 8_339_378, // + 11_633_793, 12_239_693, 12_206_053, // + 13_491_008, 14_274_597, 14_340_512, + ] + ); + } + + /// Private-tables vs key-sharded reference counting on chr20: identical + /// results, with wall times for the strategy comparison (run with + /// `-- --nocapture`). + #[test] + #[ignore] + fn chr20_private_vs_sharded_strategy() { + let shape = Shape::parse("12of19").unwrap(); + let (_, records, _) = sequence::read_records(chr20_args().reference.as_path()).unwrap(); + let (buf, _, block_len) = + sequence::pack(records.iter().map(|(n, s)| (n.as_str(), s.as_slice())), ""); + let threads = crate::run::resolve_threads(0); + let t = Instant::now(); + let a = seed::count_ref_block_parallel(&buf[..block_len], &shape, 1, threads).unwrap(); + let ms_a = t.elapsed().as_secs_f64() * 1000.0; + let t = Instant::now(); + let b = seed::count_ref_block_sharded(&buf[..block_len], &shape, 1, threads).unwrap(); + let ms_b = t.elapsed().as_secs_f64() * 1000.0; + assert_eq!(a, b); + eprintln!( + "strategy: private-tables {ms_a:.0} ms vs key-sharded {ms_b:.0} ms ({threads} threads)" + ); + } +} diff --git a/src/main.rs b/src/main.rs index 6129a7b..2342614 100644 --- a/src/main.rs +++ b/src/main.rs @@ -74,6 +74,7 @@ mod census; mod cli; mod compare; mod gpu; +mod hits_estimate; mod hsp; mod partition; mod plan; @@ -99,6 +100,7 @@ fn main() { Command::Run(args) => run::run(&args, pre_main_ms, started).map(|_| ()), Command::Benchmark(args) => benchmark::benchmark(&args, pre_main_ms), Command::Compare(args) => compare::compare(&args), + Command::HitsEstimate(args) => hits_estimate::run(&args).map(|_| ()), }; if let Err(e) = result { eprintln!("error: {e}"); diff --git a/src/plan.rs b/src/plan.rs index d391aba..a190d82 100644 --- a/src/plan.rs +++ b/src/plan.rs @@ -1067,7 +1067,6 @@ pub fn host_estimate( threads: usize, max_seeds: usize, ) -> HostEstimate { - let index = (1u64 << (2 * kmer_size)) * 4; // index_table, 4 B/entry let largest_ref = plan .reference_bins .iter() @@ -1080,6 +1079,34 @@ pub fn host_estimate( .map(|b| b.total_bp) .max() .unwrap_or(0); + host_estimate_sizes( + largest_ref, + largest_qry, + ref_bp_total, + qry_bp_total, + kmer_size, + step, + threads, + max_seeds, + ) +} + +/// As [`host_estimate`], with the largest reference and query bins supplied by +/// the caller instead of read from one plan. Batch mode holds every job in RAM +/// and picks both maxima over all jobs, so it calls this directly; the formula +/// itself exists once. +#[allow(clippy::too_many_arguments)] +pub(crate) fn host_estimate_sizes( + largest_ref: u64, + largest_qry: u64, + ref_bp_total: u64, + qry_bp_total: u64, + kmer_size: usize, + step: u32, + threads: usize, + max_seeds: usize, +) -> HostEstimate { + let index = (1u64 << (2 * kmer_size)) * 4; // index_table, 4 B/entry let pos = largest_ref / step.max(1) as u64 * 4; // pos_table, 4 B/indexed bp let counts = threads as u64 * index; // build_parallel pass-1 transient let packed_ref = largest_ref * 2; // buf + enc (no rc for a reference bin) diff --git a/src/run.rs b/src/run.rs index 678f969..4a6edd6 100644 --- a/src/run.rs +++ b/src/run.rs @@ -81,6 +81,9 @@ pub(crate) fn prepare(args: &RunArgs, phases: &mut Phases) -> Fallible .into(), ); } + if args.query_list.is_some() { + return Err("--query-list is only supported by `run`, not benchmark or --cpu-only".into()); + } let plus = args.strand == "plus" || args.strand == "both"; let minus = args.strand == "minus" || args.strand == "both"; @@ -107,7 +110,11 @@ pub(crate) fn prepare(args: &RunArgs, phases: &mut Phases) -> Fallible } else { args.seq_block_size }; - let query = Genome::load(&args.query, &args.query_prefix, seq_block_size)?; + let query_path = args + .query + .as_ref() + .ok_or("--query is required without --query-list")?; + let query = Genome::load(query_path, &args.query_prefix, seq_block_size)?; phases.add("input.query", t.elapsed()); let t = Instant::now(); let reference = Genome::load(&args.reference, &args.target_prefix, seq_block_size)?; @@ -636,11 +643,28 @@ pub(crate) struct Emitter { impl Emitter { pub(crate) fn new(args: &RunArgs) -> Fallible { - let sink: Box = match tarball_path(args) { - Some(path) => Box::new(TarGzSink::new(&path)?), - None => Box::new(DirectorySink::new(&args.output)?), + Self::new_at( + &args.output, + tarball_path(args).as_deref(), + args.diagonal_partition, + crate::census::SurvivorAudit::dump_path().as_deref(), + ) + } + + /// One emitter rooted at `output`: a directory sink, or a tarball sink at + /// `tarball` when set. `audit_path` is the `HSPZ_ANCHOR_CENSUS` dump file; + /// batch jobs pass a per-job path so emitters never share one file. + pub(crate) fn new_at( + output: &std::path::Path, + tarball: Option<&std::path::Path>, + diagonal: bool, + audit_path: Option<&std::path::Path>, + ) -> Fallible { + let sink: Box = match tarball { + Some(path) => Box::new(TarGzSink::new(path)?), + None => Box::new(DirectorySink::new(output)?), }; - let mut audit = crate::census::SurvivorAudit::dump_path() + let mut audit = audit_path .map(|path| std::fs::File::create(path).map(BufWriter::new)) .transpose()?; if let Some(out) = audit.as_mut() { @@ -652,7 +676,7 @@ impl Emitter { Ok(Emitter { sink, part: Partitioner::default(), - diagonal: args.diagonal_partition, + diagonal, partition_ms: 0.0, format_ms: 0.0, archive_ms: 0.0, @@ -783,9 +807,196 @@ impl Emitter { } } +/// Writes a `--dump-plan` membership table: one `sidebinrecordbp` +/// line per record, reference side first. Shared by the single-query path and +/// the batch executor, which calls it once per job at a per-job sibling path. +fn write_plan_dump( + path: &std::path::Path, + plan: &plan::Plan, + ref_records: &[(String, Vec)], + qry_records: &[(String, Vec)], +) -> Fallible<()> { + let mut f = std::io::BufWriter::new(std::fs::File::create(path)?); + for (side, bins, recs) in [ + ("reference", &plan.reference_bins, ref_records), + ("query", &plan.query_bins, qry_records), + ] { + for b in bins { + for &id in &b.record_ids { + let (name, seq) = &recs[id as usize]; + writeln!(f, "{side}\t{}\t{name}\t{}", b.id, seq.len())?; + } + } + } + f.flush()?; + Ok(()) +} + +/// Builds the manifest `--dump-manifest` writes for a fresh (non-replay) run. +/// `executable_hash` is passed in because the batch computes it once for every +/// job, while the single-query path calls [`plan::executable_hash`] only here +/// so an ordinary run pays no hashing overhead. +#[allow(clippy::too_many_arguments)] +fn plan_manifest( + args: &RunArgs, + contract: &crate::gpu::ExecutionContract, + sub_mat: &[i32], + res_seq: u64, + res_qry: u64, + ref_records: &[(String, Vec)], + qry_records: &[(String, Vec)], + plan: &plan::Plan, + executable_hash: u64, +) -> plan::PlanManifest { + plan::PlanManifest { + version: plan::PlanManifest::FORMAT, + hspz_version: env!("CARGO_PKG_VERSION").into(), + features: plan::compiled_features(), + max_hits: contract.max_hits, + hsp_blocks: contract.hsp_blocks, + seed: args.seed.clone(), + step: args.step, + transitions: !args.notransition, + xdrop: args.xdrop, + hspthresh: args.hspthresh, + noentropy: args.noentropy, + wga_chunk_size: args.wga_chunk_size, + lastz_interval_size: args.lastz_interval_size, + kegalign_bins: args.kegalign_bins, + seq_block_size: res_seq, + query_block_size: res_qry, + ref_hash: plan::records_hash(ref_records), + qry_hash: plan::records_hash(qry_records), + executable_hash, + sub_mat: sub_mat.to_vec(), + strand: args.strand.clone(), + target_prefix: args.target_prefix.clone(), + query_prefix: args.query_prefix.clone(), + plan: plan.clone(), + } +} + +/// Writes one `--dump-manifest` file and flushes it: `BufWriter::drop` +/// discards a failed final flush, turning a truncated manifest into a +/// silently "successful" dump. +fn write_manifest_dump(path: &std::path::Path, m: &plan::PlanManifest) -> Fallible<()> { + let mut f = std::io::BufWriter::new(std::fs::File::create(path)?); + m.write(&mut f)?; + f.flush()?; + Ok(()) +} + +/// Everything the `--time` footer reports after the executor returns. One +/// printer for the single-query and batch paths, so the two can never drift; +/// the callers pass the totals that apply to them. +struct TimeFooter<'a> { + phases: &'a Phases, + wall_ms: f64, + launches: u64, + stage_syncs: u64, + pipeline_syncs: u64, + contract: &'a crate::gpu::ExecutionContract, + ref_bins: usize, + units: usize, + lifecycle: &'a Lifecycle, + uploads: u64, + copy_stalls: u64, + files: usize, + bytes_in: u64, + bytes_out: u64, + diagonal: bool, + host_peak_est: u64, + est_shared: u64, + workers: usize, + est_prefetch: u64, + est_no_prefetch: u64, + host_budget: Option, + host_status: &'a str, + prefetch_requested: bool, + prefetch: bool, +} + +impl TimeFooter<'_> { + fn print(&self) { + eprintln!( + "\nWALL-TIME ACCOUNTING\n{}", + self.phases.report(self.wall_ms) + ); + eprintln!(" kernel launches: {}", self.launches); + // Phase 1 §12: the mechanism gate. `stage` waits are the ones stream + // ordering makes unnecessary and are 0 with --async-stages. + eprintln!( + " host syncs: {} stage, {} pipeline ({} per launch)", + self.stage_syncs, + self.pipeline_syncs, + if self.launches > 0 { + format!( + "{:.3}", + (self.stage_syncs + self.pipeline_syncs) as f64 / self.launches as f64 + ) + } else { + "-".into() + } + ); + eprintln!( + " max_hits: {} (target; resolved once; pinned unless --max-hits 0)\n hit_capacity: {} (physical; success/failure only, never output bytes)\n lifecycle: {} ref bins, \ + {} work units, {} builds, {} engines, {} ref uploads, {} query swaps", + self.contract.max_hits, + self.contract.hit_capacity, + self.ref_bins, + self.units, + self.lifecycle.seed_table_builds, + self.lifecycle.engine_creations, + self.lifecycle.reference_uploads, + self.lifecycle.query_swaps, + ); + // Phase 3 mechanism: a stalled upload is one that had not finished when + // its compute needed it, i.e. overlap that did not happen. + eprintln!( + " seed uploads: {} ({} stalled{})", + self.uploads, + self.copy_stalls, + if self.uploads > 0 { + format!( + ", {:.2}%", + self.copy_stalls as f64 / self.uploads as f64 * 100.0 + ) + } else { + String::new() + } + ); + eprintln!( + " output: {} files, {} bytes formatted, {} bytes written{}", + self.files, + self.bytes_in, + self.bytes_out, + if self.diagonal { " (-D)" } else { "" } + ); + eprintln!(" peak RSS: {:>10} KiB", timing::peak_rss_kib()); + // Phase 1 §9: the host-budget decision, in the same units as the line + // above so §11's validation is a subtraction. + let mib = |b: u64| b as f64 / 1048576.0; + eprintln!( + " host budget: estimated peak {:.0} MiB (shared {:.0} + {} worker(s), \ + {:.0} prefetching / {:.0} not), budget {}, status {}, \ + prefetch requested {} effective {}", + mib(self.host_peak_est), + mib(self.est_shared), + self.workers, + mib(self.est_prefetch), + mib(self.est_no_prefetch), + self.host_budget + .map_or("unknown".to_string(), |b| format!("{:.0} MiB", mib(b))), + self.host_status, + self.prefetch_requested, + self.prefetch, + ); + } +} + /// Builds the planner's metadata from raw records: `id` and `ordinal` are the /// input index, so `bin.record_ids` indexes straight back into `records`. -fn record_meta(records: &[(String, Vec)]) -> Vec { +pub(crate) fn record_meta(records: &[(String, Vec)]) -> Vec { records .iter() .enumerate() @@ -1094,6 +1305,30 @@ pub fn device_seeds_for(workers: usize) -> bool { } } +/// Packs one reference bin's records and builds its seed table — the host half +/// of a reference build. The single-query worker and the batch executor both +/// call it inline for the first bin and from the prefetch thread for the next, +/// so the two paths cannot drift. +fn build_ref_bin( + rbin: &plan::Bin, + ref_records: &[(String, Vec)], + target_prefix: &str, + shape: &Shape, + step: u32, + threads: usize, +) -> (PackedBin, SeedTable) { + let packed = PackedBin::build( + rbin.record_ids.iter().map(|&id| { + let (n, s) = &ref_records[id as usize]; + (n.as_str(), s.as_slice()) + }), + target_prefix, + false, + ); + let table = SeedTable::build_parallel(&packed.buf[..packed.block_len], shape, step, threads); + (packed, table) +} + /// Runs one worker's visits on `device`, streaming finished units to the /// emitter (§Phase 5: build/upload each visited reference once, reuse it /// across its slice; no GPU is shared for performance). @@ -1137,19 +1372,6 @@ fn run_bins( // while bin k's work units run, so only the first build stays exposed. // Execution order, bin identity and the lifecycle counts are untouched: this // moves *when* the data is built, not what runs or in which order. - let build_ref_bin = |rbin: &plan::Bin| -> (PackedBin, SeedTable) { - let packed = PackedBin::build( - rbin.record_ids.iter().map(|&id| { - let (n, s) = &ref_records[id as usize]; - (n.as_str(), s.as_slice()) - }), - &args.target_prefix, - false, - ); - let table = - SeedTable::build_parallel(&packed.buf[..packed.block_len], shape, args.step, threads); - (packed, table) - }; let mut pending: Option<(PackedBin, SeedTable)> = None; for (visit_index, visit) in visits.iter().enumerate() { @@ -1157,7 +1379,14 @@ fn run_bins( let t = Instant::now(); let (mut packed_ref, table) = match pending.take() { Some(built) => built, - None => build_ref_bin(rbin), + None => build_ref_bin( + rbin, + ref_records, + &args.target_prefix, + shape, + args.step, + threads, + ), }; rep.seed_table_ms += t.elapsed(); rep.lifecycle.seed_table_builds += 1; @@ -1204,7 +1433,15 @@ fn run_bins( let prefetch = next_bin.map(|nb| { scope.spawn(|| { let t = Instant::now(); - (build_ref_bin(nb), t.elapsed()) + let built = build_ref_bin( + nb, + ref_records, + &args.target_prefix, + shape, + args.step, + threads, + ); + (built, t.elapsed()) }) }); @@ -1386,6 +1623,9 @@ pub(crate) fn run(args: &RunArgs, pre_main_ms: f64, started: Instant) -> Fallibl if args.cpu_only { return run_cpu_only(args, &mut phases, pre_main_ms, started); } + if args.query_list.is_some() { + return run_batch(args, &mut phases, pre_main_ms, started); + } // Shared config, parsed once. let shape = Shape::parse(&args.seed)?; @@ -1393,8 +1633,12 @@ pub(crate) fn run(args: &RunArgs, pre_main_ms: f64, started: Instant) -> Fallibl let transitions = !args.notransition; // Load both sides as raw records — no block-size guard (§10). + let query_path = args + .query + .as_ref() + .ok_or("--query is required without --query-list")?; let t = Instant::now(); - let (_, qry_records, _) = sequence::read_records(&args.query)?; + let (_, qry_records, _) = sequence::read_records(query_path)?; phases.add("input.query", t.elapsed()); let t = Instant::now(); let (_, ref_records, _) = sequence::read_records(&args.reference)?; @@ -1551,61 +1795,30 @@ pub(crate) fn run(args: &RunArgs, pre_main_ms: f64, started: Instant) -> Fallibl phases.add("plan", t.elapsed()); if let Some(path) = &args.dump_plan { - let mut f = std::io::BufWriter::new(std::fs::File::create(path)?); - for (side, bins, recs) in [ - ("reference", &plan.reference_bins, &ref_records), - ("query", &plan.query_bins, &qry_records), - ] { - for b in bins { - for &id in &b.record_ids { - let (name, seq) = &recs[id as usize]; - writeln!(f, "{side}\t{}\t{name}\t{}", b.id, seq.len())?; - } - } - } - f.flush()?; + write_plan_dump(path, &plan, &ref_records, &qry_records)?; } if let Some(path) = &args.dump_manifest { - let mut f = std::io::BufWriter::new(std::fs::File::create(path)?); match &loaded_manifest { // Replay + dump: re-emit exactly what was loaded and validated, // never a reconstruction from this run's B/Q CLI defaults. - Some(m) => m.write(&mut f)?, + Some(m) => write_manifest_dump(path, m)?, None => { // Fresh dump only: this is the one place `plan::executable_hash` // runs, so an ordinary run pays no hashing overhead. - let m = plan::PlanManifest { - version: plan::PlanManifest::FORMAT, - hspz_version: env!("CARGO_PKG_VERSION").into(), - features: plan::compiled_features(), - max_hits: contract.max_hits, - hsp_blocks: contract.hsp_blocks, - seed: args.seed.clone(), - step: args.step, - transitions: !args.notransition, - xdrop: args.xdrop, - hspthresh: args.hspthresh, - noentropy: args.noentropy, - wga_chunk_size: args.wga_chunk_size, - lastz_interval_size: args.lastz_interval_size, - kegalign_bins: args.kegalign_bins, - seq_block_size: res_seq, - query_block_size: res_qry, - ref_hash: plan::records_hash(&ref_records), - qry_hash: plan::records_hash(&qry_records), - executable_hash: plan::executable_hash()?, - sub_mat: sub_mat.clone(), - strand: args.strand.clone(), - target_prefix: args.target_prefix.clone(), - query_prefix: args.query_prefix.clone(), - plan: plan.clone(), - }; - m.write(&mut f)?; + let m = plan_manifest( + args, + &contract, + &sub_mat, + res_seq, + res_qry, + &ref_records, + &qry_records, + &plan, + plan::executable_hash()?, + ); + write_manifest_dump(path, &m)?; } } - // `BufWriter::drop` discards a failed final flush, which would turn a - // truncated manifest into a silently "successful" dump. - f.flush()?; } // §18: reference bins to workers, deterministic LPT (whole-bin), or the @@ -2050,76 +2263,37 @@ pub(crate) fn run(args: &RunArgs, pre_main_ms: f64, started: Instant) -> Fallibl report_counts(&stats); if args.time { let wall = pre_main_ms + started.elapsed().as_secs_f64() * 1000.0; - eprintln!("\nWALL-TIME ACCOUNTING\n{}", phases.report(wall)); - eprintln!(" kernel launches: {}", launches); - // Phase 1 §12: the mechanism gate. `stage` waits are the ones stream - // ordering makes unnecessary and are 0 with --async-stages. - eprintln!( - " host syncs: {} stage, {} pipeline ({} per launch)", + TimeFooter { + phases: &phases, + wall_ms: wall, + launches, stage_syncs, pipeline_syncs, - if launches > 0 { - format!( - "{:.3}", - (stage_syncs + pipeline_syncs) as f64 / launches as f64 - ) - } else { - "-".into() - } - ); - eprintln!( - " max_hits: {} (target; resolved once; pinned unless --max-hits 0)\n hit_capacity: {} (physical; success/failure only, never output bytes)\n lifecycle: {} ref bins, \ -{} work units, {} builds, {} engines, {} ref uploads, {} query swaps", - contract.max_hits, - contract.hit_capacity, - plan.reference_bins.len(), - plan.units.len(), - lifecycle.seed_table_builds, - lifecycle.engine_creations, - lifecycle.reference_uploads, - lifecycle.query_swaps, - ); - // Phase 3 mechanism: a stalled upload is one that had not finished when - // its compute needed it, i.e. overlap that did not happen. - eprintln!( - " seed uploads: {} ({} stalled{})", + contract: &contract, + ref_bins: plan.reference_bins.len(), + units: plan.units.len(), + lifecycle: &lifecycle, uploads, copy_stalls, - if uploads > 0 { - format!(", {:.2}%", copy_stalls as f64 / uploads as f64 * 100.0) - } else { - String::new() - } - ); - eprintln!( - " output: {} files, {} bytes formatted, {} bytes written{}", - out.files, - out.bytes_in, - if out.bytes_out > 0 { + files: out.files, + bytes_in: out.bytes_in, + bytes_out: if out.bytes_out > 0 { out.bytes_out } else { out.bytes_in }, - if args.diagonal_partition { " (-D)" } else { "" } - ); - eprintln!(" peak RSS: {:>10} KiB", timing::peak_rss_kib()); - // Phase 1 §9: the host-budget decision, in the same units as the line - // above so §11's validation is a subtraction. - let mib = |b: u64| b as f64 / 1048576.0; - eprintln!( - " host budget: estimated peak {:.0} MiB (shared {:.0} + {} worker(s), \ - {:.0} prefetching / {:.0} not), budget {}, status {}, \ - prefetch requested {} effective {}", - mib(host_peak_est), - mib(est.shared), + diagonal: args.diagonal_partition, + host_peak_est, + est_shared: est.shared, workers, - mib(est.per_worker_prefetch), - mib(est.per_worker_no_prefetch), - host_budget.map_or("unknown".to_string(), |b| format!("{:.0} MiB", mib(b))), + est_prefetch: est.per_worker_prefetch, + est_no_prefetch: est.per_worker_no_prefetch, + host_budget, host_status, prefetch_requested, prefetch, - ); + } + .print(); } if let Some(a) = audit.as_ref() { eprintln!("\n(ALL REFERENCE BINS) {}", a.report()); @@ -2130,6 +2304,777 @@ pub(crate) fn run(args: &RunArgs, pre_main_ms: f64, started: Instant) -> Fallibl Ok(stats) } +// --------------------------------------------------------------------------- +// Bin-major batch: one reference × many queries (`--query-list`), W=1 only. +// +// Every job is planned independently exactly as `run -q` would (same +// `plan_within_budget` call, same once-resolved `--max-hits`), and the batch +// aborts before any GPU work unless all jobs share identical reference bins. +// Execution is bin-major: each reference bin is built and uploaded once (with +// the existing one-bin-ahead prefetch), then every job's query blocks for that +// bin run through the unchanged `seed_and_filter_all` on their original units. +// Each job owns one `Emitter` (+ `-D` history) fed in its standalone ordinal +// order via a per-job cursor in `JobRouter`. +// +// A batch is all-or-nothing: any job's failure (`?`, OOM, CUDA, router, pack) +// unwinds `run_batch`. Per-job directories are left with whatever files were +// already written, a `-Z` archive is left truncated (`TarGzSink::finish` never +// runs), later jobs on later bins never run, and `OUT/queries.tsv` is only +// written after every emitter finishes. This is documented in the +// `--query-list` CLI help as well. + +/// Batch-mode flag validation, pure so it is unit-testable without a GPU. +pub(crate) fn validate_batch_args(args: &RunArgs) -> Result<(), String> { + if args.gpus > 1 { + return Err("batch mode is W=1 in this release".into()); + } + if args.seq_block_size == 0 { + return Err("batch mode rejects -B 0 (automatic layout): give -B explicitly".into()); + } + if args.kegalign_bins { + return Err("batch mode rejects --kegalign-bins".into()); + } + if args.from_manifest.is_some() { + return Err("batch mode rejects --from-manifest".into()); + } + Ok(()) +} + +/// Reads a `--query-list` file: one query FASTA path per line, blank lines and +/// `#` comments ignored. Relative paths resolve against the list file's dir. +pub(crate) fn read_query_list(path: &std::path::Path) -> Fallible> { + let text = std::fs::read_to_string(path)?; + let base = path.parent().filter(|p| !p.as_os_str().is_empty()); + let mut out = Vec::new(); + for line in text.lines() { + let line = line.trim(); + if line.is_empty() || line.starts_with('#') { + continue; + } + let p = PathBuf::from(line); + out.push(match base { + Some(dir) if p.is_relative() => dir.join(p), + _ => p, + }); + } + if out.is_empty() { + return Err(format!("query list {} has no queries", path.display()).into()); + } + Ok(out) +} + +/// Reference-bin compatibility across batch jobs: identical count, ids, +/// membership and sizes. Anything else would silently replan one job into a +/// different dedup scope, so the batch aborts naming the job instead. +pub(crate) fn check_reference_bins_compatible( + first: &[plan::Bin], + job: &[plan::Bin], + job_index: usize, + job_path: &std::path::Path, +) -> Result<(), String> { + if first.len() != job.len() + || first + .iter() + .zip(job.iter()) + .any(|(a, b)| a.id != b.id || a.total_bp != b.total_bp || a.record_ids != b.record_ids) + { + return Err(format!( + "batch: job {:06} ({}): reference bins differ from job 000001; \ + batch only admits jobs with identical reference bins", + job_index + 1, + job_path.display() + )); + } + Ok(()) +} + +/// Output router: one next-ordinal cursor per job. A single global cursor +/// would buffer every later job while bin-major execution waits on the first, +/// so each job advances in its own standalone ordinal order. +pub(crate) struct JobRouter { + next: Vec, + totals: Vec, +} + +impl JobRouter { + pub(crate) fn new(totals: Vec) -> Self { + Self { + next: vec![0; totals.len()], + totals, + } + } + + /// Admits `(job, ordinal)` iff it is that job's next expected unit *and* + /// the bin tuple matches the job plan's unit at that ordinal. + pub(crate) fn accept( + &mut self, + job: usize, + ordinal: u32, + reference_bin: u32, + query_bin: u32, + plan: &plan::Plan, + ) -> Result<(), String> { + if job >= self.next.len() { + return Err(format!("batch router: unknown job {job}")); + } + let want = self.next[job]; + if ordinal < want { + return Err(format!( + "batch router: duplicate work unit ordinal {ordinal} for job {:06}", + job + 1 + )); + } + if ordinal > want { + return Err(format!( + "batch router: missing work unit for job {:06}: expected ordinal {want}, got {ordinal}", + job + 1 + )); + } + if ordinal >= self.totals[job] { + return Err(format!( + "batch router: work unit ordinal {ordinal} is outside job {:06}'s plan ({} units)", + job + 1, + self.totals[job] + )); + } + let pu = &plan.units[ordinal as usize]; + if (pu.reference_bin, pu.query_bin) != (reference_bin, query_bin) { + return Err(format!( + "batch router: job {:06} ordinal {ordinal} arrived as R{reference_bin} Q{query_bin} \ + but the plan has R{} Q{}", + job + 1, + pu.reference_bin, + pu.query_bin + )); + } + self.next[job] += 1; + Ok(()) + } + + /// True once every job's cursor reached its plan length. + pub(crate) fn complete(&self) -> bool { + self.next == self.totals + } +} + +/// Per-job sibling of a batch dump path: append `.{job:06}` to the full file +/// name, so `--dump-manifest out.manifest` writes `out.manifest.000001` and +/// `--dump-plan`/`--dump-raw`/the `HSPZ_ANCHOR_CENSUS` dump all agree. Never +/// `Path::with_extension`, which would replace `.manifest` with `.000001`. +/// `job` is the zero-based loop index. +pub(crate) fn job_sibling(path: &std::path::Path, job: usize) -> PathBuf { + let mut name = path.as_os_str().to_os_string(); + name.push(format!(".{:06}", job + 1)); + PathBuf::from(name) +} + +/// SHA-256 of a file via the existing `sha256sum` binary (no new crates). +fn sha256_file(path: &std::path::Path) -> Fallible { + let out = std::process::Command::new("sha256sum").arg(path).output()?; + if !out.status.success() { + return Err(format!("sha256sum failed for {}", path.display()).into()); + } + let text = String::from_utf8_lossy(&out.stdout); + text.split_whitespace() + .next() + .map(str::to_string) + .ok_or_else(|| format!("sha256sum gave no digest for {}", path.display()).into()) +} + +#[allow(clippy::too_many_lines)] +fn run_batch( + args: &RunArgs, + phases: &mut Phases, + pre_main_ms: f64, + started: Instant, +) -> Fallible { + validate_batch_args(args).map_err(|e| -> Box { e.into() })?; + let list_path = args + .query_list + .as_ref() + .ok_or("--query-list is required for batch mode")?; + let job_paths = read_query_list(list_path)?; + + let shape = Shape::parse(&args.seed)?; + let sub_mat = scoring::build_sub_mat(&args.ambiguous, args.xdrop, args.scoring.as_deref())?; + let transitions = !args.notransition; + + let t = Instant::now(); + let (_, ref_records, ref_bytes) = sequence::read_records(&args.reference)?; + phases.add("input.reference", t.elapsed()); + let ref_meta = record_meta(&ref_records); + let ref_bp_total: u64 = ref_meta.iter().map(|r| r.len).sum(); + + // Resolve `--max-hits 0` ONCE for the whole batch (design: never from + // whichever worker happens to run a unit — here there is one device). + let t = Instant::now(); + let ctx = CudaContext::new(0)?; + phases.add("CUDA context init", t.elapsed()); + // `plan` mirrors the single-query path exactly: device probe, worker budget + // and contract resolution, then each job's own `plan_within_budget`, then + // the shared hit capacity. Never the CUDA context or the query loads — the + // query reads are charged to `input.query` (once per job), as standalone. + let t = Instant::now(); + let devices = crate::gpu::device_count().max(1); + let free = crate::gpu::min_free_bytes(1)?; + let workers_upper = 1usize.min(ref_meta.len().max(1)); + let budget = plan::worker_device_budget(free, workers_upper, devices); + let mut contract = crate::gpu::ExecutionContract::resolve(&ctx, args.max_hits, args.hsp_blocks); + let resolved_max_hits = contract.max_hits; + let q_target = args + .query_block_size + .map(u64::from) + .unwrap_or(args.seq_block_size as u64); + let res_seq = args.seq_block_size as u64; + let mut plan_ms = t.elapsed(); + + // Plan EVERY job independently with the same call the single path uses. + let mut qry_records_list: Vec)>> = Vec::with_capacity(job_paths.len()); + let mut plans: Vec = Vec::with_capacity(job_paths.len()); + let mut worsts: Vec = Vec::with_capacity(job_paths.len()); + for path in &job_paths { + let tq = Instant::now(); + let (_, recs, _) = sequence::read_records(path)?; + phases.add("input.query", tq.elapsed()); + let tp = Instant::now(); + let meta = record_meta(&recs); + let (p, worst) = plan::plan_within_budget( + &ref_meta, + &meta, + res_seq, + q_target, + budget, + shape.kmer_size, + args.step, + resolved_max_hits, + false, + args.wga_chunk_size, + transitions, + ) + .map_err(|e| { + format!( + "batch: job {} ({}): {e}", + qry_records_list.len() + 1, + path.display() + ) + })?; + qry_records_list.push(recs); + plans.push(p); + worsts.push(worst); + plan_ms += tp.elapsed(); + } + let tp = Instant::now(); + // The reference bins must be identical across all jobs — abort before any + // GPU work, naming the job. + for (j, (p, path)) in plans.iter().zip(&job_paths).enumerate().skip(1) { + check_reference_bins_compatible(&plans[0].reference_bins, &p.reference_bins, j, path)?; + } + + // Physical capacity shared by the batch's engines: per-plan candidates + // never lower the semantic cap, so take the minimum that still fits every + // job (success/failure only, never output bytes). + { + let mut shared: Option = None; + for p in &plans { + let candidate = plan::max_hit_capacity( + p, + budget, + shape.kmer_size, + args.step, + resolved_max_hits, + args.wga_chunk_size, + transitions, + )?; + shared = Some(shared.map_or(candidate, |s: u32| s.min(candidate))); + } + contract.hit_capacity = crate::gpu::clamp_hit_capacity( + contract.max_hits, + shared.unwrap_or(contract.max_hits), + contract.hsp_blocks, + )?; + } + plan_ms += tp.elapsed(); + phases.add("plan", plan_ms); + + if let Some(path) = &args.dump_plan { + for (j, p) in plans.iter().enumerate() { + write_plan_dump(&job_sibling(path, j), p, &ref_records, &qry_records_list[j])?; + } + } + if let Some(path) = &args.dump_manifest { + // One executable hash for every job's manifest. Batch rejects + // `--from-manifest`, so every manifest here is a fresh dump. + let executable_hash = plan::executable_hash()?; + for (j, p) in plans.iter().enumerate() { + let m = plan_manifest( + args, + &contract, + &sub_mat, + res_seq, + q_target, + &ref_records, + &qry_records_list[j], + p, + executable_hash, + ); + write_manifest_dump(&job_sibling(path, j), &m)?; + } + } + + // Host preflight over the WHOLE batch: all queries' records stay in RAM, + // packed per unit from the held records (never re-read per bin). + let qry_bp_sum: u64 = qry_records_list + .iter() + .map(|recs| recs.iter().map(|(_, s)| s.len() as u64).sum::()) + .sum(); + let largest_ref = plans[0] + .reference_bins + .iter() + .map(|b| b.total_bp) + .max() + .unwrap_or(0); + let largest_qry = plans + .iter() + .flat_map(|p| p.query_bins.iter().map(|b| b.total_bp)) + .max() + .unwrap_or(0); + let threads = resolve_threads(args.threads); + let max_seeds = seed::max_seeds(args.wga_chunk_size, &shape, transitions); + // Shared formula with the single-query path; the two maxima span every + // job here, because the batch holds all of them in RAM at once. + let est = plan::host_estimate_sizes( + largest_ref, + largest_qry, + ref_bp_total, + qry_bp_sum, + shape.kmer_size, + args.step, + threads, + max_seeds, + ); + let ref_bins_all: Vec = plans[0].reference_bins.iter().map(|b| b.id).collect(); + let assignment: Vec> = vec![ref_bins_all]; + let prefetch_requested = !args.no_ref_prefetch; + let mut prefetch = prefetch_requested; + let mut host_budget = None; + let mut host_status = "unknown (no cgroup/meminfo reading)"; + if let Some(available) = timing::available_host_bytes().map(|b| b * 9 / 10) { + host_budget = Some(available); + let fits = plan::host_preflight(&est, &assignment, available)?; + host_status = if fits { + "fits with prefetch" + } else { + "fits without prefetch" + }; + if prefetch && !fits { + eprintln!("note: host preflight disabled reference prefetch (batch W=1)"); + } + prefetch &= fits; + } + let host_peak_est = plan::host_peak(&est, &assignment, prefetch); + + // One emitter (+ `-D` history) per job, each starting empty. + std::fs::create_dir_all(&args.output)?; + let batch_tarball = args.tarball.is_some(); + if let Some(p) = &args.tarball { + // Bare `-Z` (clap's `-` sentinel) and an empty value mean the default + // per-job archive layout; any other path cannot be honoured because a + // batch writes one archive per job. + if !p.as_os_str().is_empty() && p.as_os_str() != "-" { + eprintln!( + "note: batch mode writes one archive per job (OUT/000001.tar.gz …); \ + the -Z path value is ignored" + ); + } + } + let audit_base = crate::census::SurvivorAudit::dump_path(); + let mut emitters: Vec = Vec::with_capacity(job_paths.len()); + for (j, _) in job_paths.iter().enumerate() { + let tag = format!("{:06}", j + 1); + let (dir, tar) = if batch_tarball { + ( + args.output.join(&tag), + Some(args.output.join(format!("{tag}.tar.gz"))), + ) + } else { + (args.output.join(&tag), None) + }; + let audit_path = audit_base.as_ref().map(|p| job_sibling(p, j)); + emitters.push(Emitter::new_at( + &dir, + tar.as_deref(), + args.diagonal_partition, + audit_path.as_deref(), + )?); + } + let mut router = JobRouter::new(plans.iter().map(|p| p.units.len() as u32).collect()); + let mut job_stats: Vec = vec![Stats::default(); job_paths.len()]; + let mut job_raw: Vec)>> = + (0..job_paths.len()).map(|_| Vec::new()).collect(); + let mut job_first_ms: Vec> = vec![None; job_paths.len()]; + let mut job_done_ms: Vec = vec![0.0; job_paths.len()]; + let mut job_hit_stats: Vec = (0..job_paths.len()) + .map(|_| crate::gpu::HitStats::default()) + .collect(); + let mut job_audits: Vec> = + (0..job_paths.len()).map(|_| None).collect(); + + let total_units: usize = plans.iter().map(|p| p.units.len()).sum(); + let r_bins = plans[0].reference_bins.len(); + eprintln!( + "schedule: policy=whole-bin W=1 visits={r_bins} extra_replicas=0 (batch W=1; auto: W=1)" + ); + { + let desc = plans[0] + .reference_bins + .iter() + .map(|b| format!("R{}(Qbatch)", b.id)) + .collect::>() + .join(" "); + eprintln!("schedule: worker 0 device 0: {desc} units={total_units} visits={r_bins}"); + } + eprintln!("schedule: device mapping: worker0->device0"); + + let device_seeds = device_seeds_for(1); + if device_seeds { + eprintln!("note: generating query seeds on the device (1 worker)"); + } + + let mut pending: Option<(plan::PackedBin, SeedTable)> = None; + let mut lifecycle = Lifecycle::default(); + let mut launches = 0u64; + let (mut stage_syncs, mut pipeline_syncs) = (0u64, 0u64); + let (mut uploads, mut copy_stalls) = (0u64, 0u64); + let mut seed_table_ms = Duration::ZERO; + let mut prefetched_ms = Duration::ZERO; + + for (visit_index, rbin) in plans[0].reference_bins.iter().enumerate() { + let t = Instant::now(); + let (mut packed_ref, table) = match pending.take() { + Some(built) => built, + None => build_ref_bin( + rbin, + &ref_records, + &args.target_prefix, + &shape, + args.step, + threads, + ), + }; + seed_table_ms += t.elapsed(); + lifecycle.seed_table_builds += 1; + + let cfg = EngineConfig { + index_table: &table.index_table, + pos_table: &table.pos_table, + ref_seq: &packed_ref.enc, + sub_mat: &sub_mat, + seed_size: shape.size as u32, + xdrop: args.xdrop, + hspthresh: args.hspthresh, + noentropy: args.noentropy, + max_hits: contract.max_hits, + hit_capacity: contract.hit_capacity, + timing: args.time, + hsp_blocks: contract.hsp_blocks, + }; + let mut engine = Engine::new(&ctx, cfg, &mut *phases)?; + lifecycle.engine_creations += 1; + engine.dump_raw = args.dump_raw.is_some(); + engine.device_seeds = device_seeds; + engine.collect_hit_stats = args.hit_stats; + engine.persistent_seed_buffers = !args.no_persistent_seed_buffers; + engine.async_stages = !args.no_async_stages; + + let ref_chrs = std::mem::take(&mut packed_ref.chrs); + drop(packed_ref); + drop(table); + + // The next reference bin builds on a worker thread while this bin's + // units run on the GPU — the existing one-bin-ahead prefetch, kept at + // W=1 (`std::thread::scope` so the closure can borrow the records). + let next_bin = plans[0] + .reference_bins + .get(visit_index + 1) + .filter(|_| prefetch); + std::thread::scope(|scope| -> Fallible<()> { + let prefetch = next_bin.map(|nb| { + scope.spawn(|| { + let t = Instant::now(); + let built = build_ref_bin( + nb, + &ref_records, + &args.target_prefix, + &shape, + args.step, + threads, + ); + (built, t.elapsed()) + }) + }); + + // Bin-major: for this reference bin, every job's units at this bin + // in that job's ordinal order (ref-outer ordinals make this + // per-job ordered, so each job's `-D` history matches standalone). + for (j, plan_j) in plans.iter().enumerate() { + let units: Vec = plan_j + .units + .iter() + .filter(|u| u.reference_bin == rbin.id) + .copied() + .collect(); + for unit in units { + let qbin = &plan_j.query_bins[unit.query_bin as usize]; + let t = Instant::now(); + let packed_q = plan::PackedBin::build( + qbin.record_ids.iter().map(|&id| { + let (n, s) = &qry_records_list[j][id as usize]; + (n.as_str(), s.as_slice()) + }), + &args.query_prefix, + true, + ); + phases.add("query pack", t.elapsed()); + let intervals = sequence::intervals( + packed_q.block_len, + shape.size, + args.lastz_interval_size, + ); + let q_block_len = packed_q.block_len.saturating_sub(shape.size) as u32; + + let t = Instant::now(); + engine.swap_query(&packed_q.enc, &packed_q.enc_rc)?; + phases.add("swap_query", t.elapsed()); + let qpass = QueryPass { + fwd: &packed_q.buf[..packed_q.block_len], + rc: &packed_q.rc, + intervals: &intervals, + q_block_len, + }; + let pass = seed_and_filter_all( + &mut engine, + &qpass, + &shape, + transitions, + args, + threads, + ) + .map_err(|e| { + format!( + "batch: job {:06} unit {} ref_bin {} query_bin {}: {e}", + j + 1, + unit.ordinal, + rbin.id, + qbin.id + ) + })?; + // Per-job diagnostic attribution (byte-neutral): drain what + // this unit added into its job's accumulators. + if args.hit_stats { + job_hit_stats[j].merge(&std::mem::take(&mut engine.hit_stats)); + } + if let Some(c) = engine.census.as_mut() { + job_audits[j] + .get_or_insert_with(crate::census::SurvivorAudit::default) + .merge(&std::mem::take(c)); + } + + router + .accept(j, unit.ordinal, unit.reference_bin, unit.query_bin, plan_j) + .map_err(|e| -> Box { e.into() })?; + emitters[j].emit_unit( + unit.reference_bin, + unit.query_bin, + &ref_chrs, + &packed_q.chrs, + &packed_q.rc_chrs, + &pass, + )?; + if args.dump_raw.is_some() { + job_raw[j].extend_from_slice(&pass.raw); + } + let now_ms = started.elapsed().as_secs_f64() * 1000.0; + if job_first_ms[j].is_none() { + job_first_ms[j] = Some(now_ms); + } + job_done_ms[j] = now_ms; + job_stats[j].seeds += pass.stats.seeds; + job_stats[j].seed_hits += pass.stats.seed_hits; + job_stats[j].raw_hsps += pass.stats.raw_hsps; + job_stats[j].hsps += pass.stats.hsps; + lifecycle.work_units_executed += 1; + lifecycle.query_swaps += 1; + } + } + + if let Some(prefetch) = prefetch { + let t = Instant::now(); + let (built, standalone) = prefetch.join().expect("reference prefetch panicked"); + seed_table_ms += t.elapsed(); + prefetched_ms += standalone; + pending = Some(built); + } + Ok(()) + })?; + + let got_swaps = engine.query_swaps(); + let want_swaps: u32 = plans + .iter() + .map(|p| { + p.units + .iter() + .filter(|u| u.reference_bin == rbin.id) + .count() as u32 + }) + .sum(); + if got_swaps != want_swaps { + return Err(format!( + "batch: reference bin {}: engine query_swaps = {got_swaps}, expected {want_swaps}", + rbin.id + ) + .into()); + } + lifecycle.reference_uploads += engine.reference_uploads(); + engine.finish_bucket_autotune(); + launches += engine.launches; + stage_syncs += engine.stage_syncs(); + pipeline_syncs += engine.pipeline_syncs(); + let (u, st) = engine.seed_copy_stats(); + uploads += u; + copy_stalls += st; + if let Some(c) = engine.census.as_mut() { + eprintln!("\n(reference bin {}) {}", rbin.id, c.report()); + } + phases.merge(&engine.phases); + } + if !router.complete() { + return Err("batch: emitter did not reach every job's final ordinal".into()); + } + lifecycle.check(r_bins as u32, total_units as u32)?; + + phases.add("seed table build", seed_table_ms); + if prefetched_ms > Duration::ZERO { + phases.add_overlapped("reference bin prep (standalone)", prefetched_ms); + } + + if let Some(path) = &args.dump_raw { + for (j, raw) in job_raw.iter().enumerate() { + let mut f = std::io::BufWriter::new(std::fs::File::create(job_sibling(path, j))?); + for (strand, hsps) in raw { + for h in hsps { + writeln!( + f, + "{strand}\t{}\t{}\t{}\t{}", + h.ref_start, h.query_start, h.len, h.score + )?; + } + } + } + } + + let mut total_files = 0usize; + let mut total_bytes_in = 0u64; + let mut total_bytes_out = 0u64; + for emitter in emitters { + let out = emitter.finish(&mut *phases)?; + total_files += out.files; + total_bytes_in += out.bytes_in; + total_bytes_out += out.bytes_out.max(out.bytes_in); + } + + // OUT/queries.tsv: per-job identity and completion. + { + let mut f = + std::io::BufWriter::new(std::fs::File::create(args.output.join("queries.tsv"))?); + writeln!( + f, + "job\tpath\tbytes\tsha256\tblocks\tunits\thsps\tcompletion_ms" + )?; + for (j, path) in job_paths.iter().enumerate() { + let bytes = std::fs::metadata(path).map(|m| m.len()).unwrap_or(0); + let sha = sha256_file(path)?; + writeln!( + f, + "{:06}\t{}\t{bytes}\t{sha}\t{}\t{}\t{}\t{:.2}", + j + 1, + path.display(), + plans[j].query_bins.len(), + plans[j].units.len(), + job_stats[j].hsps, + job_done_ms[j], + )?; + } + f.flush()?; + } + + let mut stats = Stats::default(); + for s in &job_stats { + stats.seeds += s.seeds; + stats.seed_hits += s.seed_hits; + stats.raw_hsps += s.raw_hsps; + stats.hsps += s.hsps; + } + report_counts(&stats); + for (j, s) in job_stats.iter().enumerate() { + eprintln!( + "job {:06}: blocks={} units={} hsps={} first_output_ms={:.2} completed_ms={:.2}", + j + 1, + plans[j].query_bins.len(), + plans[j].units.len(), + s.hsps, + job_first_ms[j].unwrap_or(job_done_ms[j]), + job_done_ms[j], + ); + } + eprintln!( + "batch: jobs={} units={total_units} wall_ms={:.2}", + job_paths.len(), + started.elapsed().as_secs_f64() * 1000.0 + ); + if args.time { + let wall = pre_main_ms + started.elapsed().as_secs_f64() * 1000.0; + TimeFooter { + phases, + wall_ms: wall, + launches, + stage_syncs, + pipeline_syncs, + contract: &contract, + ref_bins: r_bins, + units: total_units, + lifecycle: &lifecycle, + uploads, + copy_stalls, + files: total_files, + bytes_in: total_bytes_in, + bytes_out: total_bytes_out, + diagonal: args.diagonal_partition, + host_peak_est, + est_shared: est.shared, + workers: 1, + est_prefetch: est.per_worker_prefetch, + est_no_prefetch: est.per_worker_no_prefetch, + host_budget, + host_status, + prefetch_requested, + prefetch, + } + .print(); + for (j, a) in job_audits.iter().flatten().enumerate() { + eprintln!("\n(job {:06}, ALL REFERENCE BINS) {}", j + 1, a.report()); + } + if args.hit_stats { + let mut merged = crate::gpu::HitStats::default(); + for h in &job_hit_stats { + merged.merge(h); + } + eprintln!("\nHITS PER SEED\n{}", merged.report()); + } + } + let _ = (ref_bytes, worsts); + Ok(stats) +} + /// The `--cpu-only` path: reader-attributable digests and seed/hit counts, no /// GPU. Kept on `prepare`, so it still rejects oversized single-block input — /// the whole-genome reader check is a separate concern from the executor. @@ -2561,7 +3506,7 @@ mod tests { let mut phases = Phases::new(); let mut auto = base_args(); auto.reference = PathBuf::from("/nonexistent/ref.fa"); - auto.query = PathBuf::from("/nonexistent/qry.fa"); + auto.query = Some(PathBuf::from("/nonexistent/qry.fa")); auto.seq_block_size = 0; auto.kegalign_bins = true; let err = prepare(&auto, &mut phases).err().unwrap(); @@ -2579,6 +3524,259 @@ mod tests { assert_eq!(layout_workers(2, 2, Some(0)), 2, "zero override is ignored"); } + /// Batch validation is pure: W>1, `-B 0`, `--kegalign-bins` and + /// `--from-manifest` are all rejected with clear messages. + #[test] + fn batch_rejects_unsupported_combinations() { + let mut bad_gpus = base_args(); + bad_gpus.query = None; + bad_gpus.query_list = Some(PathBuf::from("q.txt")); + bad_gpus.gpus = 2; + assert_eq!( + validate_batch_args(&bad_gpus).unwrap_err(), + "batch mode is W=1 in this release" + ); + + let mut bad_b = base_args(); + bad_b.query_list = Some(PathBuf::from("q.txt")); + bad_b.seq_block_size = 0; + assert!(validate_batch_args(&bad_b).unwrap_err().contains("-B 0")); + + let mut bad_k = base_args(); + bad_k.query_list = Some(PathBuf::from("q.txt")); + bad_k.kegalign_bins = true; + assert!( + validate_batch_args(&bad_k) + .unwrap_err() + .contains("kegalign") + ); + + let mut bad_m = base_args(); + bad_m.query_list = Some(PathBuf::from("q.txt")); + bad_m.from_manifest = Some(PathBuf::from("m.txt")); + assert!( + validate_batch_args(&bad_m) + .unwrap_err() + .contains("from-manifest") + ); + + let ok = base_args(); + assert!(validate_batch_args(&ok).is_ok()); + } + + /// `--query-list` parsing: blank lines and `#` comments are ignored; + /// relative paths resolve against the list file's directory. + #[test] + fn query_list_ignores_blanks_comments_and_resolves_relative() { + let dir = std::env::temp_dir().join(format!("hspz-qlist-{}", std::process::id())); + std::fs::create_dir_all(&dir).unwrap(); + let list = dir.join("list.txt"); + std::fs::write(&list, "# comment\n\na.fa\n \n# another\nsub/b.fa\n").unwrap(); + let got = read_query_list(&list).unwrap(); + assert_eq!(got, vec![dir.join("a.fa"), dir.join("sub/b.fa")]); + std::fs::remove_dir_all(&dir).unwrap(); + + let dir2 = std::env::temp_dir().join(format!("hspz-qlist-empty-{}", std::process::id())); + std::fs::create_dir_all(&dir2).unwrap(); + let empty = dir2.join("empty.txt"); + std::fs::write(&empty, "# only comments\n\n").unwrap(); + assert!(read_query_list(&empty).is_err()); + std::fs::remove_dir_all(&dir2).unwrap(); + } + + /// Per-job dump siblings append `.{job:06}` to the full file name, so all + /// four batch dump paths agree. The old `Path::with_extension` form dropped + /// `.manifest`/`.raw` and wrote `out.000001` instead of `out.manifest.000001`. + #[test] + fn batch_dump_siblings_append_to_the_full_file_name() { + assert_eq!( + job_sibling(PathBuf::from("out.manifest").as_path(), 0), + PathBuf::from("out.manifest.000001") + ); + assert_eq!( + job_sibling(PathBuf::from("/tmp/out.plan").as_path(), 41), + PathBuf::from("/tmp/out.plan.000042") + ); + assert_eq!( + job_sibling(PathBuf::from("dump.raw").as_path(), 2), + PathBuf::from("dump.raw.000003") + ); + // The HSPZ_ANCHOR_CENSUS dump uses the same helper, not `format!("{}.{}")`. + assert_eq!( + job_sibling(PathBuf::from("out/queries.tsv").as_path(), 9), + PathBuf::from("out/queries.tsv.000010") + ); + } + + /// Reference-bin compatibility: identical bins pass; any difference in + /// count, membership or size aborts naming the job. + #[test] + fn reference_bin_compatibility_names_the_job() { + use crate::plan::Bin; + let bins = vec![ + Bin { + id: 0, + record_ids: vec![0, 1], + total_bp: 100, + }, + Bin { + id: 1, + record_ids: vec![2], + total_bp: 50, + }, + ]; + assert!( + check_reference_bins_compatible( + &bins, + &bins.clone(), + 1, + PathBuf::from("b.fa").as_path() + ) + .is_ok() + ); + let fewer = vec![Bin { + id: 0, + record_ids: vec![0, 1], + total_bp: 100, + }]; + let err = + check_reference_bins_compatible(&bins, &fewer, 1, PathBuf::from("b.fa").as_path()) + .unwrap_err(); + assert!(err.contains("job 000002") && err.contains("b.fa"), "{err}"); + let mut moved = bins.clone(); + moved[0].record_ids = vec![0, 2]; + let err = + check_reference_bins_compatible(&bins, &moved, 2, PathBuf::from("c.fa").as_path()) + .unwrap_err(); + assert!(err.contains("job 000003"), "{err}"); + let mut resized = bins.clone(); + resized[1].total_bp = 51; + assert!( + check_reference_bins_compatible(&bins, &resized, 1, PathBuf::from("b.fa").as_path()) + .is_err() + ); + } + + /// The batch router keeps one next-ordinal cursor per job: a bin-major + /// interleave (R0J0, R0J1, R1J0, R1J1) advances each job in its standalone + /// order, while duplicates, skips, swapped bins and unknown jobs fail. + #[test] + fn job_router_routes_by_job_and_rejects_bad_identity() { + use crate::plan::{Bin, Plan, WorkUnit}; + // Two jobs sharing R=2 bins: job 0 has Q=1 (2 units), job 1 has Q=2 + // (4 units). Ordinals are ref-outer per job. + let bins = || { + vec![ + Bin { + id: 0, + record_ids: vec![0], + total_bp: 10, + }, + Bin { + id: 1, + record_ids: vec![1], + total_bp: 10, + }, + ] + }; + let plan0 = Plan { + reference_bins: bins(), + query_bins: vec![Bin { + id: 0, + record_ids: vec![0], + total_bp: 5, + }], + units: vec![ + WorkUnit { + ordinal: 0, + reference_bin: 0, + query_bin: 0, + }, + WorkUnit { + ordinal: 1, + reference_bin: 1, + query_bin: 0, + }, + ], + }; + let plan1 = Plan { + reference_bins: bins(), + query_bins: vec![ + Bin { + id: 0, + record_ids: vec![0], + total_bp: 5, + }, + Bin { + id: 1, + record_ids: vec![1], + total_bp: 5, + }, + ], + units: vec![ + WorkUnit { + ordinal: 0, + reference_bin: 0, + query_bin: 0, + }, + WorkUnit { + ordinal: 1, + reference_bin: 0, + query_bin: 1, + }, + WorkUnit { + ordinal: 2, + reference_bin: 1, + query_bin: 0, + }, + WorkUnit { + ordinal: 3, + reference_bin: 1, + query_bin: 1, + }, + ], + }; + let mut r = JobRouter::new(vec![2, 4]); + // Bin-major traversal the executor uses. + r.accept(0, 0, 0, 0, &plan0).unwrap(); + r.accept(1, 0, 0, 0, &plan1).unwrap(); + r.accept(1, 1, 0, 1, &plan1).unwrap(); + r.accept(0, 1, 1, 0, &plan0).unwrap(); + r.accept(1, 2, 1, 0, &plan1).unwrap(); + r.accept(1, 3, 1, 1, &plan1).unwrap(); + assert!(r.complete()); + + // Duplicate: job 0 ordinal 1 was already consumed. + let mut r = JobRouter::new(vec![2, 4]); + r.accept(0, 0, 0, 0, &plan0).unwrap(); + r.accept(0, 1, 1, 0, &plan0).unwrap(); + let err = r.accept(0, 1, 1, 0, &plan0).unwrap_err(); + assert!(err.contains("duplicate"), "{err}"); + + // Missing: skipping ordinal 0 for job 1. + let mut r = JobRouter::new(vec![2, 4]); + let err = r.accept(1, 1, 0, 1, &plan1).unwrap_err(); + assert!(err.contains("missing"), "{err}"); + + // Swapped bin tuple at the expected ordinal. + let mut r = JobRouter::new(vec![2, 4]); + let err = r.accept(0, 0, 1, 0, &plan0).unwrap_err(); + assert!(err.contains("R1 Q0") && err.contains("R0 Q0"), "{err}"); + + // Unknown ordinal past the plan end, and unknown job. + let mut r = JobRouter::new(vec![2, 4]); + r.accept(0, 0, 0, 0, &plan0).unwrap(); + r.accept(0, 1, 1, 0, &plan0).unwrap(); + let err = r.accept(0, 2, 0, 0, &plan0).unwrap_err(); + assert!(err.contains("outside"), "{err}"); + assert!( + r.accept(7, 0, 0, 0, &plan0) + .unwrap_err() + .contains("unknown job") + ); + assert!(!r.complete()); + } + /// `prepare` is the GPU-free path `benchmark` and `--cpu-only` share; neither /// implements frozen-plan replay or manifest dumps, so both flags must fail /// fast here rather than being silently dropped. Both paths point at inputs @@ -2589,14 +3787,14 @@ mod tests { let mut from_manifest = base_args(); from_manifest.reference = PathBuf::from("/nonexistent/ref.fa"); - from_manifest.query = PathBuf::from("/nonexistent/qry.fa"); + from_manifest.query = Some(PathBuf::from("/nonexistent/qry.fa")); from_manifest.from_manifest = Some(PathBuf::from("/nonexistent/plan.manifest")); let err = prepare(&from_manifest, &mut phases).err().unwrap(); assert!(err.to_string().contains("--from-manifest"), "{err}"); let mut dump_manifest = base_args(); dump_manifest.reference = PathBuf::from("/nonexistent/ref.fa"); - dump_manifest.query = PathBuf::from("/nonexistent/qry.fa"); + dump_manifest.query = Some(PathBuf::from("/nonexistent/qry.fa")); dump_manifest.dump_manifest = Some(PathBuf::from("/nonexistent/out.manifest")); let err = prepare(&dump_manifest, &mut phases).err().unwrap(); assert!(err.to_string().contains("--dump-manifest"), "{err}"); diff --git a/src/seed.rs b/src/seed.rs index 554ff46..36fdfbb 100644 --- a/src/seed.rs +++ b/src/seed.rs @@ -465,6 +465,316 @@ pub fn max_seeds(span: u32, shape: &Shape, transitions: bool) -> usize { span as usize * per_pos } +// --------------------------------------------------------------------------- +// R92 PR1: exact per-unit seed-hit estimator (host histograms, no GPU). +// +// `hits(r,q) = Σ_key refcount_r[key] × qrycount_q[key]` over the pipeline's +// compact seeds, equal to the `unit ledger:` hits column EXACTLY (u64). The +// reference counts use the same strided starts as `SeedTable::build`; the +// query counts use the same `kmer_at` validity over the same chunked LASTZ +// interval schedule (`seed::chunks`) on both strands, counting base k-mers +// only. Transition variants are folded into the reference side +// (`fold_transitions`): each valid query window emits +// its base k-mer plus one single-transition variant per transition care +// position (`push_seeds`, identical XOR mapping), so dotting base query +// counts against folded reference weights is exact. + +/// Dense-table guard: one `u32` table is `4^k` entries; `k > 12` needs 1 GiB+ +/// per table and is out of scope (12of19 default). +pub fn estimator_table_size(shape: &Shape) -> Result { + if shape.kmer_size > 12 { + return Err(format!( + "hits-estimate needs k<=12 (dense 4^k table); seed has k={}", + shape.kmer_size + )); + } + Ok(1usize << (2 * shape.kmer_size)) +} + +/// Strided start sequence shared with `SeedTable::build`: positions +/// `start + i*step` for `i in 0..num` over a block of `len` bases. +fn ref_starts(len: usize, shape: &Shape, step: u32) -> (usize, usize) { + let step = step.max(1) as usize; + let offset = (shape.size + 1) % step; + let num = (len + offset).saturating_sub(shape.size) / step; + (step - offset, num) +} + +/// Serial reference histogram over a packed block. +pub fn count_ref_block(buf: &[u8], shape: &Shape, step: u32) -> Result, String> { + let nkeys = estimator_table_size(shape)?; + let (start, num) = ref_starts(buf.len(), shape, step); + let step = step.max(1) as usize; + let mut hist = vec![0u32; nkeys]; + for i in 0..num { + let k = shape.kmer_at(buf, start + i * step); + if k != INVALID_KMER { + hist[k as usize] += 1; + } + } + Ok(hist) +} + +/// Adds `part` into `acc` (the histogram merge; counts commute). +fn merge_hist(acc: &mut [u32], part: &[u32]) { + for (a, b) in acc.iter_mut().zip(part.iter()) { + *a += *b; + } +} + +/// Parallel reference histogram: the step-index range is split into `threads` +/// contiguous pieces, each counted into a private table, then summed. +/// Exact — the same positions, the same `kmer_at`, addition commutes. +pub fn count_ref_block_parallel( + buf: &[u8], + shape: &Shape, + step: u32, + threads: usize, +) -> Result, String> { + let nkeys = estimator_table_size(shape)?; + let (start, num) = ref_starts(buf.len(), shape, step); + let step = step.max(1) as usize; + let threads = threads.max(1).min(num.max(1)); + if threads <= 1 || num < 1 << 16 { + return count_ref_block(buf, shape, step as u32); + } + let per = num.div_ceil(threads); + let parts: Vec> = std::thread::scope(|s| { + let hs: Vec<_> = (0..threads) + .map(|t| { + let (lo, hi) = (t * per, ((t + 1) * per).min(num)); + s.spawn(move || { + let mut h = vec![0u32; nkeys]; + for i in lo..hi { + let k = shape.kmer_at(buf, start + i * step); + if k != INVALID_KMER { + h[k as usize] += 1; + } + } + h + }) + }) + .collect(); + hs.into_iter() + .map(|h| h.join().expect("ref count worker panicked")) + .collect() + }); + let mut out = vec![0u32; nkeys]; + for p in &parts { + merge_hist(&mut out, p); + } + Ok(out) +} + +/// Key-sharded reference histogram (the alternative strategy, test-only): +/// the key space is split into `threads` shards and every worker scans all +/// positions, counting only keys in its shard — one final table, no merge, +/// but `kmer_at` is recomputed per worker. Used by the strategy-comparison +/// test; production uses the private-tables path above. +#[cfg(test)] +pub fn count_ref_block_sharded( + buf: &[u8], + shape: &Shape, + step: u32, + threads: usize, +) -> Result, String> { + let nkeys = estimator_table_size(shape)?; + let (start, num) = ref_starts(buf.len(), shape, step); + let step = step.max(1) as usize; + let threads = threads.max(1).min(nkeys); + if threads <= 1 || num < 1 << 16 { + return count_ref_block(buf, shape, step as u32); + } + let kper = nkeys.div_ceil(threads); + let mut out = vec![0u32; nkeys]; + { + let ptr = out.as_mut_ptr() as usize; + std::thread::scope(|s| { + for t in 0..threads { + let (lo, hi) = (t * kper, ((t + 1) * kper).min(nkeys)); + s.spawn(move || { + for i in 0..num { + let k = shape.kmer_at(buf, start + i * step); + if k != INVALID_KMER && (k as usize) >= lo && (k as usize) < hi { + // SAFETY: shard `[lo,hi)` is disjoint across workers. + unsafe { *(ptr as *mut u32).add(k as usize) += 1 }; + } + } + }); + } + }); + } + Ok(out) +} + +/// Base-k-mer query histogram over one packed query bin: `fwd[..block_len]` +/// plus its reverse complement `rc`, over the pipeline's LASTZ intervals on +/// the strands `plus`/`minus` select. Counts base k-mers only (variants are +/// folded on the reference side). `stride` S counts window starts `j % S == 0` +/// only (S=1 is exact); `lastz_interval` is the pipeline's `-I` tiling and +/// `wga_chunk` its `-C` chunk size. Interval ranges are split across `threads` +/// workers with private tables, then merged. +/// +/// The window starts are exactly the ones `seed_and_filter_all` batches: each +/// interval is cut with [`chunks`] (whose exclusive upper bound includes the +/// interval end only when the chunk walk does not land on it — an interval +/// whose length is a multiple of `wga_chunk` leaves its end to the next +/// interval), and the minus strand chunks its RC mirror +/// `(qbl - end, qbl - start)`, whose own chunk alignment decides the mirror's +/// duplicate windows. Counting raw interval ends as inclusive, as the first +/// version of this estimator did, double-counts every shared endpoint. +pub fn count_query_block( + fwd: &[u8], + rc: &[u8], + block_len: usize, + shape: &Shape, + plus: bool, + minus: bool, + lastz_interval: u32, + wga_chunk: u32, + stride: u32, + threads: usize, +) -> Result, String> { + let nkeys = estimator_table_size(shape)?; + if block_len <= shape.size { + return Ok(vec![0u32; nkeys]); + } + let stride = stride.max(1); + let qbl = (block_len - shape.size) as u32; + let ivs = crate::sequence::intervals(block_len, shape.size, lastz_interval); + // (is_rc, lo, hi-exclusive) window-start ranges, exactly the pipeline's. + let mut ranges: Vec<(bool, u32, u32)> = Vec::new(); + if plus { + for &(s, e) in &ivs { + ranges.extend( + chunks(s, e, wga_chunk) + .into_iter() + .map(|(lo, hi)| (false, lo, hi)), + ); + } + } + if minus { + for &(s, e) in &ivs { + let (lo, hi) = (qbl - e, qbl - s); + ranges.extend( + chunks(lo, hi, wga_chunk) + .into_iter() + .map(|(l, h)| (true, l, h)), + ); + } + } + let threads = threads.max(1).min(ranges.len().max(1)); + if threads <= 1 { + let mut hist = vec![0u32; nkeys]; + count_query_ranges(&mut hist, fwd, rc, &ranges, shape, stride); + return Ok(hist); + } + let chunk = ranges.len().div_ceil(threads); + let parts: Vec> = std::thread::scope(|s| { + let hs: Vec<_> = ranges + .chunks(chunk) + .map(|ch| { + s.spawn(move || { + let mut h = vec![0u32; nkeys]; + count_query_ranges(&mut h, fwd, rc, ch, shape, stride); + h + }) + }) + .collect(); + hs.into_iter() + .map(|h| h.join().expect("query count worker panicked")) + .collect() + }); + let mut out = vec![0u32; nkeys]; + for p in &parts { + merge_hist(&mut out, p); + } + Ok(out) +} + +/// Counts base k-mers over `ranges` (exclusive upper bounds, as [`chunks`] +/// returns them), skipping window starts `j % stride != 0`. +fn count_query_ranges( + hist: &mut [u32], + fwd: &[u8], + rc: &[u8], + ranges: &[(bool, u32, u32)], + shape: &Shape, + stride: u32, +) { + for &(is_rc, lo, hi) in ranges { + let seq = if is_rc { rc } else { fwd }; + for j in lo..hi { + if j % stride != 0 { + continue; + } + let k = shape.kmer_at(seq, j as usize); + if k != INVALID_KMER { + hist[k as usize] += 1; + } + } + } +} + +/// Folds transition variants into reference weights (exact): each valid query +/// window with base k-mer K scores `ref[K] + Σ_t ref[K ^ (2 << 2t)]` over the +/// transition care positions — the same XOR mapping `push_seeds` emits — so +/// the dot product over base query counts reproduces the pipeline's compact +/// seed stream. `transitions=false` (`--notransition`) is the identity. u64: +/// a 3 Gbp bin holds < 2^32 positions per key, but the 13-term sum can exceed +/// u32 (up to 13× the bin length in the adversarial case). +pub fn fold_transitions(ref_hist: &[u32], shape: &Shape, transitions: bool) -> Vec { + let mut out: Vec = ref_hist.iter().map(|&c| c as u64).collect(); + if !transitions { + return out; + } + let n = out.len(); + for t in 0..shape.kmer_size { + if !shape.transition[t] { + continue; + } + let flip = 2u32 << (2 * t); + for k in 0..n { + out[k] += ref_hist[(k as u32 ^ flip) as usize] as u64; + } + } + out +} + +/// Row-major `r*Q+q` exact hit counts: each reference histogram is folded +/// once, then dotted against every query histogram (base counts). One scoped +/// worker per reference bin (R is small); the histogram passes above carry +/// the thread budget. +pub fn unit_hits( + ref_hists: &[Vec], + qry_hists: &[Vec], + shape: &Shape, + transitions: bool, +) -> Vec { + let q = qry_hists.len(); + if ref_hists.is_empty() || q == 0 { + return Vec::new(); + } + let rows: Vec> = std::thread::scope(|s| { + let hs: Vec<_> = ref_hists + .iter() + .map(|rh| { + s.spawn(move || { + let w = fold_transitions(rh, shape, transitions); + qry_hists + .iter() + .map(|qh| w.iter().zip(qh.iter()).map(|(&a, &b)| a * b as u64).sum()) + .collect::>() + }) + }) + .collect(); + hs.into_iter() + .map(|h| h.join().expect("dot worker panicked")) + .collect() + }); + rows.into_iter().flatten().collect() +} + /// The `[i, e)` chunk bounds `seeder.cpp` walks for one interval. `end` is /// inclusive there, hence the `+ 1`. pub fn chunks(start: u32, end: u32, wga_chunk: u32) -> Vec<(u32, u32)> { From 34d30f650b8571b4300a4ce26d03ed77f658e2ae Mon Sep 17 00:00:00 2001 From: alejandrogzi Date: Fri, 11 Sep 2026 10:21:29 +0200 Subject: [PATCH 2/3] =?UTF-8?q?release:=20v0.0.6=20=E2=80=94=20--query-lis?= =?UTF-8?q?t=20batch=20mode,=20hits-estimate?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- Cargo.lock | 2 +- Cargo.toml | 2 +- assets/changelog/changelog.md | 20 ++++++++++++++++++++ 3 files changed, 22 insertions(+), 2 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index e8e2b20..644d958 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -385,7 +385,7 @@ checksum = "2304e00983f87ffb38b55b444b5e3b60a884b5d30c0fca7d82fe33449bbe55ea" [[package]] name = "hspz" -version = "0.0.5" +version = "0.0.6" dependencies = [ "clap", "cuda-core", diff --git a/Cargo.toml b/Cargo.toml index 54ace64..761a9d8 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,7 +1,7 @@ [package] name = "hspz" authors = ["Alejandro Gonzales-Irribarren "] -version = "0.0.5" +version = "0.0.6" edition = "2024" license = "GPL-3.0" description = "GPU-accelerated high-scoring ungapped alignment pair backend" diff --git a/assets/changelog/changelog.md b/assets/changelog/changelog.md index d11be47..f856e5b 100644 --- a/assets/changelog/changelog.md +++ b/assets/changelog/changelog.md @@ -36,6 +36,26 @@ All notable changes to `hspZ` are documented here, newest first. +## [0.0.6] — 2026-09-11 + +One reference, many queries, and a CPU-only seed-hit estimator. + +- **One reference, many queries: `--query-list`.** `hspz run --reference ref.fa --query-list queries.txt --output OUT` + runs every listed query FASTA in one process: each reference bin is built and uploaded once instead of once per + query, while every query keeps its own plan, cap, seed batches, dedup scopes and `-D` history, so its output under + `OUT/000001/`, `OUT/000002/`, … is byte-identical to a standalone run (verified file for file on 100 queries). + `-Z` writes one archive per job; `OUT/queries.tsv` maps job → path, size, sha256, blocks, units, HSPs, completion. + Measured on an L4 against hg38: 100 queries of 5 Mbp took 10,400–11,600 s as separate runs and 1,890–1,900 s + batched (−82%); three 5 Mbp queries −62%; whole-chromosome queries (61–121 Mbp) −2.6% (GPU-bound); whole-genome + queries gain under 1%. This release: W=1 (`--gpus 1`), all jobs must plan to identical reference bins, `-B 0`, + `--kegalign-bins` and `--from-manifest` are rejected in batch mode, a failing job aborts the batch (partial outputs + may remain), dumps become per-job siblings (`out.manifest.000001`). Bin-major order delays the first result. +- **`hspz hits-estimate` (CPU only).** Prints, before any GPU work, the exact number of seed hits each work unit + (reference bin × query block) will produce, from dense k-mer histograms with the transition variants folded in; + matches the `--time` ledger exactly (whole genome: 42 units, 6.12e12 hits). GPU time per unit is ~0.37 ns per hit on + a 4090, so the rows are the unit costs of a multi-GPU schedule. 56–62 s on 32 cores for hg38 × mm39 (10 GiB RSS); + `--stride N` samples the query side. A planning/diagnostic tool; `run` does not call it. + ## [0.0.5] — 2026-09-10 Unit-level static partitioning behind a device-class guard, with a per-unit ledger From 0f2d9fc9e84d6053203891896ee3fbef0d09b30e Mon Sep 17 00:00:00 2001 From: alejandrogzi Date: Fri, 11 Sep 2026 10:23:20 +0200 Subject: [PATCH 3/3] docs: --query-list batch mode, hits-estimate subcommand, regenerated CLI reference --- docs/scripts/generate-cli-reference.sh | 2 +- docs/src/generated/cli-reference.txt | 61 ++++++++++++++++++++++---- docs/src/pages/docs/cli.astro | 2 +- docs/src/pages/docs/gpus.mdx | 25 +++++++++++ 4 files changed, 80 insertions(+), 10 deletions(-) diff --git a/docs/scripts/generate-cli-reference.sh b/docs/scripts/generate-cli-reference.sh index 59a86f7..f437777 100644 --- a/docs/scripts/generate-cli-reference.sh +++ b/docs/scripts/generate-cli-reference.sh @@ -18,7 +18,7 @@ emit() { { emit '' - for command in run benchmark compare; do + for command in run benchmark compare hits-estimate; do printf '\n%s\n\n' "$divider" emit "$command" done diff --git a/docs/src/generated/cli-reference.txt b/docs/src/generated/cli-reference.txt index 62de37a..9f0708f 100644 --- a/docs/src/generated/cli-reference.txt +++ b/docs/src/generated/cli-reference.txt @@ -5,10 +5,11 @@ GPU-accelerated high-scoring ungapped alignment pair backend Usage: hspZ Commands: - run Seed, extend and filter one reference/query pair - benchmark Cold/warm benchmark: the optimization oracle (PLAN.md Milestone 3) - compare Run the C++ CUDA reference and this implementation back to back and compare runtime and HSP output - help Print this message or the help of the given subcommand(s) + run Seed, extend and filter one reference/query pair + benchmark Cold/warm benchmark: the optimization oracle (PLAN.md Milestone 3) + compare Run the C++ CUDA reference and this implementation back to back and compare runtime and HSP output + hits-estimate Exact per-unit seed-hit estimates before any GPU work (CPU-only) + help Print this message or the help of the given subcommand(s) Options: -h, --help Print help @@ -20,13 +21,15 @@ $ hspZ run --help Seed, extend and filter one reference/query pair -Usage: hspZ run [OPTIONS] --reference --query +Usage: hspZ run [OPTIONS] --reference Options: -r, --reference -q, --query + --query-list + One query FASTA path per line (blank lines and `#` comments ignored). Bin-major batch: one reference × many queries, W=1 only. All-or-nothing: a failing job aborts the batch, so partial per-job outputs (and a truncated `-Z` archive) may remain on disk -o, --output Directory to write `tmp.block.r.{plus,minus}.segments` into [default: .] -S, --strand @@ -96,7 +99,7 @@ Options: -D, --diagonal-partition Split each output file along its diagonal, KegAlign-style, emitting `*.split1.segments` onward instead of one file. Partitions the HSP structs directly — no unsplit file is written first -Z, --tarball [] - Write output into one `.tar.gz` instead of a directory. Without a path, `.tar.gz` is used. Archived directly from the formatted bytes + Write output into one `.tar.gz` instead of a directory. Without a path, `.tar.gz` is used. Archived directly from the formatted bytes. In batch mode (`--query-list`) one archive per job is written (`OUT/000001.tar.gz` …) and a custom path value is ignored (a note is printed) -h, --help Print help @@ -106,13 +109,15 @@ $ hspZ benchmark --help Cold/warm benchmark: the optimization oracle (PLAN.md Milestone 3) -Usage: hspZ benchmark [OPTIONS] --reference --query +Usage: hspZ benchmark [OPTIONS] --reference Options: -r, --reference -q, --query + --query-list + One query FASTA path per line (blank lines and `#` comments ignored). Bin-major batch: one reference × many queries, W=1 only. All-or-nothing: a failing job aborts the batch, so partial per-job outputs (and a truncated `-Z` archive) may remain on disk -o, --output Directory to write `tmp.block.r.{plus,minus}.segments` into [default: .] -S, --strand @@ -182,7 +187,7 @@ Options: -D, --diagonal-partition Split each output file along its diagonal, KegAlign-style, emitting `*.split1.segments` onward instead of one file. Partitions the HSP structs directly — no unsplit file is written first -Z, --tarball [] - Write output into one `.tar.gz` instead of a directory. Without a path, `.tar.gz` is used. Archived directly from the formatted bytes + Write output into one `.tar.gz` instead of a directory. Without a path, `.tar.gz` is used. Archived directly from the formatted bytes. In batch mode (`--query-list`) one archive per job is written (`OUT/000001.tar.gz` …) and a custom path value is ignored (a note is printed) -i, --iterations Timed Seed + Filter iterations after the warmup [default: 20] -w, --warmup @@ -231,3 +236,43 @@ Options: -h, --help Print help + +------------------------------------------------------------------------------ + +$ hspZ hits-estimate --help + +Exact per-unit seed-hit estimates before any GPU work (CPU-only) + +Usage: hspZ hits-estimate [OPTIONS] --reference --query + +Options: + -r, --reference + + -q, --query + + -B, --seq-block-size + Target bin size in bases, as in `run` (`-B 0` resolves to the default) [default: 500000000] + --query-block-size + Bin target for the query side; defaults to `--seq-block-size` + --kegalign-bins + Bin records the way KegAlign does (sequential fill, not LPT) + --seed + 12of19, 14of22, or an arbitrary pattern of 1s, 0s and Ts (k<=12) [default: 12of19] + -e, --step + [default: 1] + --notransition + Don't allow one transition in a seed hit + -S, --strand + plus / minus / both [default: both] + -t, --threads + Worker threads for histogram passes; 0 uses available parallelism [default: 0] + -I, --lastz-interval-size + LASTZ interval size, as in `run` (`-I`): interval boundaries decide which shared window starts the pipeline seeds once and which twice [default: 10000000] + -C, --wga-chunk-size + Query seed chunk size, as in `run` (`-C`): the interval schedule (`seed::chunks`) is part of the count, so this must match `run` [default: 250000] + --stride + Count every S-th query window start and scale by S (approximation) [default: 1] + --dump-plan + Write the plan's bin membership to this file and continue + -h, --help + Print help diff --git a/docs/src/pages/docs/cli.astro b/docs/src/pages/docs/cli.astro index 097acad..31187aa 100644 --- a/docs/src/pages/docs/cli.astro +++ b/docs/src/pages/docs/cli.astro @@ -4,6 +4,6 @@ import reference from "../../generated/cli-reference.txt?raw" --- -

This reference covers `hspZ`, `run`, `benchmark`, and `compare`. CI regenerates it from the release binary and fails if the snapshot drifts.

+

This reference covers `hspZ`, `run`, `benchmark`, `compare`, and `hits-estimate`. CI regenerates it from the release binary and fails if the snapshot drifts.

{reference}
 diff --git a/docs/src/pages/docs/gpus.mdx b/docs/src/pages/docs/gpus.mdx index 0b4d57c..3d06a4e 100644 --- a/docs/src/pages/docs/gpus.mdx +++ b/docs/src/pages/docs/gpus.mdx @@ -38,6 +38,31 @@ Every knob is exact: every mode produces the same HSP set and the same bytes. | `HSPZ_CHUNK_WALK=full` | sparse | Restores the pre-0.0.4 chunk walk that copied the whole cumulative hit array per seed batch (rollback only). | | `HSPZ_UNIT_PARTITION=auto\|0\|1` | `auto` | Splits the frozen plan's work units across workers by count quotas when every device the run uses is the same class and each worker has its own device (W=1, time-sliced workers and mixed classes keep whole-bin ownership); output is unchanged. −11.9% seed-and-filter wall on 2× RTX 4090 (1,400 → 1,233 s) and −5.1% on 4× RTX 4090 (718 → 681 s) at `MAX_HITS=16,711,680`. | +## Many queries against one reference + +`--query-list` runs every listed query FASTA in one process: + +```bash +hspZ run \ + --reference ref.2bit \ + --query-list queries.txt \ + --output OUT \ + --gpus 1 +``` + +`queries.txt` holds one query path per line (blank lines and `#` comments are +ignored). Each reference bin is built and uploaded once instead of once per +query, while every query keeps its own plan, cap, seed batches, dedup scopes +and `-D` history, so its output under `OUT/000001/`, `OUT/000002/`, … is +byte-identical to a standalone run. `-Z` writes one archive per job; +`OUT/queries.tsv` maps job → path, size, sha256, blocks, units, HSPs and +completion. A failing job aborts the batch (partial outputs may remain). + +In 0.0.6 batch mode requires W=1 (`--gpus 1`), all jobs must plan to identical +reference bins, `-B 0` is rejected as a layout, and `--kegalign-bins` and +`--from-manifest` are rejected. Bin-major order delays the first result: check +`OUT/queries.tsv` rather than waiting on the first job directory. + ## Choosing `-B` The default 500 Mbp preserves the KegAlign-matched one-GPU digest; hspZ never