From 3ba26b3dc1d9593c8dd78811c7dae585687b48a9 Mon Sep 17 00:00:00 2001 From: Eric Coissac Date: Sat, 15 Aug 2026 19:19:48 +0200 Subject: [PATCH] feat: add sparse on-disk format for presence matrices The index packing API now accepts a `sparse` parameter to generate `PersistentSparseBitMatrix` files alongside existing dense matrices. The sibling cache automatically detects this format via an `is_multi.prsb` marker file and routes queries identically to the dense variant. A new `--sparse` CLI flag exposes the option, with tests verifying end-to-end pipeline correctness and storage equivalence. --- .gitignore | 1 + docmd/architecture/siblings.md | 56 ++++++++++++++++++++++ src/obicompactvec/src/bitmatrix/mod.rs | 2 +- src/obicompactvec/src/bitmatrix/sparse.rs | 34 +++++++++++-- src/obicompactvec/src/lib.rs | 2 +- src/obicompactvec/src/traits.rs | 8 ++++ src/obikindex/src/index.rs | 20 ++++++-- src/obikindex/src/merge.rs | 2 +- src/obikindex/src/select.rs | 2 +- src/obikmer/src/cmd/pack/mod.rs | 10 +++- src/obikphylo/src/siblings/cache.rs | 22 ++++++++- src/obikphylo/src/siblings/tests.rs | 41 ++++++++++++++++ src/obilayeredmap/src/layer.rs | 45 ++++++++++++++---- src/obilayeredmap/src/tests/layer.rs | 58 ++++++++++++++++++++++- 14 files changed, 279 insertions(+), 24 deletions(-) diff --git a/.gitignore b/.gitignore index bfe28adc..5bb83c6d 100644 --- a/.gitignore +++ b/.gitignore @@ -16,6 +16,7 @@ benchmark/simulated_data benchmark/specimen_index_presence benchmark/specimen_index_count benchmark/global_index_presence +benchmark/global_index_presence_sav benchmark/all_specific benchmark/global_index_count benchmark/stats diff --git a/docmd/architecture/siblings.md b/docmd/architecture/siblings.md index eb032890..2a7d26c6 100644 --- a/docmd/architecture/siblings.md +++ b/docmd/architecture/siblings.md @@ -496,3 +496,59 @@ Also still deferred, unchanged from the implementation plan: full Alanko et al. subset-hierarchy compression (only the exact-duplicate special case is built), a sparse `PersistentCompactIntMatrix` (count matrices), and BRWT-style column-correlation exploitation. + +## Wired into `pack` and the sibling-annex build path (2026-08-15) + +`PersistentSparseBitMatrix` went from a validated but unused type to a +real, selectable on-disk format: + +- **Generic `Layer`**: `obilayeredmap::Layer`'s presence-only methods + (`n_cols`, `sub_matrix`, `fill_sub_matrix`) are generic over any + `D: LayerData> + BinaryMatrix`, not hardcoded to + `PersistentBitMatrix` — `PersistentSparseBitMatrix` implements + `LayerData` (`open`/`read`) the same way. `find_slot`/`index_batch` were + already generic over any `D: LayerData`, so they needed no change. + Verified by `obilayeredmap`'s + `presence_layer_generic_over_sparse_matches_dense` test: build a dense + presence layer, convert it to sparse via `build_from_dense`, open both + as `Layer`/`Layer` on + the same directory, assert `n_cols`/`sub_matrix`/`find_slot` agree. + (This test must stay at `k=4` with mutually non-colliding canonical + 4-mers across its input sequences — `K`/`M` are process-wide + `AtomicUsize`s in test builds, not thread-local, so a test using a + different `k` races every other test in the same crate binary; a k=11 + version of this test passed alone but failed under the full + `obilayeredmap` suite for exactly that reason before being fixed.) +- **`obikphylo::siblings::cache::Mat`** gained a third variant, + `SparsePresence(Layer)`, alongside `Count` + and `Presence` — every method (`find_slot`, `index_batch`, + `iter_minorants_batch`, `n_cols`, `fill_sub_matrix_carries`) dispatches + to it identically to `Presence`, since both go through the same generic + `Layer` code. `PartitionCache::build` picks the variant per layer by + checking for `presence/is_multi.prsb` (the sparse format's own marker + file, see the design section above) before falling back to the dense + open path. +- **`pack_sparse_bit_matrix`** (new, `obicompactvec::bitmatrix::sparse`): + `pack --sparse`'s entry point. Idempotent (checks `is_multi.prsb` + first); packs to dense `matrix.pbmx` first if that hasn't happened yet + (the dense→sparse transpose needs random row access, which only the + packed/columnar dense forms give), then `build_from_dense`s the sparse + form into the same directory and deletes `matrix.pbmx` — old-format + files are removed only after the new format is fully written, mirroring + `pack_bit_matrix`'s own crash-safety convention. +- **CLI**: `obikmer pack --sparse` threads a `sparse: bool` through + `KmerIndex::pack_matrices` (all other call sites — `select`, `merge`, + `finalize_indexed` — pass `false`, unchanged dense behaviour). Count + matrices are untouched by `--sparse` (no sparse `PersistentCompactIntMatrix` + — see "still deferred" above). +- **End-to-end coverage**: `obikphylo::siblings::tests:: + sibling_annex_works_after_pack_sparse` builds a two-genome index, packs + it `--sparse`, asserts `is_multi.prsb` exists, then runs + `build_sibling_annex` and checks the resulting `FamilyMask`s match the + dense-path test (`sibling_annex_one_sibling_each`) exactly — proves the + sparse format round-trips through the real build pipeline + (`PartitionCache` sparse-detection included), not just the + `obicompactvec`/`obilayeredmap` unit layers below it. + +Full workspace `cargo test` (all crates, unit + doc tests) green after +this change. diff --git a/src/obicompactvec/src/bitmatrix/mod.rs b/src/obicompactvec/src/bitmatrix/mod.rs index 0523da87..6435ca2b 100644 --- a/src/obicompactvec/src/bitmatrix/mod.rs +++ b/src/obicompactvec/src/bitmatrix/mod.rs @@ -22,7 +22,7 @@ mod sparse; pub use builder::PersistentBitMatrixBuilder; pub use packed::pack_bit_matrix; pub use persistent::PersistentBitMatrix; -pub use sparse::{PersistentSparseBitMatrix, PersistentSparseBitMatrixBuilder}; +pub use sparse::{PersistentSparseBitMatrix, PersistentSparseBitMatrixBuilder, pack_sparse_bit_matrix}; pub(crate) use pairwise::{pairwise_matrix, pairwise2_matrix}; diff --git a/src/obicompactvec/src/bitmatrix/sparse.rs b/src/obicompactvec/src/bitmatrix/sparse.rs index 09dea336..873a0554 100644 --- a/src/obicompactvec/src/bitmatrix/sparse.rs +++ b/src/obicompactvec/src/bitmatrix/sparse.rs @@ -5,7 +5,7 @@ //! used by any production code path yet — a new type, not a replacement. //! //! On-disk layout (a directory, mirroring [`super::PersistentBitMatrix`]'s -//! own directory-of-files convention): `meta.json` plus four components — +//! own directory-of-files convention): `sparse_meta.json` plus four components — //! `is_multi.prsb` (rank-capable flag: singleton row vs. multi-genome //! row), `singleton.pfiv` (genome index, one entry per singleton row, //! `ceil(log2(n_cols))` bits each), `multi.pfiv` (`dict_id`, one entry per @@ -33,13 +33,14 @@ use crate::meta::field; use crate::rankselect::{PersistentRankSelectBitVec, PersistentRankSelectBitVecBuilder}; use super::PersistentBitMatrix; +use super::packed::PackedBitMatrix; fn is_multi_path(dir: &Path) -> PathBuf { dir.join("is_multi.prsb") } fn singleton_path(dir: &Path) -> PathBuf { dir.join("singleton.pfiv") } fn multi_path(dir: &Path) -> PathBuf { dir.join("multi.pfiv") } fn dict_offsets_base(dir: &Path) -> PathBuf { dir.join("dict_offsets") } fn dict_values_path(dir: &Path) -> PathBuf { dir.join("dict_values.bin") } -fn meta_path(dir: &Path) -> PathBuf { dir.join("meta.json") } +fn meta_path(dir: &Path) -> PathBuf { dir.join("sparse_meta.json") } struct SparseMeta { n: usize, @@ -53,7 +54,7 @@ impl SparseMeta { fn load(dir: &Path) -> io::Result { let s = fs::read_to_string(meta_path(dir))?; let get = |name: &str| { - field(&s, name).ok_or_else(|| io::Error::new(io::ErrorKind::InvalidData, format!("bad meta.json: missing {name}"))) + field(&s, name).ok_or_else(|| io::Error::new(io::ErrorKind::InvalidData, format!("bad sparse_meta.json: missing {name}"))) }; Ok(Self { n: get("n")?, @@ -370,3 +371,30 @@ impl PersistentSparseBitMatrixBuilder { Ok(builder) } } + +/// `pack --sparse`'s entry point: converts a presence directory into the +/// sparse on-disk format in place, mirroring [`super::pack_bit_matrix`]'s +/// convention (old-format files removed only after the new format is +/// fully written, so a crash mid-conversion leaves the previous, still +/// valid format rather than a half-written one). Idempotent — does +/// nothing if `is_multi.prsb` already exists. Packs to dense +/// (`matrix.pbmx`) first via [`super::pack_bit_matrix`] if that hasn't +/// happened yet, since the dense→sparse transpose (`build_from_dense`) +/// needs random row access, which only the packed/columnar dense forms +/// give — not a further reason to keep the dense file around afterward. +pub fn pack_sparse_bit_matrix(dir: &Path) -> io::Result<()> { + if is_multi_path(dir).exists() { + return Ok(()); + } + + let packed_path = dir.join("matrix.pbmx"); + if !packed_path.exists() { + super::pack_bit_matrix(dir)?; + } + + let dense = PersistentBitMatrix::Packed(PackedBitMatrix::open(&packed_path)?); + PersistentSparseBitMatrixBuilder::build_from_dense(&dense, dir)?.close()?; + drop(dense); + fs::remove_file(&packed_path)?; + Ok(()) +} diff --git a/src/obicompactvec/src/lib.rs b/src/obicompactvec/src/lib.rs index 58450c47..a7ce5df3 100644 --- a/src/obicompactvec/src/lib.rs +++ b/src/obicompactvec/src/lib.rs @@ -19,7 +19,7 @@ pub use bitvec::{BitIter, PersistentBitVec, PersistentBitVecBuilder}; pub use fixedintvec::{PersistentFixedIntVec, PersistentFixedIntVecBuilder, bit_width_for_range}; pub use rankselect::{PersistentRankSelectBitVec, PersistentRankSelectBitVecBuilder}; pub use eliasfano::{EliasFano, EliasFanoBuilder}; -pub use bitmatrix::{PersistentBitMatrix, PersistentBitMatrixBuilder, PersistentSparseBitMatrix, PersistentSparseBitMatrixBuilder, pack_bit_matrix}; +pub use bitmatrix::{PersistentBitMatrix, PersistentBitMatrixBuilder, PersistentSparseBitMatrix, PersistentSparseBitMatrixBuilder, pack_bit_matrix, pack_sparse_bit_matrix}; pub use builder::PersistentCompactIntVecBuilder; pub use colgroup::{ColGroup, FilterMask, MatrixGroupOps, eval_filter_mask}; pub use intmatrix::{PersistentCompactIntMatrix, PersistentCompactIntMatrixBuilder, pack_compact_int_matrix}; diff --git a/src/obicompactvec/src/traits.rs b/src/obicompactvec/src/traits.rs index 5437b7a1..62692031 100644 --- a/src/obicompactvec/src/traits.rs +++ b/src/obicompactvec/src/traits.rs @@ -26,6 +26,14 @@ pub trait BinaryMatrix { fn fill_sub_matrix(&self, slots: &[usize], out: &mut [Vec]); /// Per-genome k-mer totals. fn count_ones(&self) -> Array1; + + /// Allocating variant of [`fill_sub_matrix`](Self::fill_sub_matrix) — + /// provided once here so implementers only need the fill-in-place form. + fn sub_matrix(&self, slots: &[usize]) -> Vec> { + let mut out: Vec> = (0..self.n_cols()).map(|_| Vec::new()).collect(); + self.fill_sub_matrix(slots, &mut out); + out + } } /// Convert a Jaccard distance matrix (`1 - J`) into a Mash distance matrix, per diff --git a/src/obikindex/src/index.rs b/src/obikindex/src/index.rs index 6b7e117d..c8291060 100644 --- a/src/obikindex/src/index.rs +++ b/src/obikindex/src/index.rs @@ -122,7 +122,7 @@ impl KmerIndex { fs::File::create(output.join(SENTINEL_INDEXED)).map_err(OKIError::Io)?; let idx = KmerIndex::open(output)?; let t_pack = Stage::start("pack"); - idx.pack_matrices()?; + idx.pack_matrices(false)?; rep.push(t_pack.stop()); Ok(idx) } @@ -274,8 +274,14 @@ impl KmerIndex { /// /// Reduces per-query file-open overhead from O(n_genomes) to O(1) per partition. /// Column files are kept in place; packed files take priority when opening. - pub fn pack_matrices(&self) -> OKIResult<()> { - use obicompactvec::{pack_bit_matrix, pack_compact_int_matrix}; + /// + /// If `sparse` is set, presence matrices go one step further, from the + /// dense `.pbmx` form into `obicompactvec::PersistentSparseBitMatrix`'s + /// on-disk format (see `docmd/architecture/siblings.md`) — count + /// matrices are unaffected, sparse count matrices aren't implemented + /// (see the sparse-matrix design plan's "Explicitly deferred"). + pub fn pack_matrices(&self, sparse: bool) -> OKIResult<()> { + use obicompactvec::{pack_bit_matrix, pack_compact_int_matrix, pack_sparse_bit_matrix}; use obilayeredmap::meta::PartitionMeta; let n = self.n_partitions(); @@ -292,7 +298,13 @@ impl KmerIndex { let layer_dir = index_dir.join(format!("layer_{l}")); let presence_dir = layer_dir.join("presence"); let counts_dir = layer_dir.join("counts"); - if presence_dir.exists() { pack_bit_matrix(&presence_dir).map_err(OKIError::Io)?; } + if presence_dir.exists() { + if sparse { + pack_sparse_bit_matrix(&presence_dir).map_err(OKIError::Io)?; + } else { + pack_bit_matrix(&presence_dir).map_err(OKIError::Io)?; + } + } if counts_dir.exists() { pack_compact_int_matrix(&counts_dir).map_err(OKIError::Io)?; } } Ok(()) diff --git a/src/obikindex/src/merge.rs b/src/obikindex/src/merge.rs index 8d7840e9..5139d72e 100644 --- a/src/obikindex/src/merge.rs +++ b/src/obikindex/src/merge.rs @@ -249,7 +249,7 @@ impl KmerIndex { let pb = spinner("pack"); pb.set_message("consolidating column files …"); let dst2 = KmerIndex::open(output)?; - dst2.pack_matrices()?; + dst2.pack_matrices(false)?; pb.finish_and_clear(); rep.push(t.stop()); } diff --git a/src/obikindex/src/select.rs b/src/obikindex/src/select.rs index 4639d3b3..bc42a3a1 100644 --- a/src/obikindex/src/select.rs +++ b/src/obikindex/src/select.rs @@ -118,7 +118,7 @@ impl KmerIndex { self.meta.write(&self.root_path)?; let t_pack = Stage::start("pack"); - self.pack_matrices()?; + self.pack_matrices(false)?; rep.push(t_pack.stop()); Ok(()) } diff --git a/src/obikmer/src/cmd/pack/mod.rs b/src/obikmer/src/cmd/pack/mod.rs index 4f34c52d..10621772 100644 --- a/src/obikmer/src/cmd/pack/mod.rs +++ b/src/obikmer/src/cmd/pack/mod.rs @@ -9,6 +9,14 @@ use tracing::info; pub struct PackArgs { /// Index directory to pack pub index: PathBuf, + + /// Pack presence matrices into the sparse, deduplicated on-disk format + /// instead of the dense one — see `docmd/architecture/siblings.md`. + /// Smaller and faster for single-row access on real, sparse data; + /// column-oriented access (`--metric` distance matrices) is much + /// slower on the sparse format. + #[arg(long)] + pub sparse: bool, } pub fn run(args: PackArgs) { @@ -33,7 +41,7 @@ pub fn run(args: PackArgs) { let mut rep = Reporter::new(); let t = Stage::start("pack"); - idx.pack_matrices().unwrap_or_else(|e| { + idx.pack_matrices(args.sparse).unwrap_or_else(|e| { eprintln!("pack error: {e}"); std::process::exit(1); }); diff --git a/src/obikphylo/src/siblings/cache.rs b/src/obikphylo/src/siblings/cache.rs index 99933c6d..1ef1a672 100644 --- a/src/obikphylo/src/siblings/cache.rs +++ b/src/obikphylo/src/siblings/cache.rs @@ -1,6 +1,6 @@ use rayon::prelude::*; -use obicompactvec::{PersistentBitMatrix, PersistentCompactIntMatrix}; +use obicompactvec::{PersistentBitMatrix, PersistentCompactIntMatrix, PersistentSparseBitMatrix}; use obikpartitionner::KmerPartition; use obikseq::CanonicalKmer; use obilayeredmap::Layer; @@ -40,6 +40,13 @@ use super::{olm_to_ok, SiblingAnnex, INDEX_SUBDIR}; pub(super) enum Mat { Count(Layer), Presence(Layer), + /// Same role as `Presence`, over a layer packed by `pack --sparse` + /// (`PersistentSparseBitMatrix`) instead of the dense form — see + /// `docmd/architecture/siblings.md`'s sparse-matrix section. Every + /// `Mat` method below dispatches to the same `Layer` generic code + /// (`D: LayerData`) as `Presence`, so this arm is purely a storage + /// choice, not a behavioural difference. + SparsePresence(Layer), } impl Mat { @@ -47,6 +54,7 @@ impl Mat { match self { Mat::Count(l) => l.find_slot(kmer), Mat::Presence(l) => l.find_slot(kmer), + Mat::SparsePresence(l) => l.find_slot(kmer), } } @@ -60,6 +68,7 @@ impl Mat { match self { Mat::Count(l) => l.index_batch(kmers), Mat::Presence(l) => l.index_batch(kmers), + Mat::SparsePresence(l) => l.index_batch(kmers), } } @@ -74,6 +83,7 @@ impl Mat { match self { Mat::Count(l) => l.iter_minorants_batch(annex, batch_size), Mat::Presence(l) => l.iter_minorants_batch(annex, batch_size), + Mat::SparsePresence(l) => l.iter_minorants_batch(annex, batch_size), } } @@ -81,6 +91,7 @@ impl Mat { match self { Mat::Count(l) => l.n_cols(), Mat::Presence(l) => l.n_cols(), + Mat::SparsePresence(l) => l.n_cols(), } } @@ -100,6 +111,7 @@ impl Mat { pub(super) fn fill_sub_matrix_carries(&self, slots: &[usize], out: &mut [Vec]) { match self { Mat::Presence(l) => l.fill_sub_matrix(slots, out), + Mat::SparsePresence(l) => l.fill_sub_matrix(slots, out), Mat::Count(l) => { let mut counts: Vec> = out.iter().map(|_| Vec::new()).collect(); l.fill_sub_matrix(slots, &mut counts); @@ -145,8 +157,16 @@ impl PartitionCache { for l in 0..meta.n_layers { let layer_dir = index_dir.join(format!("layer_{l}")); let use_counts = with_counts && layer_dir.join("counts").exists(); + // `pack --sparse` converts a layer's `presence/` directory + // in place, leaving `is_multi.prsb` as the on-disk marker + // that distinguishes the sparse format + // (`obicompactvec::PersistentSparseBitMatrix`) from the + // dense one — see `bitmatrix/sparse.rs`'s module docs. + let is_sparse = layer_dir.join("presence").join("is_multi.prsb").exists(); let mat = if use_counts { Layer::::open(&layer_dir, &meta.mode).ok().map(Mat::Count) + } else if is_sparse { + Layer::::open(&layer_dir, &meta.mode).ok().map(Mat::SparsePresence) } else { Layer::::open(&layer_dir, &meta.mode).ok().map(Mat::Presence) }; diff --git a/src/obikphylo/src/siblings/tests.rs b/src/obikphylo/src/siblings/tests.rs index 513316ca..9979678c 100644 --- a/src/obikphylo/src/siblings/tests.rs +++ b/src/obikphylo/src/siblings/tests.rs @@ -148,6 +148,47 @@ fn sibling_annex_one_sibling_each() { assert!(!b.is_minorant(), "g2's stored minorant flag should not be set"); } +#[test] +fn sibling_annex_works_after_pack_sparse() { + // Same fixture as `sibling_annex_one_sibling_each`, but with + // `pack --sparse` (`KmerIndex::pack_matrices(true)`) run on the merged + // index before building the sibling annex — proves `PartitionCache`'s + // sparse-detection (`Mat::SparsePresence`, gated on `presence/ + // is_multi.prsb`) and the generic `Layer` methods it relies on + // actually round-trip through the real build pipeline, not just the + // unit-level `Layer` tests in + // `obilayeredmap`. + let dir = tempdir().unwrap(); + let g1 = build_single_genome_index(dir.path(), "g1", b"AACCGCTTAAG"); + let g2 = build_single_genome_index(dir.path(), "g2", b"AACCGGTTAAG"); + let merged = merge_two(dir.path(), &g1, &g2); + merged.pack_matrices(true).expect("pack_matrices(sparse)"); + + let index_dir = merged.partition().part_dir(0).join(INDEX_SUBDIR); + assert!( + index_dir.join("layer_0").join("presence").join("is_multi.prsb").exists(), + "pack_matrices(true) must leave the sparse marker file behind" + ); + + merged.build_sibling_annex().expect("build_sibling_annex"); + + let g1_kmer = canonical(b"AACCGCTTAAG"); + let g2_kmer = canonical(b"AACCGGTTAAG"); + let expected_mask = FamilyMask::EMPTY.with(1).with(2); + + let a = annex_info_for(&merged, g1_kmer); + assert_eq!(a.bits(), expected_mask.bits(), "AACCGCTTAAG"); + assert_eq!(a.siblings(), 1); + assert!(is_minorant(g1_kmer, a, K), "g1 should be the minorant"); + assert!(a.is_minorant(), "g1's stored minorant flag should be set at build time"); + + let b = annex_info_for(&merged, g2_kmer); + assert_eq!(b.bits(), expected_mask.bits(), "AACCGGTTAAG"); + assert_eq!(b.siblings(), 1); + assert!(!is_minorant(g2_kmer, b, K), "g2 should not be the minorant"); + assert!(!b.is_minorant(), "g2's stored minorant flag should not be set"); +} + #[test] fn sibling_annex_zero_siblings_when_identical_and_no_variant_exists() { // Same k-mer in both genomes, no other genome around to carry a diff --git a/src/obilayeredmap/src/layer.rs b/src/obilayeredmap/src/layer.rs index d50fd7b8..4aa86365 100644 --- a/src/obilayeredmap/src/layer.rs +++ b/src/obilayeredmap/src/layer.rs @@ -3,8 +3,10 @@ use std::fs; use std::path::Path; use obicompactvec::{ + BinaryMatrix, PersistentBitMatrix, PersistentBitMatrixBuilder, PersistentCompactIntMatrix, PersistentCompactIntMatrixBuilder, + PersistentSparseBitMatrix, }; use obikseq::CanonicalKmer; use obiskio::{UnitigFileReader, UnitigFileWriter}; @@ -47,6 +49,14 @@ impl LayerData for PersistentBitMatrix { fn read(&self, slot: usize) -> Box<[bool]> { self.row(slot) } } +impl LayerData for PersistentSparseBitMatrix { + type Item = Box<[bool]>; + fn open(layer_dir: &Path) -> OLMResult { + PersistentSparseBitMatrix::open(&layer_dir.join(PRESENCE_DIR)).map_err(OLMError::Io) + } + fn read(&self, slot: usize) -> Box<[bool]> { self.row(slot) } +} + // ── Structures ──────────────────────────────────────────────────────────────── pub struct Layer { @@ -223,15 +233,18 @@ impl Layer { // ── Mode 3 — presence/absence matrix ───────────────────────────────────────── -impl Layer { - pub fn append_genome_column( - layer_dir: &Path, - value_of: impl Fn(usize) -> bool, - ) -> OLMResult<()> { - PersistentBitMatrix::append_column(&layer_dir.join(PRESENCE_DIR), value_of) - .map_err(OLMError::Io) - } - +// ── Mode 3 — presence/absence matrix, generic over dense/sparse storage ────── +// +// `n_cols`/`sub_matrix`/`fill_sub_matrix` are identical in shape for any +// `D` satisfying `BinaryMatrix` — genuinely generic (not just two +// near-identical impl blocks) so `Layer` gets +// them for free, matching `Layer` at the same call +// sites (see `obikphylo::siblings::cache::Mat`). Construction +// (`append_genome_column`/`build_presence`) stays `PersistentBitMatrix`-only +// below: sparse matrices aren't built column-by-column, they're built +// row-by-row from an already-built dense layer +// (`PersistentSparseBitMatrixBuilder::build_from_dense`). +impl> + BinaryMatrix> Layer { /// Number of genome columns in this layer's presence matrix — see /// `PersistentBitMatrix::n_cols`'s docs for the `Implicit` mono-genome /// special case (always reports `1`, regardless of the index's real @@ -244,7 +257,9 @@ impl Layer { /// /// Returns a column-first `Vec>`: one inner `Vec` per genome /// column, containing the presence bits for the requested slots in order. - /// Column access is sequential for cache efficiency. + /// Column access is sequential for cache efficiency on the dense + /// storage — on sparse storage it's a naive row-by-row decode, see + /// `docmd/architecture/siblings.md`'s sparse-matrix section. pub fn sub_matrix(&self, slots: &[usize]) -> Vec> { self.data.sub_matrix(slots) } @@ -258,6 +273,16 @@ impl Layer { pub fn fill_sub_matrix(&self, slots: &[usize], out: &mut [Vec]) { self.data.fill_sub_matrix(slots, out) } +} + +impl Layer { + pub fn append_genome_column( + layer_dir: &Path, + value_of: impl Fn(usize) -> bool, + ) -> OLMResult<()> { + PersistentBitMatrix::append_column(&layer_dir.join(PRESENCE_DIR), value_of) + .map_err(OLMError::Io) + } pub fn build_presence( out_dir: &Path, diff --git a/src/obilayeredmap/src/tests/layer.rs b/src/obilayeredmap/src/tests/layer.rs index 3ae37376..033e819d 100644 --- a/src/obilayeredmap/src/tests/layer.rs +++ b/src/obilayeredmap/src/tests/layer.rs @@ -1,5 +1,6 @@ use super::*; -use obikseq::{set_k, Kmer, Sequence as _, Unitig}; +use obicompactvec::PersistentSparseBitMatrixBuilder; +use obikseq::{set_k, Unitig}; use obiskio::DEFAULT_BLOCK_BITS; use tempfile::tempdir; @@ -42,6 +43,61 @@ fn canonical_kmer_iter_matches_reader() { assert_eq!(from_iter, from_reader, "CanonicalKmerIter and UnitigFileReader disagree"); } +// ── Generic `Layer` over dense vs. sparse presence storage ─────────────── + +#[test] +fn presence_layer_generic_over_sparse_matches_dense() { + // k=4, matching every other test in this crate's binary: `K`/`M` are + // process-wide (not thread-local) in test builds — see + // `obikseq::params` — so a test using a different k here would race + // against `map.rs`'s k=4 tests running concurrently in the same + // binary. These three sequences were chosen (randomised search) to + // have no shared canonical 4-mer among them, required for `ptr_hash` + // MPHF construction (duplicate keys make it panic). + set_k(4); + let dir = tempdir().unwrap(); + write_unitigs(dir.path(), &[b"AAATCTA", b"CTTCGCC", b"TGATACG"]); + let n_genomes = 3; + let mode = IndexMode::Exact; + + // Deterministic, arbitrary presence function — genome `g` carries a + // kmer iff `(kmer's raw bits + g)` is even. Doesn't need to be + // biologically meaningful, just the same on both sides of the + // dense/sparse comparison below. + Layer::::build_presence(dir.path(), DEFAULT_BLOCK_BITS, &mode, n_genomes, |kmer, g| { + (kmer.raw().wrapping_add(g as u64)) % 2 == 0 + }).unwrap(); + + let dense_layer = Layer::::open(dir.path(), &mode).unwrap(); + assert!(dense_layer.n_cols() >= 1); + + // Build the sparse form directly into the same `presence/` dir the + // dense form already lives in (matches how `pack --sparse` will work: + // same directory, distinct filenames — see + // `docmd/architecture/siblings.md`'s sparse-matrix section). + let dense_matrix = obicompactvec::PersistentBitMatrix::open(dir.path()).unwrap(); + PersistentSparseBitMatrixBuilder::build_from_dense(&dense_matrix, &dir.path().join(PRESENCE_DIR)) + .unwrap() + .close() + .unwrap(); + + let sparse_layer = Layer::::open(dir.path(), &mode).unwrap(); + + // Same generic methods, same results, different concrete `D`. + assert_eq!(dense_layer.n_cols(), sparse_layer.n_cols()); + let n = dense_layer.n(); + assert_eq!(n, sparse_layer.n()); + let slots: Vec = (0..n).collect(); + assert_eq!(dense_layer.sub_matrix(&slots), sparse_layer.sub_matrix(&slots)); + + // The matrix-agnostic, MPHF-only surface (from the generic + // `impl Layer`) must also agree: same kmer set, + // looked up through either concrete `D`. + for kmer in all_canonical_kmers(dir.path()) { + assert_eq!(dense_layer.find_slot(kmer), sparse_layer.find_slot(kmer)); + } +} + #[test] fn enumerate_kmers_is_stable_across_calls() { set_k(4);