@@ -2319,6 +2369,83 @@ scan the full unsampled index — never threaded --subsample/
out of scope here since the reported problem was specifically about the
--sankoff/--tnt pipeline's redundant/inconsistent scans, not these
two standalone flags.
+
query never benefits from sparse row-major access (found 2026-08-19, not implemented)
+
Benchmarked obikmer query against global_index_presence (dense-packed)
+vs. global_index_presence_sparse (pack --sparse), 100k simulated reads
+× 2 specimens (benchmark/, see
+benchmark_query_testing.md).
+Correctness: 0 mismatches — sparse and dense return bit-identical query
+results. Performance: sparse consistently slower than dense (~30-50%,
+reproducible across two runs with warm disk cache), the opposite of
+pack --sparse's stated intent ("faster for single-row access... like
+query").
+
Root cause, read from source, not measured in isolation:
+KmerPartition::query_partition_with (obikpartitionner/src/query_layer.rs:155-220)
+is architecturally column-major: stage 2 walks for g in 0..n_cols { for
+slot in hit_slots { layer.col_value(g, slot) } }, documented (correctly)
+as the right locality strategy for the packed/columnar formats, where
+col_value → PersistentBitMatrix::get is a genuine O(1) mmap'd column
+read (persistent.rs:110-113).
+
For Self::Sparse, that same get(c, slot) (persistent.rs:114-118)
+allocates a full n_cols-wide buffer and calls fill_row — materializing
+the entire row — just to return one cell. Called from inside the
+column-major double loop, this reconstructs the same row once per genome
+column touched: O(hits × n_cols) full-row rebuilds instead of O(hits).
+PersistentSparseBitMatrix's own native row-major decode
+(for_each_genome_in_row, sparse.rs:164-177, used correctly by its own
+row/fill_row/fill_sub_matrix) is never reached from the query path
+at all.
+
fill_sub_matrix (the existing BinaryMatrix trait primitive,
+traits.rs:13-37) is not the right replacement for query either, even
+once its own dispatch bug is fixed (see next section) — its output shape
+is inherently column-dense: out[col] gets an entry for every column,
+including columns with zero hits among the requested slots. On real
+sparse data (a hit typically touching a handful of genomes out of dozens)
+that's still O(n_cols) output regardless of true sparsity. What query
+actually wants is the sparse triple stream (slot, col, value) it already
+consumes as QueryHit::Value — not a materialized sub-matrix.
+
Proposed primitive (design only, not implemented — explicit ask: keep
+count matrices not excluded, even though effort right now is
+presence/absence only):
+
/// Visit every nonzero cell among `slots`. Order unspecified.
+fnfor_each_nonzero(&self,slots:&[usize],f:implFnMut(usize/*idx into slots*/,usize/*col*/,u32/*value*/));
+
+
+
On PersistentSparseBitMatrix: native override, one pass per slot via
+ the existing (currently private) for_each_genome_in_row — O(Σ row
+ nnz), zero n_cols-wide allocation. This is the whole point: expose code
+ that already exists rather than write anything new for the sparse side.
+
On PersistentBitMatrix::{Packed,Columnar,Implicit}: provided
+ default, derived from fill_sub_matrix (materialize, then filter to
+ true cells) — reuses the already-optimal column-major/mmap path for
+ those formats, no new code needed there either.
+
On PersistentCompactIntMatrix (counts): same provided-default
+ treatment, derived from its own existing fill_sub_matrix (u32-typed
+ already, intmatrix.rs:387) — not hand-optimized (no sparse count
+ format exists — "Explicitly deferred" per traits.rs:9-12), but not
+ excluded either: it gets a working, not-pathological implementation for
+ free today, on the same trait, ready for a native override the day a
+ sparse count format lands. This is why the signature carries u32
+ rather than bool — presence is 0/1, counts are u32, one trait
+ covers both without a bool/u32 split forcing counts out of the design.
+
+
Would let query_partition_with's stage 2 collapse to one
+layer.matrix().for_each_nonzero(&hit_slots, |i, g, v| on_event(...))
+call per layer, format-agnostic, with each backend's existing (or
+default-derived) implementation deciding the actual access pattern.
+
Separately, an existing bug in the generic path (found while tracing
+this, itself not yet fixed): PersistentBitMatrix::fill_sub_matrix
+(persistent.rs:190-215, the enum wrapper backing BinaryMatrix's
+default trait impl) does not delegate to
+PersistentSparseBitMatrix::fill_sub_matrix for Self::Sparse — it
+reimplements the same naive per-(column, slot) fill_row_bool loop
+instead, bypassing the efficient native method one file over
+(sparse.rs:249-258). obikphylo::siblings::cache::Mat
+(cache.rs:138-145) independently built its own parallel enum wrapper
+that dispatches correctly — a sign this was worked around rather than
+fixed at the source. Any future for_each_nonzero work should fix this
+dispatch too (or route through it), rather than adding a third
+independently-dispatching wrapper.
benchmark/Makefile exercises indexing, merge, and phylo distance
+reconstruction against simulated bacterial genomes. It now also covers
+obikmer query — the read-matching path — and the sparse packed
+presence-matrix format (obikmer pack --sparse), previously untested by
+this pipeline.
+
Motivation
+
+
query had no end-to-end coverage. A regression there would not be caught
+ by verify_presence/verify_merge_presence, which only check index
+ content against the .npz truth, never the query API.
+
pack --sparse produces a presence-matrix format documented (see
+ siblings.md) as faster for single-row
+ access (query) and slower for column-oriented access (phylo --metric).
+ global_index_presence/ built by merge_presence.sh is always packed
+ dense (packing is a stage inside merge, not a separate pack
+ invocation) — there was no dense/sparse regression check.
+
+
Query read source
+
Query reads are independent of simulated_data/ (which is folded into the
+index being queried): reusing those reads would test against the exact
+error draw the index was built from. query_data/<species>/<strain>/ holds
+a second, independent iss generate run against the same reference
+genome, via simulate_query_one.sh — unseeded, so a second draw picks up
+different sequencing errors than simulate_one.sh's draw for the same
+genome. Fixed at 100,000 read pairs per genome (not coverage-proportional
+like the 15x used for simulated_data/), so wall/RSS numbers stay
+comparable across genomes of very different sizes.
+
Two query-source specimens, hardcoded as QUERY_SPECIMENS in
+make_deps.py: Escherichia_coli--K-12_MG1655 (common, well-represented
+bacterium) and Saccharolobus_islandicus--M.16.4 (the only archaeon in
+SPECIES — distant lineage, stresses the query path differently from a
+close-relative match). Two is enough to catch a dense/sparse regression
+without duplicating the exhaustive per-specimen coverage
+verify_merge_presence already provides across all SPECIMENS.
+
Sparse global index
+
global_index_presence_sparse/ is built by pack_sparse.sh: copy
+global_index_presence/ wholesale, then obikmer pack --sparse in place.
+This works directly because merge's pack stage (merge.rs:252,
+pack_matrices(false)) keeps the per-genome column files on disk after
+dense-packing — pack_sparse_bit_matrix (obicompactvec/src/bitmatrix/sparse.rs:447)
+reads those, is idempotent, and removes matrix.pbmx once the sparse form
+is written, so Persistent::open falls through to the sparse format
+afterward. No separate merge run needed.
+
Query runs
+
query_one.sh dense|sparse SPECIMEN runs obikmer query --count-missing
+against global_index_presence or global_index_presence_sparse, output
+gzipped to query_{dense,sparse}/SPECIMEN.fasta.gz, Reporter wall/RSS
+captured to stats/query_{dense,sparse}/SPECIMEN.stats (same
+stderr-parsing convention as merge_presence.sh).
+
Flags: --count-missing only. --mismatch is a no-op today
+(query/mod.rs:212-213, prints "not yet implemented, ignored") — left off
+rather than tested for a feature that doesn't exist yet.
+
Dense/sparse regression
+
verify_query.py compares the two query outputs per specimen, matched by
+read id (not stream position — the query pipeline chunks input across
+worker threads and doesn't guarantee output order). Compares kmer_count,
+kmer_missing, and the full kmer_strict_matches map per read. Any
+mismatch is a real regression: dense and sparse must be content-identical,
+only I/O access pattern differs. .stats → stats/verify_query/,
+aggregated by aggregate_stats.sh query|verify_query-style cases
+(query_dense, query_sparse, verify_query).
+
Performance comparison
+
No dedicated script: the wall/RSS columns from the query_dense and
+query_sparse aggregated .stats CSVs are the dense-vs-sparse performance
+comparison — the expected win for query on sparse, per the pack --sparse
+help text.
+
Scope
+
count track excluded from the sparse branch: pack --sparse targets
+presence matrices only (per CLI help); pack_matrices leaves count
+matrices untouched regardless of the sparse flag
+(obikindex/src/index.rs:308).
+
New Makefile targets
+
simulate_query, pack_sparse, query_dense, query_sparse,
+aggregate_query_dense, aggregate_query_sparse, verify_query,
+aggregate_verify_query — the last three folded into all.
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
\ No newline at end of file
diff --git a/DevDoc/implementation/chunkreader.refs/index.html b/DevDoc/implementation/chunkreader.refs/index.html
index 70a9ae3a..3be4d1b2 100644
--- a/DevDoc/implementation/chunkreader.refs/index.html
+++ b/DevDoc/implementation/chunkreader.refs/index.html
@@ -944,6 +944,34 @@
+
+
+
+
+
+
+
Both matrix types are enums behind a transparent API — the caller never matches on the variant. PersistentCompactIntMatrix has two variants (Columnar, Packed). PersistentBitMatrix has four:
-
-
Columnar format
-
Packed format
+
Variant
+
Storage
+
When
-
Bit
-
PersistentBitMatrix (Columnar variant)
-
PersistentBitMatrix (Packed variant)
+
Columnar
+
one .pbiv/.pciv file per column + meta.json
+
build-time default (*Builder::new)
-
Int
-
PersistentCompactIntMatrix (Columnar variant)
-
PersistentCompactIntMatrix (Packed variant)
+
Packed
+
single matrix.pbmx mmap file
+
query-optimised, produced by pack_bit_matrix/pack_compact_int_matrix
pack --sparse; see siblings.md for the sparse-vs-dense access-pattern trade-off
+
+
+
Implicit (bit only)
+
no file at all
+
mono-genome presence layers — n_cols is always reported as 1, every value is true
-
Both matrix types are enums (Columnar / Packed / Implicit for bit) behind a transparent API. col_view(c) returns the appropriate view directly:
+
PersistentBitMatrix::open(layer_dir) auto-detects the variant, in order: matrix.pbmx → Packed, presence/meta.json → Columnar, presence/sparse_meta.json → Sparse, layer_meta.json (no presence dir at all) → Implicit. col_view/col/sub_matrix panic on Sparse/Implicit where the operation has no direct-slice equivalent (Sparse is k-mer-major, not column-major; Implicit has no backing storage) — callers needing per-column data on those variants go through row/fill_row.
+
col_view(c) returns the appropriate view directly:
Mode 3 (PersistentBitMatrix) has no push_layer on LayeredMap; callers build directly via Layer<PersistentBitMatrix>::build_presence.
+
Layer\<D> — raw mapping, iteration, and batch access
+
Beyond query/find (membership-checked), Layer<D> exposes lower-level access used by consumers that already know a kmer is in the layer (e.g. cross-partition sibling resolution) or that need to sweep every kmer/slot without paying for a membership check each time.
Pure MPHF mapping, no evidence/fingerprint check — equivalent to MphfOnly::index. Only meaningful when the caller already knows kmer belongs to the layer; on an absent kmer the MPHF still returns some slot (undefined, not None).
+
Kmer iteration
+
Four iterators, all built from unitigs.bin (physical layout order, not correlated with MPHF slot numbers):
+
pubfniter_kmers(&self)->KmerIter<'_>
+pubfnenumerate_kmers(&self)->Enumerate<KmerIter<'_>>// (order_index, kmer)
+pubfniter_kmers_batch(&self,n:usize)->KmerBatchIter<'_>// Vec<CanonicalKmer> of size ≤ n
+pubfnenumerate_kmers_batch(&self,n:usize)->implIterator<Item=(usize,Vec<CanonicalKmer>)>+Send+'static
+
+
KmerIter/KmerBatchIter own a clone of the underlying Arc<UnitigFileReader> rather than borrowing self — Send + 'static, streamed from disk one kmer at a time, never materialised as a whole. Multiple instances can coexist concurrently, each with its own cursor. enumerate_kmers_batch's index is the batch's starting offset in iteration order (a multiple of n except for the final, possibly shorter, batch).
+
Batch lookup on payload vectors/views
+
PersistentCompactIntVec, PersistentBitVec, IntSliceView, BitSliceView all expose:
Both sort slots internally for sequential mmap access, then reorder results back to the caller's original order. fill_batch fills a caller-provided buffer, avoiding the Vec allocation.
+
sub_matrix / fill_sub_matrix
+
// Layer<PersistentCompactIntMatrix>
+pubfnsub_matrix(&self,slots:&[usize])->Vec<Vec<u32>>// column-first
+pubfnfill_sub_matrix(&self,slots:&[usize],out:&mut[Vec<u32>])
+
+// Layer<PersistentBitMatrix> (and any D: BinaryMatrix, e.g. PersistentSparseBitMatrix)
+pubfnsub_matrix(&self,slots:&[usize])->Vec<Vec<bool>>
+pubfnfill_sub_matrix(&self,slots:&[usize],out:&mut[Vec<bool>])
+
+
Column-first to match the on-disk column-major layout. fill_sub_matrix sorts slots once, then calls each column's fill_batch in turn — no redundant per-column sort. On PersistentSparseBitMatrix (k-mer-major, no column method) this degrades to a row-by-row decode; see siblings.md.
+
LayeredStore\<S> and aggregation traits
LayeredStore<S> is a generic aggregation wrapper over Vec<S>. It propagates three traits from obicompactvec::traits up the hierarchy via blanket impls:
+
+
+
+
diff --git a/DevDocMD/architecture/siblings.md b/DevDocMD/architecture/siblings.md
index 1247cbdb..56396eb7 100644
--- a/DevDocMD/architecture/siblings.md
+++ b/DevDocMD/architecture/siblings.md
@@ -756,3 +756,129 @@ scan the full unsampled index — never threaded `--subsample`/`--entropy`,
out of scope here since the reported problem was specifically about the
`--sankoff`/`--tnt` pipeline's redundant/inconsistent scans, not these
two standalone flags.
+
+## `query` never benefits from sparse row-major access (found 2026-08-19, not implemented)
+
+Benchmarked `obikmer query` against `global_index_presence` (dense-packed)
+vs. `global_index_presence_sparse` (`pack --sparse`), 100k simulated reads
+× 2 specimens (`benchmark/`, see
+[benchmark_query_testing.md](../implementation/benchmark_query_testing.md)).
+Correctness: 0 mismatches — sparse and dense return bit-identical query
+results. Performance: sparse consistently *slower* than dense (~30-50%,
+reproducible across two runs with warm disk cache), the opposite of
+`pack --sparse`'s stated intent ("faster for single-row access... like
+query").
+
+**Root cause, read from source, not measured in isolation:**
+`KmerPartition::query_partition_with` (`obikpartitionner/src/query_layer.rs:155-220`)
+is architecturally column-major: stage 2 walks `for g in 0..n_cols { for
+slot in hit_slots { layer.col_value(g, slot) } }`, documented (correctly)
+as the right locality strategy for the packed/columnar formats, where
+`col_value` → `PersistentBitMatrix::get` is a genuine O(1) mmap'd column
+read (`persistent.rs:110-113`).
+
+For `Self::Sparse`, that same `get(c, slot)` (`persistent.rs:114-118`)
+allocates a full `n_cols`-wide buffer and calls `fill_row` — materializing
+the *entire row* — just to return one cell. Called from inside the
+column-major double loop, this reconstructs the same row once per genome
+column touched: O(hits × n_cols) full-row rebuilds instead of O(hits).
+`PersistentSparseBitMatrix`'s own native row-major decode
+(`for_each_genome_in_row`, `sparse.rs:164-177`, used correctly by its own
+`row`/`fill_row`/`fill_sub_matrix`) is never reached from the query path
+at all.
+
+**`fill_sub_matrix` (the existing `BinaryMatrix` trait primitive,
+`traits.rs:13-37`) is not the right replacement for `query` either**, even
+once its own dispatch bug is fixed (see next section) — its output shape
+is inherently column-dense: `out[col]` gets an entry for every column,
+including columns with zero hits among the requested slots. On real
+sparse data (a hit typically touching a handful of genomes out of dozens)
+that's still O(n_cols) output regardless of true sparsity. What `query`
+actually wants is the sparse triple stream `(slot, col, value)` it already
+consumes as `QueryHit::Value` — not a materialized sub-matrix.
+
+**Proposed primitive** (design only, not implemented — explicit ask: keep
+count matrices *not excluded*, even though effort right now is
+presence/absence only):
+
+Not a closure-driven `for_each` — a real `Iterator`, one concrete struct
+per matrix format, so the traversal state (current position in the sorted
+slot list, current column, permutation, sparse-row decode cursor…) lives
+in named struct fields instead of being threaded implicitly through
+recursion or a captured closure. RPITIT (stable since 1.75, and this
+workspace is edition 2024) means the trait method can return it without
+naming or boxing the concrete type:
+
+```rust
+/// Yields every nonzero cell among `slots`, in implementation-defined order.
+fn nonzero_iter<'a>(&'a self, slots: &'a [usize]) -> impl Iterator + 'a;
+// item: (idx into `slots`, col, value)
+```
+
+This is the one foundational primitive per format — both `fill_sub_matrix`
+and a `for_each`-style callback become trivial, free consumers of it
+(`.for_each(f)` is already `Iterator::for_each` from std; `fill_sub_matrix`
+becomes "drain the iterator, scatter into `out[][]`"), instead of two
+independently-maintained traversals that can silently diverge (see the bug
+below — this is exactly how it happened).
+
+- **On `PersistentSparseBitMatrix`**: the struct is nearly free to write —
+ it wraps the existing (currently private) `for_each_genome_in_row`
+ per-row decode, advancing to the next `slots` entry on exhaustion. O(Σ
+ row nnz), zero `n_cols`-wide allocation.
+- **On `PersistentBitMatrix::{Packed,Columnar}`**: revised — cheaper than
+ first thought, by reusing the same split already used for
+ `fill_matrix`'s own implementation instead of hand-writing a resumable
+ state machine at the matrix level. The base-vector layer
+ (`BitSliceView`, `views.rs`) already separates the two concerns:
+ `fill_batch_sorted` (`views.rs:55-60`, sorted-slot batch lookup) and a
+ genuine per-bit `Iterator` (`BitSliceIter`, `views.rs:94+`) sit side by
+ side there, one level below the matrix. Adding a
+ "positions among `sorted_slots` where the bit is set" iterator at that
+ same vector level is a `filter` over the existing `get()` — no new
+ state machine, since `std::iter::Filter` already *is* one, generated by
+ the compiler:
+ ```rust
+ fn nonzero_among_sorted<'s>(&'s self, sorted_slots: &'s [usize]) -> impl Iterator + 's {
+ sorted_slots.iter().copied().filter(move |&slot| self.get(slot))
+ }
+ ```
+ The matrix-level `nonzero_iter` then composes these per column with
+ `flat_map` over `0..n_cols` (each column's hits, tagged with `c`,
+ slot mapped back through the sort permutation `fill_batch`/
+ `fill_batch_sorted` already carry) — again a combinator chain, not a
+ hand-rolled struct. Same algorithm, same mmap/sort locality as today's
+ `fill_sub_matrix`; just assembled from `std` iterator adaptors instead
+ of a loop body writing into a buffer, mirroring the vector/matrix split
+ the codebase already uses for `fill_batch_sorted` rather than
+ introducing a new shape.
+- **`Implicit`**: trivial (`slots.iter().map(|&i| (i, 0, 1))`, one column,
+ always present).
+- **On `PersistentCompactIntMatrix` (counts)**: same treatment as
+ `Packed`/`Columnar` — no sparse count format exists yet ("Explicitly
+ deferred" per `traits.rs:9-12`), so no native low-effort case the way
+ `Sparse` has one, but not excluded either: the iterator's `Item` is
+ already `(usize, usize, u32)`, not `bool`, specifically so presence
+ (`0`/`1`) and counts (arbitrary `u32`) share one primitive instead of a
+ bool/u32 split forcing counts out of the design. Ready for a native
+ sparse-count struct later without a signature change.
+
+Would let `query_partition_with`'s stage 2 collapse to one
+`for (i, g, v) in layer.matrix().nonzero_iter(&hit_slots) { on_event(...) }`
+per layer, format-agnostic, each backend's struct deciding the actual
+traversal.
+
+**This also closes the existing dispatch bug for free, by construction**:
+`PersistentBitMatrix::fill_sub_matrix` (`persistent.rs:190-215`, the enum
+wrapper backing `BinaryMatrix`'s trait impl) today does *not* delegate to
+`PersistentSparseBitMatrix::fill_sub_matrix` for `Self::Sparse` — it
+reimplements the same naive per-(column, slot) `fill_row_bool` loop
+instead, bypassing the efficient native method one file over
+(`sparse.rs:249-258`). `obikphylo::siblings::cache::Mat`
+(`cache.rs:138-145`) independently built its own parallel enum wrapper
+that dispatches correctly — evidence this was worked around rather than
+fixed at the source: two hand-written traversals for the same format,
+free to drift apart, and they did. If `fill_sub_matrix` itself is
+rewritten as "drain `nonzero_iter`, scatter into `out[][]`", there is only
+one traversal per format left to get right — the bug class doesn't just
+get fixed once, it stops being possible to reintroduce.
diff --git a/doc/sitemap.xml.gz b/doc/sitemap.xml.gz
index d240d53a38bf83f4e1976065c5477c6a68e9751f..425040cd954e2b1b059a7ef7e70cb82e7f9be68e 100644
GIT binary patch
delta 13
Ucmb=gXP58h;An7eo5)@P02}fI82|tP
delta 13
Ucmb=gXP58h;9yu&Kasrx02<>1`v3p{