obikmer/docmd/architecture/query.md

Query system

Goal

Given a set of query sequences, determine for each sequence how many of its k-mers are found in the index and, for each indexed genome, how many k-mers match. The query system is the foundation for read classification and sequence-to-genome mapping.

---

Input

• Query sequences in FASTA or FASTQ format (gzip supported, streaming stdin supported). GenBank flat files are not supported at query time (only at index time).
• Sequences shorter than k bases are silently skipped.
• Non-ACGT characters are handled by the superkmer decomposition layer: they act as hard breaks, producing shorter superkmers (identical to the behaviour at indexing time).

---

Algorithm

The query follows the same superkmer-based partitioning strategy used at indexing time. Everything below happens inside process_chunk (query.rs); there is no separate per-stage function — dereplication, MPHF/matrix lookup, and the Findere window are three phases of one function, operating on a single flat allocation per chunk.

for each chunk of sequences (parallel workers via obipipeline, one call to process_chunk):
    build QueryBatch (QueryBatch::from_records):
        decompose all sequences into superkmers (SuperKmerIter)
        deduplicate identical superkmers into
            map: HashMap<RoutableSuperKmer, Vec<SKDesc>>   ← SKDesc = (seq_idx, kmer_offset)
    allocate KmerResults (KmerResults::new): one flat allocation for the whole chunk —
        data: Vec<u32>  sized total_smers × n_genomes  (row-major, zero-initialised)
        in_index: Vec<bool>  sized total_smers
    split the batch's unique superkmers by partition, via minimiser hash
    for each partition p:
        query_partition_with(p, superkmers_routed_to_p, on_hit):
            → load QueryLayer(s) for p
            → for each s-mer of each superkmer: try each layer's MphfLayer::find in turn,
              stop at the first hit (a k-mer belongs to at most one layer)
            → on_hit(sk_idx, kmer_idx, row): broadcast `row` into KmerResults, once per
              SKDesc referencing this superkmer (all (seq_idx, position) occurrences)
    for each sequence, for each genome (inlined sliding window, no separate function):
        monotone-deque scan over s-mer positions → win_min[pos][g] = min over the
        z-window [pos, pos+z), "not in index" treated as 0
    for each sequence, for each output position (0..n_kmers_out):
        accumulate confirmed k_user-mer results into acc (kmer_count / kmer_missing /
        per-genome totals) and, if --detail, into cov
    emit annotated sequences (emit_batch)

Superkmers that appear more than once in the batch (same sequence or across sequences) are deduplicated: each unique RoutableSuperKmer is queried once per partition, and the result is broadcast to every SKDesc entry that references it. See Future work for a discussion of pushing this dereplication down to k-mer granularity.

Findere requires full-sequence aggregation. The sliding window is run once per (sequence, genome) pair on the complete s-mer result range, after all partitions have contributed to KmerResults. Applying it per superkmer would produce false negatives at superkmer boundaries, where the z-window spans two superkmers.

Batches are processed in parallel via obipipeline workers; the --threads flag controls the number of worker threads.

---

Findere z-window filter

For approximate index modes, the index physically stores s-mers of size s = k_user − z + 1; idx.kmer_size() (bound to k in process_chunk) is this physically-indexed s-mer size, so decomposing the query at k naturally produces s-mer results.

The z-window aggregation is inlined per (sequence, genome) in process_chunk — a monotone-deque sliding-window minimum, not a separate function operating on an intermediate Option<Box<[u32]>> vector:

for each sequence, for each genome g independently:
    dq: VecDeque<(pos, value)>          // reused across all (sequence, genome) pairs
    for i in 0..n_smers:
        v = if results.is_in_index(seq, i) { results.val(seq, i, g) } else { 0 }
        evict dq entries that left the window [i-z+1, i]
        maintain dq monotone non-decreasing (pop back while back.value >= v), push (i, v)
        if i + 1 >= z:
            pos = i + 1 - z
            win_min[pos][g] = dq.front().value     // minimum over the z-window

Input length n_smers, output length n_kmers_out = n_smers − z + 1 (k_user-mer positions). The scan is O(n_smers) per genome (each position enters and leaves the deque once).

Output position pos is confirmed for genome g iff win_min[pos][g] > 0 — equivalent to "all z consecutive s-mer values in the window are nonzero for g", since a monotone-deque minimum is 0 as soon as any element in the window is 0.

The value reported per confirmed position is the window minimum, not the leftmost s-mer's raw value. For presence indexes (0/1 values) this is equivalent to a logical AND either way. For count indexes it is not: the accumulated count for genome g at position pos is min(results[pos..pos+z][g]), the weakest link across the window — not results[pos][g].

kmer_missing bookkeeping is separate from the window value and does look specifically at the leftmost s-mer: a position with no genome confirmed (win_min[pos][*] == 0 for all g) counts as kmer_missing iff the leftmost s-mer of that window, results.is_in_index(seq, pos), is absent from the index entirely (see kmer_missing semantics).

Short sequences: when the s-mer count is less than z, n_kmers_out is 0 (process_chunk skips the sequence: if out_n == 0 { continue; }) — no k_user-mer is emitted for that sequence.

