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.
This commit is contained in:
Eric Coissac
2026-08-16 21:51:02 +02:00
parent 50f4820cb9
commit 3ba26b3dc1
14 changed files with 279 additions and 24 deletions
+1
View File
@@ -16,6 +16,7 @@ benchmark/simulated_data
benchmark/specimen_index_presence benchmark/specimen_index_presence
benchmark/specimen_index_count benchmark/specimen_index_count
benchmark/global_index_presence benchmark/global_index_presence
benchmark/global_index_presence_sav
benchmark/all_specific benchmark/all_specific
benchmark/global_index_count benchmark/global_index_count
benchmark/stats benchmark/stats
+56
View File
@@ -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 et al. subset-hierarchy compression (only the exact-duplicate special case
is built), a sparse `PersistentCompactIntMatrix` (count matrices), and is built), a sparse `PersistentCompactIntMatrix` (count matrices), and
BRWT-style column-correlation exploitation. 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<D>`**: `obilayeredmap::Layer<D>`'s presence-only methods
(`n_cols`, `sub_matrix`, `fill_sub_matrix`) are generic over any
`D: LayerData<Item = Box<[bool]>> + 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<PersistentBitMatrix>`/`Layer<PersistentSparseBitMatrix>` 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<PersistentSparseBitMatrix>)`, 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<D>` 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.
+1 -1
View File
@@ -22,7 +22,7 @@ mod sparse;
pub use builder::PersistentBitMatrixBuilder; pub use builder::PersistentBitMatrixBuilder;
pub use packed::pack_bit_matrix; pub use packed::pack_bit_matrix;
pub use persistent::PersistentBitMatrix; 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}; pub(crate) use pairwise::{pairwise_matrix, pairwise2_matrix};
+31 -3
View File
@@ -5,7 +5,7 @@
//! used by any production code path yet — a new type, not a replacement. //! used by any production code path yet — a new type, not a replacement.
//! //!
//! On-disk layout (a directory, mirroring [`super::PersistentBitMatrix`]'s //! 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 //! `is_multi.prsb` (rank-capable flag: singleton row vs. multi-genome
//! row), `singleton.pfiv` (genome index, one entry per singleton row, //! row), `singleton.pfiv` (genome index, one entry per singleton row,
//! `ceil(log2(n_cols))` bits each), `multi.pfiv` (`dict_id`, one entry per //! `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 crate::rankselect::{PersistentRankSelectBitVec, PersistentRankSelectBitVecBuilder};
use super::PersistentBitMatrix; use super::PersistentBitMatrix;
use super::packed::PackedBitMatrix;
fn is_multi_path(dir: &Path) -> PathBuf { dir.join("is_multi.prsb") } fn is_multi_path(dir: &Path) -> PathBuf { dir.join("is_multi.prsb") }
fn singleton_path(dir: &Path) -> PathBuf { dir.join("singleton.pfiv") } fn singleton_path(dir: &Path) -> PathBuf { dir.join("singleton.pfiv") }
fn multi_path(dir: &Path) -> PathBuf { dir.join("multi.pfiv") } fn multi_path(dir: &Path) -> PathBuf { dir.join("multi.pfiv") }
fn dict_offsets_base(dir: &Path) -> PathBuf { dir.join("dict_offsets") } fn dict_offsets_base(dir: &Path) -> PathBuf { dir.join("dict_offsets") }
fn dict_values_path(dir: &Path) -> PathBuf { dir.join("dict_values.bin") } 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 { struct SparseMeta {
n: usize, n: usize,
@@ -53,7 +54,7 @@ impl SparseMeta {
fn load(dir: &Path) -> io::Result<Self> { fn load(dir: &Path) -> io::Result<Self> {
let s = fs::read_to_string(meta_path(dir))?; let s = fs::read_to_string(meta_path(dir))?;
let get = |name: &str| { 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 { Ok(Self {
n: get("n")?, n: get("n")?,
@@ -370,3 +371,30 @@ impl PersistentSparseBitMatrixBuilder {
Ok(builder) 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(())
}
+1 -1
View File
@@ -19,7 +19,7 @@ pub use bitvec::{BitIter, PersistentBitVec, PersistentBitVecBuilder};
pub use fixedintvec::{PersistentFixedIntVec, PersistentFixedIntVecBuilder, bit_width_for_range}; pub use fixedintvec::{PersistentFixedIntVec, PersistentFixedIntVecBuilder, bit_width_for_range};
pub use rankselect::{PersistentRankSelectBitVec, PersistentRankSelectBitVecBuilder}; pub use rankselect::{PersistentRankSelectBitVec, PersistentRankSelectBitVecBuilder};
pub use eliasfano::{EliasFano, EliasFanoBuilder}; 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 builder::PersistentCompactIntVecBuilder;
pub use colgroup::{ColGroup, FilterMask, MatrixGroupOps, eval_filter_mask}; pub use colgroup::{ColGroup, FilterMask, MatrixGroupOps, eval_filter_mask};
pub use intmatrix::{PersistentCompactIntMatrix, PersistentCompactIntMatrixBuilder, pack_compact_int_matrix}; pub use intmatrix::{PersistentCompactIntMatrix, PersistentCompactIntMatrixBuilder, pack_compact_int_matrix};
+8
View File
@@ -26,6 +26,14 @@ pub trait BinaryMatrix {
fn fill_sub_matrix(&self, slots: &[usize], out: &mut [Vec<bool>]); fn fill_sub_matrix(&self, slots: &[usize], out: &mut [Vec<bool>]);
/// Per-genome k-mer totals. /// Per-genome k-mer totals.
fn count_ones(&self) -> Array1<u64>; fn count_ones(&self) -> Array1<u64>;
/// 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<Vec<bool>> {
let mut out: Vec<Vec<bool>> = (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 /// Convert a Jaccard distance matrix (`1 - J`) into a Mash distance matrix, per
+16 -4
View File
@@ -122,7 +122,7 @@ impl KmerIndex {
fs::File::create(output.join(SENTINEL_INDEXED)).map_err(OKIError::Io)?; fs::File::create(output.join(SENTINEL_INDEXED)).map_err(OKIError::Io)?;
let idx = KmerIndex::open(output)?; let idx = KmerIndex::open(output)?;
let t_pack = Stage::start("pack"); let t_pack = Stage::start("pack");
idx.pack_matrices()?; idx.pack_matrices(false)?;
rep.push(t_pack.stop()); rep.push(t_pack.stop());
Ok(idx) Ok(idx)
} }
@@ -274,8 +274,14 @@ impl KmerIndex {
/// ///
/// Reduces per-query file-open overhead from O(n_genomes) to O(1) per partition. /// 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. /// 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; use obilayeredmap::meta::PartitionMeta;
let n = self.n_partitions(); let n = self.n_partitions();
@@ -292,7 +298,13 @@ impl KmerIndex {
let layer_dir = index_dir.join(format!("layer_{l}")); let layer_dir = index_dir.join(format!("layer_{l}"));
let presence_dir = layer_dir.join("presence"); let presence_dir = layer_dir.join("presence");
let counts_dir = layer_dir.join("counts"); 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)?; } if counts_dir.exists() { pack_compact_int_matrix(&counts_dir).map_err(OKIError::Io)?; }
} }
Ok(()) Ok(())
+1 -1
View File
@@ -249,7 +249,7 @@ impl KmerIndex {
let pb = spinner("pack"); let pb = spinner("pack");
pb.set_message("consolidating column files …"); pb.set_message("consolidating column files …");
let dst2 = KmerIndex::open(output)?; let dst2 = KmerIndex::open(output)?;
dst2.pack_matrices()?; dst2.pack_matrices(false)?;
pb.finish_and_clear(); pb.finish_and_clear();
rep.push(t.stop()); rep.push(t.stop());
} }
+1 -1
View File
@@ -118,7 +118,7 @@ impl KmerIndex {
self.meta.write(&self.root_path)?; self.meta.write(&self.root_path)?;
let t_pack = Stage::start("pack"); let t_pack = Stage::start("pack");
self.pack_matrices()?; self.pack_matrices(false)?;
rep.push(t_pack.stop()); rep.push(t_pack.stop());
Ok(()) Ok(())
} }
+9 -1
View File
@@ -9,6 +9,14 @@ use tracing::info;
pub struct PackArgs { pub struct PackArgs {
/// Index directory to pack /// Index directory to pack
pub index: PathBuf, 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) { pub fn run(args: PackArgs) {
@@ -33,7 +41,7 @@ pub fn run(args: PackArgs) {
let mut rep = Reporter::new(); let mut rep = Reporter::new();
let t = Stage::start("pack"); 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}"); eprintln!("pack error: {e}");
std::process::exit(1); std::process::exit(1);
}); });
+21 -1
View File
@@ -1,6 +1,6 @@
use rayon::prelude::*; use rayon::prelude::*;
use obicompactvec::{PersistentBitMatrix, PersistentCompactIntMatrix}; use obicompactvec::{PersistentBitMatrix, PersistentCompactIntMatrix, PersistentSparseBitMatrix};
use obikpartitionner::KmerPartition; use obikpartitionner::KmerPartition;
use obikseq::CanonicalKmer; use obikseq::CanonicalKmer;
use obilayeredmap::Layer; use obilayeredmap::Layer;
@@ -40,6 +40,13 @@ use super::{olm_to_ok, SiblingAnnex, INDEX_SUBDIR};
pub(super) enum Mat { pub(super) enum Mat {
Count(Layer<PersistentCompactIntMatrix>), Count(Layer<PersistentCompactIntMatrix>),
Presence(Layer<PersistentBitMatrix>), Presence(Layer<PersistentBitMatrix>),
/// 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<D>` generic code
/// (`D: LayerData`) as `Presence`, so this arm is purely a storage
/// choice, not a behavioural difference.
SparsePresence(Layer<PersistentSparseBitMatrix>),
} }
impl Mat { impl Mat {
@@ -47,6 +54,7 @@ impl Mat {
match self { match self {
Mat::Count(l) => l.find_slot(kmer), Mat::Count(l) => l.find_slot(kmer),
Mat::Presence(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 { match self {
Mat::Count(l) => l.index_batch(kmers), Mat::Count(l) => l.index_batch(kmers),
Mat::Presence(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 { match self {
Mat::Count(l) => l.iter_minorants_batch(annex, batch_size), Mat::Count(l) => l.iter_minorants_batch(annex, batch_size),
Mat::Presence(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 { match self {
Mat::Count(l) => l.n_cols(), Mat::Count(l) => l.n_cols(),
Mat::Presence(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<bool>]) { pub(super) fn fill_sub_matrix_carries(&self, slots: &[usize], out: &mut [Vec<bool>]) {
match self { match self {
Mat::Presence(l) => l.fill_sub_matrix(slots, out), Mat::Presence(l) => l.fill_sub_matrix(slots, out),
Mat::SparsePresence(l) => l.fill_sub_matrix(slots, out),
Mat::Count(l) => { Mat::Count(l) => {
let mut counts: Vec<Vec<u32>> = out.iter().map(|_| Vec::new()).collect(); let mut counts: Vec<Vec<u32>> = out.iter().map(|_| Vec::new()).collect();
l.fill_sub_matrix(slots, &mut counts); l.fill_sub_matrix(slots, &mut counts);
@@ -145,8 +157,16 @@ impl PartitionCache {
for l in 0..meta.n_layers { for l in 0..meta.n_layers {
let layer_dir = index_dir.join(format!("layer_{l}")); let layer_dir = index_dir.join(format!("layer_{l}"));
let use_counts = with_counts && layer_dir.join("counts").exists(); 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 { let mat = if use_counts {
Layer::<PersistentCompactIntMatrix>::open(&layer_dir, &meta.mode).ok().map(Mat::Count) Layer::<PersistentCompactIntMatrix>::open(&layer_dir, &meta.mode).ok().map(Mat::Count)
} else if is_sparse {
Layer::<PersistentSparseBitMatrix>::open(&layer_dir, &meta.mode).ok().map(Mat::SparsePresence)
} else { } else {
Layer::<PersistentBitMatrix>::open(&layer_dir, &meta.mode).ok().map(Mat::Presence) Layer::<PersistentBitMatrix>::open(&layer_dir, &meta.mode).ok().map(Mat::Presence)
}; };
+41
View File
@@ -148,6 +148,47 @@ fn sibling_annex_one_sibling_each() {
assert!(!b.is_minorant(), "g2's stored minorant flag should not be set"); 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<D>` methods it relies on
// actually round-trip through the real build pipeline, not just the
// unit-level `Layer<PersistentSparseBitMatrix>` 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] #[test]
fn sibling_annex_zero_siblings_when_identical_and_no_variant_exists() { fn sibling_annex_zero_siblings_when_identical_and_no_variant_exists() {
// Same k-mer in both genomes, no other genome around to carry a // Same k-mer in both genomes, no other genome around to carry a
+35 -10
View File
@@ -3,8 +3,10 @@ use std::fs;
use std::path::Path; use std::path::Path;
use obicompactvec::{ use obicompactvec::{
BinaryMatrix,
PersistentBitMatrix, PersistentBitMatrixBuilder, PersistentBitMatrix, PersistentBitMatrixBuilder,
PersistentCompactIntMatrix, PersistentCompactIntMatrixBuilder, PersistentCompactIntMatrix, PersistentCompactIntMatrixBuilder,
PersistentSparseBitMatrix,
}; };
use obikseq::CanonicalKmer; use obikseq::CanonicalKmer;
use obiskio::{UnitigFileReader, UnitigFileWriter}; use obiskio::{UnitigFileReader, UnitigFileWriter};
@@ -47,6 +49,14 @@ impl LayerData for PersistentBitMatrix {
fn read(&self, slot: usize) -> Box<[bool]> { self.row(slot) } fn read(&self, slot: usize) -> Box<[bool]> { self.row(slot) }
} }
impl LayerData for PersistentSparseBitMatrix {
type Item = Box<[bool]>;
fn open(layer_dir: &Path) -> OLMResult<Self> {
PersistentSparseBitMatrix::open(&layer_dir.join(PRESENCE_DIR)).map_err(OLMError::Io)
}
fn read(&self, slot: usize) -> Box<[bool]> { self.row(slot) }
}
// ── Structures ──────────────────────────────────────────────────────────────── // ── Structures ────────────────────────────────────────────────────────────────
pub struct Layer<D: LayerData = ()> { pub struct Layer<D: LayerData = ()> {
@@ -223,15 +233,18 @@ impl Layer<PersistentCompactIntMatrix> {
// ── Mode 3 — presence/absence matrix ───────────────────────────────────────── // ── Mode 3 — presence/absence matrix ─────────────────────────────────────────
impl Layer<PersistentBitMatrix> { // ── Mode 3 — presence/absence matrix, generic over dense/sparse storage ──────
pub fn append_genome_column( //
layer_dir: &Path, // `n_cols`/`sub_matrix`/`fill_sub_matrix` are identical in shape for any
value_of: impl Fn(usize) -> bool, // `D` satisfying `BinaryMatrix` — genuinely generic (not just two
) -> OLMResult<()> { // near-identical impl blocks) so `Layer<PersistentSparseBitMatrix>` gets
PersistentBitMatrix::append_column(&layer_dir.join(PRESENCE_DIR), value_of) // them for free, matching `Layer<PersistentBitMatrix>` at the same call
.map_err(OLMError::Io) // 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<D: LayerData<Item = Box<[bool]>> + BinaryMatrix> Layer<D> {
/// Number of genome columns in this layer's presence matrix — see /// Number of genome columns in this layer's presence matrix — see
/// `PersistentBitMatrix::n_cols`'s docs for the `Implicit` mono-genome /// `PersistentBitMatrix::n_cols`'s docs for the `Implicit` mono-genome
/// special case (always reports `1`, regardless of the index's real /// special case (always reports `1`, regardless of the index's real
@@ -244,7 +257,9 @@ impl Layer<PersistentBitMatrix> {
/// ///
/// Returns a column-first `Vec<Vec<bool>>`: one inner `Vec` per genome /// Returns a column-first `Vec<Vec<bool>>`: one inner `Vec` per genome
/// column, containing the presence bits for the requested slots in order. /// 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<Vec<bool>> { pub fn sub_matrix(&self, slots: &[usize]) -> Vec<Vec<bool>> {
self.data.sub_matrix(slots) self.data.sub_matrix(slots)
} }
@@ -258,6 +273,16 @@ impl Layer<PersistentBitMatrix> {
pub fn fill_sub_matrix(&self, slots: &[usize], out: &mut [Vec<bool>]) { pub fn fill_sub_matrix(&self, slots: &[usize], out: &mut [Vec<bool>]) {
self.data.fill_sub_matrix(slots, out) self.data.fill_sub_matrix(slots, out)
} }
}
impl Layer<PersistentBitMatrix> {
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( pub fn build_presence(
out_dir: &Path, out_dir: &Path,
+57 -1
View File
@@ -1,5 +1,6 @@
use super::*; use super::*;
use obikseq::{set_k, Kmer, Sequence as _, Unitig}; use obicompactvec::PersistentSparseBitMatrixBuilder;
use obikseq::{set_k, Unitig};
use obiskio::DEFAULT_BLOCK_BITS; use obiskio::DEFAULT_BLOCK_BITS;
use tempfile::tempdir; use tempfile::tempdir;
@@ -42,6 +43,61 @@ fn canonical_kmer_iter_matches_reader() {
assert_eq!(from_iter, from_reader, "CanonicalKmerIter and UnitigFileReader disagree"); assert_eq!(from_iter, from_reader, "CanonicalKmerIter and UnitigFileReader disagree");
} }
// ── Generic `Layer<D>` 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::<PersistentBitMatrix>::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::<PersistentBitMatrix>::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::<PersistentSparseBitMatrix>::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<usize> = (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<D: LayerData> Layer<D>`) 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] #[test]
fn enumerate_kmers_is_stable_across_calls() { fn enumerate_kmers_is_stable_across_calls() {
set_k(4); set_k(4);