v0.0.6 — one reference × many queries (--query-list), hits-estimate

Cargo.toml

@@ -1,7 +1,7 @@
 [package]
 name = "hspz"
 authors = ["Alejandro Gonzales-Irribarren <alejandrxgzi@gmail.com>"]
-version = "0.0.5"
+version = "0.0.6"
 edition = "2024"
 license = "GPL-3.0"
 description = "GPU-accelerated high-scoring ungapped alignment pair backend"

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

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

docs/src/generated/cli-reference.txt

@@ -5,10 +5,11 @@ GPU-accelerated high-scoring ungapped alignment pair backend
 Usage: hspZ <COMMAND>
 
 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 <REFERENCE> --query <QUERY>
+Usage: hspZ run [OPTIONS] --reference <REFERENCE>
 
 Options:
   -r, --reference <REFERENCE>
 
   -q, --query <QUERY>
 
+      --query-list <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 <OUTPUT>
           Directory to write `tmp<n>.block<q>.r<r>.{plus,minus}.segments` into [default: .]
   -S, --strand <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 [<TARBALL>]
-          Write output into one `.tar.gz` instead of a directory. Without a path, `<output>.tar.gz` is used. Archived directly from the formatted bytes
+          Write output into one `.tar.gz` instead of a directory. Without a path, `<output>.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 <REFERENCE> --query <QUERY>
+Usage: hspZ benchmark [OPTIONS] --reference <REFERENCE>
 
 Options:
   -r, --reference <REFERENCE>
 
   -q, --query <QUERY>
 
+      --query-list <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 <OUTPUT>
           Directory to write `tmp<n>.block<q>.r<r>.{plus,minus}.segments` into [default: .]
   -S, --strand <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 [<TARBALL>]
-          Write output into one `.tar.gz` instead of a directory. Without a path, `<output>.tar.gz` is used. Archived directly from the formatted bytes
+          Write output into one `.tar.gz` instead of a directory. Without a path, `<output>.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 <ITERATIONS>
           Timed Seed + Filter iterations after the warmup [default: 20]
   -w, --warmup <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 <REFERENCE> --query <QUERY>
+
+Options:
+  -r, --reference <REFERENCE>
+
+  -q, --query <QUERY>
+
+  -B, --seq-block-size <SEQ_BLOCK_SIZE>
+          Target bin size in bases, as in `run` (`-B 0` resolves to the default) [default: 500000000]
+      --query-block-size <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 <SEED>
+          12of19, 14of22, or an arbitrary pattern of 1s, 0s and Ts (k<=12) [default: 12of19]
+  -e, --step <STEP>
+          [default: 1]
+      --notransition
+          Don't allow one transition in a seed hit
+  -S, --strand <STRAND>
+          plus / minus / both [default: both]
+  -t, --threads <THREADS>
+          Worker threads for histogram passes; 0 uses available parallelism [default: 0]
+  -I, --lastz-interval-size <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 <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 <STRIDE>
+          Count every S-th query window start and scale by S (approximation) [default: 1]
+      --dump-plan <DUMP_PLAN>
+          Write the plan's bin membership to this file and continue
+  -h, --help
+          Print help

docs/src/pages/docs/cli.astro

@@ -4,6 +4,6 @@ import reference from "../../generated/cli-reference.txt?raw"
 ---
 
 <DocsLayout title="CLI reference" description="A checked-in snapshot generated directly from the production Clap definitions.">
-  <p>This reference covers `hspZ`, `run`, `benchmark`, and `compare`. CI regenerates it from the release binary and fails if the snapshot drifts.</p>
+  <p>This reference covers `hspZ`, `run`, `benchmark`, `compare`, and `hits-estimate`. CI regenerates it from the release binary and fails if the snapshot drifts.</p>
   <pre><code>{reference}</code></pre>
 </DocsLayout>

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

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<PathBuf>,
+    /// 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<PathBuf>,
     /// Directory to write `tmp<n>.block<q>.r<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,
     /// `<output>.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<PathBuf>,
 }
@@ -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<u32>,
+    /// 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<PathBuf>,
+}
+
 #[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"
+        );
     }
 }

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(),