Exact indexes: z = 1, the window degenerates to a single element — the "minimum" is just that element's value, and pos == i, i.e. a passthrough.

Effective z at query time

effective_z is resolved at the start of run():

let effective_z = args.findere_z.unwrap_or_else(|| match idx.meta().config.evidence {
    IndexMode::Approx { z, .. } | IndexMode::Hybrid { z, .. } => z as usize,
    IndexMode::Exact => 1,
});

The -z CLI option overrides the index metadata value. A higher z increases stringency (lower FP, some true positives may be discarded at sequence ends); a lower z increases sensitivity.

---

Layer lookup: MphfLayer::find

MphfLayer::open(dir, mode: &IndexMode) receives the mode from PartitionMeta — no per-layer file is read. The caller (QueryLayer) never chooses the dispatch path: it is fixed at open time by LayerEvidence. See obilayeredmap for the full find / find_strict API.

QueryLayer variant selection

QueryLayer::open (obikpartitionner/src/query_layer.rs:28-45) only ever returns two variants — Presence or Count, checked in this order:

Order	Condition	Variant	Data returned per k-mer
1	with_counts=true and counts/ exists	Count	raw count per genome
2	(else) presence/ exists, or counts/ doesn't exist at all	Presence	see below
3	(else — counts/ exists, presence/ doesn't, with_counts=false)	Count	counts used as-is

There is no QueryLayer::SetOnly variant. The "no on-disk matrix at all" case is handled one level down: Presence wraps PersistentBitMatrix, whose own open() (obicompactvec/src/bitmatrix.rs:260-288) auto-detects among three internal representations — Packed (presence/matrix.pbmx), Columnar (presence/meta.json), or Implicit { n_rows, n_cols } when neither file exists (built from layer_meta.json, fill_row returning all-1s without touching disk). This is where "1 for every genome" actually happens — not at the QueryLayer level.

Worth double-checking, not confirmed as a bug: PersistentBitMatrix::open's Implicit branch constructs Implicit { n_rows: meta.n, n_cols: 1 } — n_cols is hardcoded to 1, not to the layer's actual n_genomes. fill_row for Implicit only writes buf[..1], leaving the rest of a longer n_genomes-sized buffer untouched (zeroed by the caller beforehand). If this path is ever reached for a layer covering more than one genome, only genome index 0 would read as present. Whether that's reachable in practice (layers might always be single-genome when they fall back to Implicit) wasn't verified here — flagging for follow-up, not fixing.

---

Presence / count mode at query time

The --force-presence flag and --presence-threshold control how per-genome values are accumulated, independently of what the index stores:

genome_totals[g] += if presence { u32::from(v >= threshold) } else { v }

presence is true when --force-presence is set or when the index has no counts (!with_counts). The default presence_threshold is 1, so any nonzero count counts as a match.

---

Coverage vectors (--detail)

When --detail is requested, a 3-D accumulator cov[seq_idx][genome][kmer_pos] is allocated after all partitions are queried, with dimensions derived from n_kmers_out = n_smers − z + 1 (k_user-mer positions, not s-mer positions):

cov[seq_idx][g][pos] += contribution
where pos is the k_user-mer index in the filtered (post-Findere) vector

Coverage reflects confirmed k_user-mers only. The vectors are emitted in the JSON annotation under the key "coverage".

---

kmer_missing semantics

kmer_missing counts k_user-mer positions where the leftmost s-mer of the window (results.is_in_index(seq_idx, pos), KmerResults) is false — i.e. absent from the index entirely. K-mers where the z-window fails because a later s-mer is absent or zero (but the leftmost one is present) are not counted as missing — the leftmost s-mer being present is used as proxy for index membership.

---

Output format

Output sequences are written in OBITools4 format: the original sequence with a JSON annotation map in the title line.

>read_id {"kmer_count":59,"kmer_strict_matches":{"genome_a":42,"genome_b":7}}
ATCGATCG...

With --detail:

>read_id {"kmer_count":59,"kmer_strict_matches":{...},"coverage":{"genome_a":[0,1,2,...],...}}
ATCGATCG...

Genome keys follow the iteration order of meta.genomes.

---

Annotation schema

Key	Type	Condition	Semantics
kmer_count	int	always	k-mers confirmed (post-Findere) with at least one genome match
kmer_missing	int	--count-missing	k-mers absent from the index entirely (leftmost s-mer of the window not found)
kmer_strict_matches	object	always	per-genome accumulated value, non-zero entries only (label → count or 0/1)
coverage	object	--detail	per-genome array of per-position contributions (label → [u32])

kmer_count + kmer_missing ≤ total k_user-mers in the sequence. The gap corresponds to k_user-mers whose z-window was not fully confirmed (at least one s-mer absent or zero for all genomes) but whose first s-mer was present in the index.

---

CLI

obikmer query <index> [--detail] [--mismatch] [--count-missing]
              [--force-presence] [--presence-threshold <n>]
              [-z <z>] [-T <threads>] [--chunk-size <MiB>]
              <query.fa> [<query2.fa> ...]

Option	Default	Semantics
-z / --findere-z	from index metadata	Override Findere z parameter
--detail	off	Emit per-position coverage vectors in JSON
--count-missing	off	Add kmer_missing field to JSON
--force-presence	off	Report 0/1 per genome regardless of index counts
--presence-threshold	1	Minimum count to declare genome present
-T / --threads	all CPUs	Worker threads
--chunk-size	auto (from available RAM and thread count)	I/O chunk size in MiB — see Future work, point 3 for why the auto-sizing formula currently under-estimates memory on indexes with many genomes

--mismatch is accepted but currently ignored with a warning on stderr.

---

Future work

• --mismatch: 1-mismatch approximate matching — generate 3·k single-substitution variants per k-mer, look each up independently.
• Read classification (--classify): assign each read to the genome with the highest match score.
• Whitelist / blacklist filtering: threshold-based accept/reject on per-genome match scores.

Throughput & parallelism — identified potential (not yet implemented)

Observed on a 192-core (8×24 NUMA) machine: query uses ~10 cores or fewer, and the default chunk size gets the process OOM-killed. Root causes and candidate fixes, in dependency order:

1. Single-threaded I/O source (main core-utilization bottleneck). run() builds all_chunks via paths.into_iter().flat_map(read_sequence_chunks_sized(...)) and passes it directly as the input iterator to pipe.apply(). In obipipeline::Pipe::apply (scheduler.rs), input.next() is called exclusively from the dedicated source thread — so file opening, decompression, and FASTA/FASTQ chunk-boundary parsing for all input files run serially in one thread, regardless of --threads. Compare with steps::scatter (used by index) and cmd/superkmer.rs: there, file opening + streaming is itself a Flat pipeline stage (||?), executed across the n_workers pool, with obipipeline::throttle(paths, max_open) bounding concurrently-open files in the source thread. That pattern parallelises I/O across files (and NUMA nodes); query.rs cannot. Fix direction: restructure query's pipe with an initial Flat stage analogous to scatter's, opening/chunking files across workers instead of in flat_map.

2. Gzip decompression is inherently single-threaded per file. niffler/flate2 (used by xopen) do standard DEFLATE, which has no parallel-decodable structure for an arbitrary stream. Fix (1) parallelises across files but not within one large gzip file. Parking a possible fix (rapidgzip-rs) is tracked in chunkreader.md.

3. Chunk-size memory formula ignores n_genomes. chunk_bytes = available_memory_bytes() / (n_workers * 16) (query.rs:407-414) assumes a fixed ~8–16× overhead per raw input byte. But KmerResults::new (query.rs:165-179) allocates data: Vec<u32> sized total_kmers_in_chunk × n_genomes — dense, for every k-mer position in the chunk, hit or not — plus win_min and (with --detail) cov, same scaling. Real per-chunk memory is O(n_genomes), not constant; the formula doesn't know n_genomes at all. This is the direct cause of the OOM kill on indexes with many reference genomes.

4. MPHF lookup and matrix-row fetch are fused, not staged. QueryLayer::find_into (obikpartitionner/src/query_layer.rs:48-67) does the MPHF find and the fill_row matrix read in one call per k-mer, inside a single-threaded loop (query_partition_with). There is no separation between "is this k-mer indexed" (cheap, O(1), independent of n_genomes) and "what are its per-genome values" (the expensive, n_genomes-scaling part).

5. Dereplication should happen at k-mer granularity, directly — not via an intermediate superkmer-level dedup. QueryBatch::from_records currently dereplicates at the superkmer level (HashMap<RoutableSuperKmer, Vec<SKDesc>>, query.rs:112). This misses redundancy between k-mers shared by different superkmers (read overlaps, repeats, a SNP splitting an otherwise-identical run). Superkmer construction (SuperKmerIter) stays mandatory — it is the mechanism that computes minimizers/partition routing, not an optional dedup layer — but the dedup structure built on top of it should key directly on CanonicalKmer, in the same pass: HashMap<CanonicalKmer, Vec<(seq_idx, pos)>>. This also means the MPHF find itself runs once per distinct k-mer instead of once per occurrence — a win independent of the matrix-fetch cost below.

6. Stage 1 output: bucket confirmed hits by layer, keyed by MPHF slot. For each unique canonical k-mer, MPHF lookup across a partition's layers stops at the first match (query_partition_with:105-111) — a k-mer belongs to at most one layer. So stage 1's output can be reshaped directly into:

HashMap<layer_idx, HashMap<slot, Vec<(seq_idx, pos)>>>

replacing the CanonicalKmer key by the resolved slot (compact integer, and exactly what stage 2 needs to address the matrix). K-mers matching no layer simply have no entry here (they still count toward in_index/kmer_missing bookkeeping, which stays O(1) per position, independent of n_genomes).

7. Partition-level parallelism is currently absent — and a NUMA-aware mechanism for exactly this already exists, unused, in obikindex. process_chunk's partition loop (query.rs:250-278, for (part_idx, part_sks) in by_part.iter().enumerate()) processes every partition of a chunk sequentially on the single worker thread that owns that chunk. This is a parallelism axis on its own, independent of the column question below. More importantly: docmd/architecture/numa_partition_runner.md and numa_worker_pools.md document PartitionRunner (obikindex/src/numa.rs), already implemented and already used by merge.rs, index.rs (build_layers), select.rs, reindex.rs, rebuild.rs — one controller thread per NUMA node, a Rayon pool pinned to that node's CPUs (hwlocality, numa feature, default-on in obikindex/Cargo.toml), adaptive worker activation driven by both a CPU-efficiency signal and an I/O-throughput signal (CpuSample/IoSample, /proc/self/io on Linux). It exists precisely because a naive into_par_iter() on the global Rayon pool measurably degrades ×60 on this codebase's own 192-core/8-NUMA reference machine (numa_worker_pools.md, § Problem) once workers contend for cross-socket memory bandwidth on shared mmap'd/hashed structures — exactly the shape of the matrix-column scan in point 8 below. obikmer already depends on obikindex (obikmer/Cargo.toml, for KmerIndex), so PartitionRunner is directly reachable from cmd/query.rs — no new dependency. Both the partition-level loop and (see point 8) the genome-column scan should be driven through it rather than through ad-hoc rayon::into_par_iter(), to avoid reproducing the already-measured-and-fixed contention problem. Also relevant: the "CPU-only signal stalls on I/O-bound stages" issue documented for pack_matrices (mmap-heavy, page-fault-bound) applies just as much to a column-major mmap scan over persistent matrices — reuse the existing dual CPU/IO activation signal rather than re-deriving one.

8. Stage 2: column-major matrix fetch, parallel across genome columns — via PartitionRunner, not naive rayon. Both persistent matrix formats are column-oriented on disk: ColumnarCompactIntMatrix/ColumnarBitMatrix (obicompactvec/src/{intmatrix,bitmatrix}.rs) mmap one file per genome column; PackedCompactIntMatrix/PackedBitMatrix mmap one region-offset per column in a single file. fill_row(slot, buf) as used today (query.rs:262-272 via on_hit) reads one slot across all n_genomes columns per hit — the worst possible access pattern for this layout (up to n_genomes scattered mmap regions touched per single k-mer). Better: for each layer, walk the matrix column by column (genome by genome): for each genome, scan the slot keys collected in step 6 for that layer and call col.get(slot), keeping only nonzero results, and broadcast to the associated (seq_idx, pos) list. Total get() calls are unchanged (n_hits × n_genomes in the worst case) — the win is locality (sequential access within one mmap'd column at a time, not scattered across all columns per hit), not fewer operations. Columns are independent (read-only, disjoint mmap regions) → embarrassingly parallel across genomes, but — per point 7 — obicompactvec's existing into_par_iter() over 0..n_cols (sum(), count_nonzero(), pairwise distance matrices) is the naive, unpinned pattern the rest of the codebase is actively migrating away from, not a model to copy here. Route this through PartitionRunner (or the same NUMA-pool machinery) instead. Two things to settle when this is designed: how the partition axis (point 7), the column axis, and the existing chunk-level n_workers obipipeline pool compose without oversubscribing the machine (three different concurrency mechanisms — raw-thread pipe workers, PartitionRunner's pinned Rayon pools, and whatever drives the column scan — need a single reconciled thread budget, not three independent ones); and the threshold below which per-column dispatch overhead outweighs the gain (small n_genomes or small per-layer hit counts) — to be measured, not assumed.

9. Sparse per-genome representation, fed directly to Findere. Stage 2's output should be HashMap<genome_idx, Vec<(seq_idx, position, count)>>, sorted by (seq_idx, position) once collected, instead of a dense KmerResults-style matrix — the key must carry seq_idx, not just genome_idx, because a chunk batches many sequences and position is only meaningful within one; a plain Vec<(position, count)> per genome would silently mix positions from different sequences and corrupt the sliding-window scan. This bounds retained memory by actual nonzero hits on both axes (position sparsity from non-matching k-mers, genome sparsity from a matched k-mer typically belonging to only a handful of genomes out of possibly many). The Findere sliding-window (process_chunk, the win_min/deque loop) would need reworking to run per (sequence, genome) over its sparse, sorted (position, count) list — detect runs of ≥z consecutive positions, window-min within each run — instead of today's dense O(total_kmers × n_genomes) scan. This is also a genuine complexity win (O(hits log hits) per genome vs. dense scan), not just memory. Not covered by this sparsification: --detail's cov accumulator (query.rs:304-308) has the identical n_genomes-dense scaling problem and wasn't folded into points above. It doesn't need to be retained densely throughout processing, though — only the final JSON serialization (emit_batch) requires a dense [u32] per (seq, genome), and only for the sequences actually being output with --detail. Densification can stay a late, output-time-only step, reconstructed from the sparse per-genome lists.

Secondary patterns available from scatter.rs/superkmer.rs, not yet in query.rs:

• throttle() + CommonArgs::effective_max_open() to bound concurrently-open input files (query.rs defines its own QueryArgs, doesn't reuse this).
• Progress bar with EMA throughput + live active-worker gauges (obisys::spinner, flat_active/transform_active counters) — diagnostic value for locating the bottleneck.
• obisys::Reporter/Stage::start/stop timing per phase (used by index, filter; absent from query).

None of this is implemented yet — parked here as a coherent roadmap while the design is discussed further. Suggested dependency order: (1) I/O parallelism → (3) genome-aware chunk sizing → (4)–(9) staged/k-mer-deduped/NUMA-aware-partition-and-column-major/sparse query engine (larger refactor, biggest structural payoff — reuses PartitionRunner rather than inventing a new parallelism mechanism) → (2) parallel gzip (separate, orthogonal, tracked in chunkreader.md) → secondary diagnostics patterns.

---

Implementation plan

Concrete, phased translation of the roadmap above. Phases 0–2 are small, independent, low-risk, and each individually testable against current query output — land them first, in order, and measure on the reference 192-core/8-NUMA machine before deciding whether phases 3–5 (the staged/sparse engine, the larger structural payoff) are still worth their cost. Phases 3–5 are one coordinated change spanning obikmer, obikpartitionner, and obicompactvec — they should not be split across releases mid-way, because the intermediate state (e.g. k-mer-level dedup feeding the old dense KmerResults) has no correctness or performance benefit on its own. Phase 6 is unrelated to phases 0–5 and can happen any time, independently, if rapidgzip-rs is validated (see chunkreader.md).

Instrumentation is deliberately sequenced before the I/O fix (reordering the roadmap's own listed order), because every later phase's justification rests on a measurement ("to be measured, not assumed" appears throughout the roadmap above) — without it, phases 3–5 would be undertaken on faith.

Performance measurement on the reference 192-core/8-NUMA machine is done by the project owner, not from this development environment (macOS, 16 cores — PartitionRunner's NUMA pinning is Linux-only, so even phase 4's mechanism can't be functionally exercised for its actual purpose here). Each phase below is therefore written to be self-measuring: the debug-level logging it adds must be enough, on its own, to judge whether that phase's algorithmic choice paid off from a cluster run's logs, without needing to attach a profiler.

Conventions applied to every phase below

Debug logging. Every phase that changes an algorithmic choice (not phase 0, which is the logging) adds tracing::debug!/trace! at points that let a cluster run's logs answer "did this help": counts, ratios, and timings that quantify the specific claim that phase makes — e.g. phase 3 must log how many MPHF find calls were saved by k-mer-level dedup (the whole justification for that phase), phase 4 must log per-column scan timings, phase 5 must log actual retained-memory / sparsity ratios achieved. Prefer one structured debug! per chunk (fields, not prose) over free-text — the cluster logs will be the only evidence available for judging these choices, so they need to be grep/awk-able, not just readable.

Unit tests. This project's convention (obiread, obikseq, obidebruinj, obicompactvec, obilayeredmap, obiskio, obifastwrite) is #[cfg(test)] #[path = "tests/<name>.rs"] mod tests; at the bottom of the source file, with the actual test code in a sibling src/tests/<name>.rs. Neither obikmer nor obikpartitionner (the two crates phases 3 and 5 touch most) currently have a src/tests/ directory at all — this needs creating, following the existing pattern exactly, not inventing a new one.

Workflow (jj). Work happens in a fresh jj commit, easy to abandon. jj new between phases is reasonable where it helps isolate a phase for review, but only when the working copy compiles at that point (project convention) — phase 3's internal sub-steps (batch dedup change, then query_layer.rs split, then the new return shape) will likely not each compile independently since they're one coupled change, so treat "commit boundary" and "plan phase boundary" as related but not forced to match 1:1; use judgement per phase rather than mechanically splitting on every bullet.

Phase 0 — Instrumentation (prerequisite for measuring every later phase)

Goal: make core utilization, throughput, and per-stage timing visible on a real run, so phases 1–5 can be justified with numbers instead of assumption.

• obikmer/src/cmd/query.rs: wrap run()'s main loop with obisys::Reporter/Stage::start("query")/.stop(), printed at the end via rep.print() — same pattern as index.rs/filter.rs.
• Add an obisys::spinner("query") progress bar around the pipe.apply(...) loop, with an EMA throughput readout (bases/s or k-mers/s, mirroring steps::scatter's ema_rate computation, scatter.rs:88-118) and live gauges for "chunks in flight" / "workers busy" — reuse the AtomicU32 counter pattern from scatter.rs (flat_active, transform_active) rather than inventing a new one.
• Add max_open_files: Option<usize> to QueryArgs and a effective_max_open() method mirroring CommonArgs::effective_max_open() (obikmer/src/cli.rs:90-94) — needed by phase 1's throttle() call. (QueryArgs can't just embed CommonArgs — it doesn't take kmer_size/minimizer_size/partitions/level_max/theta from the CLI, those come from the index metadata — so this is a small standalone addition, not a flatten.)
• Add one structured debug! per process_chunk call: chunk byte size, sequence count, s-mer count, wall time, and (once later phases exist) the fields they add — this single log line is the baseline every later phase's own logging gets compared against.
• Validation: none needed beyond "the numbers appear and look sane" — this phase changes no query logic or output.
• Deliverable used by every phase below: a before/after throughput and core-utilization measurement on the reference machine.

Phase 1 — Parallel per-file I/O (fixes root cause of low core utilization)

Goal: file opening, decompression, and chunk-boundary parsing run across the n_workers pool instead of serially in the pipe's dedicated source thread.

• obikmer/src/cmd/query.rs:
  • Replace the paths.into_iter().flat_map(read_sequence_chunks_sized(...)) construction (current run(), building all_chunks) with obipipeline::throttle(paths.into_iter(), args.effective_max_open()), passed as the pipe's input.
  • Add a new QueryData::Path(PathBuf) variant (alongside Chunk/Output) to carry the throttled path through the pipe's type-erasure mechanism.
  • Add a new first pipe stage, Flat/fallible (||?), modeled on scatter.rs:60-86 and superkmer.rs:54-65: given a Throttled<PathBuf>, call read_sequence_chunks_sized(path, chunk_bytes) and yield each Rope chunk, keeping pw.guard alive until the file's iterator is exhausted (reuse or adapt scatter.rs's GuardedIter wrapper — same lifetime problem, same fix).
  • The existing process_chunk transform stage becomes the pipe's second stage, unchanged in its own logic — it still receives one Rope chunk at a time, just no longer all coming from one serial source.
  • make_pipe! invocation grows from one stage (Chunk => Output) to two (Path => Chunk => Output).
• Log, per file: time spent waiting on the throttle() slot (queueing due to max_open), and time spent opening/decompressing/producing the first chunk — this is what directly proves (or disproves) that I/O is now spread across workers instead of serialized.
• Validation: run query on a small multi-file input, diff output against the pre-change version — content must be identical; record order across files is not guaranteed to be preserved even before this change (chunk-level dispatch across n_workers already reorders completions), so the diff must be order-insensitive (sort by read id, or compare as sets) if it wasn't already.
• Measure: core utilization on the reference machine with several large input files, compare against phase 0's baseline.

Phase 2 — Genome-aware chunk-size formula (fixes OOM)

Goal: chunk_bytes reflects actual per-chunk memory (O(n_genomes)), not a fixed multiplier.

• obikmer/src/cmd/query.rs, run(): n_genomes and args.detail are already computed above the chunk_bytes calculation (n_genomes at the top of run(), before line 407 in the current file) — reorder if needed, then replace:

  let computed = avail / (n_workers as u64 * 16);
  

  with a formula that scales the divisor by n_genomes (and roughly doubles it when --detail is set, since cov duplicates the per-genome accumulation): e.g. per_chunk_multiplier = base_overhead + n_genomes as u64 * BYTES_PER_KMER_PER_GENOME * if detail { 2 } else { 1 }, replacing the flat 16. BYTES_PER_KMER_PER_GENOME should be derived from KmerResults's actual layout (4 bytes per u32 entry in data, plus the bool in in_index, plus win_min's equal-sized buffer) rather than guessed.
  • args.chunk_size (manual --chunk-size override) keeps taking priority, unchanged.
• Log the resolved chunk_bytes, n_genomes, and the estimated peak per-chunk memory (chunk_bytes × the same multiplier used to derive it) once at startup — lets a cluster run confirm the estimate was actually respected, not just that the process didn't get OOM-killed (which could also happen to be true for the wrong reason).
• Validation: build a test index with a large n_genomes (e.g. hundreds), run query with default chunk sizing under a memory limit (ulimit -v or a cgroup), confirm it no longer gets OOM-killed and that memory scales as predicted when n_genomes grows.
• Note: this phase is superseded once phase 5 lands (sparse retained memory no longer scales with n_genomes × total_kmers at all) — but it's needed immediately regardless, since phases 3–5 are a bigger, riskier change and users need a working query in the meantime.

Phase 3 — K-mer-level dereplication, staged MPHF/matrix lookup

Goal: replace superkmer-level dedup with k-mer-level dedup (roadmap point 5), and split the fused MPHF-find/matrix-fetch (point 4) so stage 1's output is bucketed by layer and MPHF slot (point 6).

• obikmer/src/cmd/query.rs:
  • Replace QueryBatch::from_records's dedup map (HashMap<RoutableSuperKmer, Vec<SKDesc>>, current query.rs:112) with a per-partition HashMap<CanonicalKmer, Vec<(seq_idx: u32, pos: u32)>>, built in the same SuperKmerIter pass: superkmer construction and partition routing (part_idx from the superkmer's minimizer hash) are unchanged, only the granularity of what gets deduplicated changes — each CanonicalKmer within a superkmer is inserted individually instead of the whole superkmer being the dedup key.
  • Verified: CanonicalKmer (obikseq/src/kmer.rs:390, pub type CanonicalKmer = CanonicalKmerOf<KLen>) — the underlying CanonicalKmerOf<L> derives Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash (kmer.rs:269). Usable as a HashMap/HashSet key as-is, no change needed.
• obikpartitionner/src/query_layer.rs:
  • Split QueryLayer::find_into (query_layer.rs:48-67) into two methods: find_slot(&self, kmer: CanonicalKmer) -> Option<usize> (MPHF only, no matrix touch) and keep fill_row as-is for phase 4 to call later.
  • Replace query_partition_with's inner loop (query_layer.rs:103-113) with a version that, for each unique CanonicalKmer, calls find_slot across the partition's layers (stopping at first hit, same as today), and instead of immediately filling a row, records (layer_idx, slot).
  • New return shape for the partition-level query, replacing today's on_hit(sk_idx, kmer_idx, row) callback: HashMap<layer_idx, HashMap<slot, Vec<(seq_idx, pos)>>> (roadmap point 6) — built directly from the k-mer dedup map's Vec<(seq_idx,pos)> values, keyed by the resolved slot instead of the k-mer.
• This phase alone has no throughput benefit yet (matrix fetch still happens, just deferred) beyond the k-mer-level dedup itself (fewer MPHF calls when queries have overlapping/repeated k-mers) — its purpose is to produce the input phase 4 needs. Land phase 3+4 together, not phase 3 alone, per the "don't split 3–5 across releases" note above.
• Log, per chunk: total k-mer occurrences vs. unique CanonicalKmer count (the dedup ratio — the entire justification for this phase) and the resulting MPHF find call count. If the dedup ratio is close to 1.0 on real query data (little redundancy), that's the cluster run telling us this phase wasn't worth it — the logging needs to be able to say that, not just confirm the happy path.
• Unit tests: create obikmer/src/cmd/tests/query.rs (new src/tests/ dir for this crate, following the project's #[cfg(test)] #[path = "tests/query.rs"] mod tests; convention) and obikpartitionner/src/tests/query_layer.rs (likewise new for this crate). Cover: the k-mer-level dedup map construction on synthetic sequences with known repeated/overlapping k-mers (assert unique-kmer count and occurrence lists); the find_slot/bucket-by-layer-and-slot construction against a small hand-built QueryLayer fixture, asserting the (layer_idx, slot, seq_idx, pos) tuples match what the old per-occurrence loop would have produced.

Phase 4 — Column-major matrix fetch, NUMA-aware (roadmap points 7–8)

Goal: replace fill_row-per-hit (row-major, worst-case mmap locality) with a column-major scan, driven through PartitionRunner rather than ad-hoc rayon, and extend the same mechanism to the partition loop itself.

• obicompactvec: add public per-column accessors alongside the existing pub(crate) ones — e.g. PersistentCompactIntMatrix::col_get(c: usize, slot: usize) -> u32 and equivalent for PersistentBitMatrix/Packed* variants (the underlying col(c)/col_view(c) + .get(slot) machinery already exists internally, per intmatrix.rs/bitmatrix.rs — this phase only needs to expose it).
• obikpartitionner: new function, e.g. fetch_layer_columns(layer: &QueryLayer, hits: &HashMap<slot, Vec<(seq_idx,pos)>>, n_genomes: usize) -> HashMap<genome_idx, Vec<(seq_idx, pos, count)>> — for each genome column g in 0..n_genomes, scan hits's slot keys, call the new col_get, keep nonzero results, append to the per-genome accumulator.
• obikmer/src/cmd/query.rs (or wherever this ends up being orchestrated): drive both the per-partition loop (today query.rs:250-278) and, within each partition/layer, the per-genome column loop, through obikindex::numa::PartitionRunner::run() (obikindex/src/numa.rs) instead of a plain loop or rayon::into_par_iter(). obikmer already depends on obikindex, so this needs no new dependency.
• Composition with existing parallelism — open design decision, resolve empirically before finishing this phase (flagged as unresolved in the roadmap, point 8): three concurrency mechanisms now coexist — obipipeline's raw-thread chunk-worker pool, PartitionRunner's NUMA-pinned Rayon pools, and whichever axis (partition or column) is chosen as PartitionRunner's unit of work. Proposed starting point to measure, not a final answer:
  1. Reduce obipipeline's n_workers for the chunk pipe to something small and NUMA-node-shaped (e.g. one chunk in flight per NUMA node) once PartitionRunner is doing the fine-grained fan-out within process_chunk — avoids the two mechanisms both trying to claim all 192 cores independently.
  2. Start with PartitionRunner driving the column axis only (point 8) inside each partition/layer, keep the partition loop (point 7) sequential within that; add partition-level PartitionRunner usage only if profiling still shows idle cores with many small partitions.
  3. Reuse PartitionRunner's existing dual CPU/IO activation signal as-is (numa_partition_runner.md's IoSample/CpuSample) rather than tuning new thresholds up front — the column scan's mmap-page-fault-bound shape is close enough to pack_matrices's documented profile to start from the same constants.
• Log, per layer processed: number of genome columns scanned, number of get() calls, wall time, and (via PartitionRunner's own on_done) active-worker count over time per NUMA node — this is the only way to see, from cluster logs alone, whether the composition question above actually got resolved well or whether one mechanism is starving the other.
• Unit tests: obikpartitionner/src/tests/query_layer.rs (same file as phase 3, extended) — cover fetch_layer_columns on a small hand-built multi-column matrix fixture (mix of Columnar/Packed/Implicit), asserting the sparse (seq_idx, pos, count) output matches a naive row-major reference implementation computed in the test itself. obicompactvec/src/tests/ (existing tests/mod.rs umbrella, per its established convention) gets the new col_get accessors' own unit tests alongside its current matrix tests.
• Validation: microbenchmark fetch_layer_columns against today's per-hit fill_row loop on a large-n_genomes layer, measuring wall time and (on Linux) page-fault counters (/proc/self/stat minor/major faults) — the claimed win is locality, so it should show up there even before end-to-end query numbers move.

Phase 5 — Sparse Findere rework (roadmap point 9)

Goal: replace the dense KmerResults/win_min sliding-window scan with one operating on phase 4's sparse per-genome output.

• obikmer/src/cmd/query.rs, process_chunk:
  • Remove KmerResults (query.rs:157-202) and the dense win_min allocation (query.rs:290-291, sized max_n_kmers × n_genomes).
  • Keep a lightweight dense in_index: Vec<bool> per chunk (sized total_kmers, independent of n_genomes) from phase 3's stage 1 — still needed for kmer_missing bookkeeping (leftmost-s-mer-of-window membership test), which phase 4's sparse structure doesn't carry (a k-mer with no genome hit has no entry there at all).
  • New per-(seq_idx, genome) scan: for each genome's Vec<(seq_idx, pos, count)> (sorted, per phase 4), group by seq_idx (contiguous after sort), then within each sequence's positions detect runs of pos, pos+1, pos+2, ... of length ≥ z; within each run, the existing monotone-deque window-minimum logic (query.rs's current dq loop, conceptually unchanged) applies — but the deque now only scans real entries in the run, never zero-filled gaps.
  • Update SeqAcc accumulation and emit_batch to consume this per-genome sparse iteration instead of results.val/results.is_in_index.
  • --detail/cov: build sparsely during the same scan (only positions with a confirmed contribution get an entry), densify into the [u32] JSON array only in emit_batch, only for genomes/sequences actually being serialized (per roadmap point 9's note, query.rs:304-308's current dense allocation goes away).
• Log, per chunk: total sparse entries retained vs. what the old dense KmerResults would have allocated (total_smers × n_genomes) — the sparsity ratio is this phase's entire reason for existing, so it must be directly visible in the logs, not inferred from process RSS. Also log the run-detection stats (number of runs found, average run length) — a low average run length relative to z would mean most positions still fail to form a full window, worth knowing.
• Unit tests: obikmer/src/cmd/tests/query.rs (extended from phase 3) — the property test described below is the primary deliverable here, not an afterthought; write it as an actual #[test] (or a small internal fuzz/property-style loop over randomized fixtures if a property-testing crate isn't already a dependency — check before adding one, per this project's dependency-approval rule) rather than a one-off manual comparison.
• Validation — this is the correctness-critical phase: property-test comparing old (dense, pre-phase-3) and new (sparse) implementations on the same randomized input/index fixtures, asserting identical kmer_count, kmer_missing, kmer_strict_matches, and (with --detail) coverage for every sequence. Keep both implementations compiled side by side (behind a debug-only flag or a temporary parallel code path) only for the duration of this validation; delete the dense path once parity is confirmed — per this project's own convention, superseded code is not kept "just in case."
• Update docmd/architecture/query.md itself: once this phase lands, the "Findere z-window filter" section (which currently — correctly — describes the dense deque-over-0..n_smers scan) needs another pass to describe the sparse run-detection algorithm instead, as already flagged when this phase was discussed.

Phase 6 — Parallel gzip decompression (independent, optional)

Tracked separately in chunkreader.md; parked pending validation of rapidgzip-rs on real data. Not a dependency of, or a dependency for, phases 0–5 — xopen is shared infrastructure (obiread), phase 1 benefits from it but doesn't require it (phase 1 parallelises across files; this phase would additionally parallelise within one large file).

Cross-cutting risks

• Thread-budget oversubscription (phase 4): the single biggest unresolved design question in this whole plan — see phase 4's composition note. Should be settled with real measurements early in phase 4, not assumed from the design alone.
• obicompactvec API surface growth (phase 4): new public per-column accessors are additive (existing fill_row/row stay for other callers — dump, select, distance computations) — no breaking change expected, but worth checking obicompactvec's other callers aren't already relying on fill_row being the only/cheapest access path in a way that would make maintaining two access patterns (row-major and column-major) a real maintenance cost rather than a one-off addition.
• PersistentBitMatrix::Implicit's hardcoded n_cols: 1 — resolved, not a bug. LayerMeta's own doc comment (obicompactvec/src/layer_meta.rs:1-9) states it is written "alongside mphf.bin" and read by PersistentBitMatrix::open "to determine n_rows for the implicit (mono-genome presence/absence) case" — i.e. Implicit is a documented single-genome fast path (no presence matrix needed when there is trivially one genome), not a generic "no matrix built yet" fallback. n_cols: 1 is correct by design for the case it's meant to handle. Phase 4's column loop is safe as planned — this was worth checking once, doesn't need further action.