feat: add sparse compact integer matrix implementation and utilities

Introduces `PersistentSparseCompactIntMatrix`, a row-major sparse integer matrix combining a bitmatrix support layer with dual value streams for singleton and multi partitions. Adds builder and packing utilities to convert dense count matrices into the new on-disk format, featuring idempotency guards and crash-safe persistence workflows. Exposes core accessors, iterators, and column weight calculations while adjusting internal visibility for companion modules. Includes unit tests and diagnostic examples to validate equivalence with dense representations and measure compaction metrics.
This commit is contained in:
Eric Coissac
2026-08-26 14:05:54 +02:00
parent 294f132a0a
commit b7a8b5e6cf
9 changed files with 1051 additions and 2 deletions
@@ -0,0 +1,212 @@
//! Diagnostic: build a `PersistentSparseCompactIntMatrix` from a real dense
//! `PersistentCompactIntMatrix` (batched `build_from_dense`), verify the
//! result is cell-for-cell identical on an actual count matrix, and report
//! the on-disk compaction ratio, value-distribution stats, and effective
//! bits/value.
//!
//! Usage: cargo run --release --example compare_sparse_dense_count -p obicompactvec -- <layer_dir>
//! (layer_dir is the directory containing `counts/matrix.pcmx` or
//! `counts/col_*.pciv`, e.g.
//! `benchmark/global_index_count/partitions/part_00003/index/layer_1`)
use std::error::Error;
use std::path::Path;
use std::time::Instant;
use obicompactvec::{PersistentCompactIntMatrix, PersistentSparseCompactIntMatrix, PersistentSparseCompactIntMatrixBuilder};
fn main() -> Result<(), Box<dyn Error>> {
let layer_dir = std::env::args().nth(1).expect("usage: compare_sparse_dense_count <layer_dir>");
let layer_dir = Path::new(&layer_dir);
let t0 = Instant::now();
let dense = PersistentCompactIntMatrix::open(layer_dir)?;
println!("dense ouverte en {:?} ({} lignes x {} colonnes)", t0.elapsed(), dense.n(), dense.n_cols());
let out_dir = tempfile::tempdir()?;
let t0 = Instant::now();
let sparse = PersistentSparseCompactIntMatrixBuilder::build_from_dense(&dense, out_dir.path())?.finish()?;
println!("build_from_dense: {:?}", t0.elapsed());
compare(&dense, &sparse)?;
content_stats(&dense);
compaction_stats(layer_dir, out_dir.path(), &dense, &sparse)?;
Ok(())
}
fn compare(dense: &PersistentCompactIntMatrix, sparse: &PersistentSparseCompactIntMatrix) -> Result<(), Box<dyn Error>> {
let n = dense.n();
let n_cols = dense.n_cols();
assert_eq!(n, sparse.n(), "n mismatch");
assert_eq!(n_cols, sparse.n_cols(), "n_cols mismatch");
let mut dense_row = vec![0u32; n_cols];
let mut sparse_row = vec![0u32; n_cols];
let mut total_cells = 0usize;
let mut mismatched_cells = 0usize;
let mut first_mismatch = None;
let t0 = Instant::now();
for slot in 0..n {
dense.fill_row(slot, &mut dense_row);
sparse.fill_row(slot, &mut sparse_row);
for c in 0..n_cols {
total_cells += 1;
if dense_row[c] != sparse_row[c] {
mismatched_cells += 1;
if first_mismatch.is_none() {
first_mismatch = Some((slot, c, dense_row[c], sparse_row[c]));
}
}
}
}
println!("comparaison case-a-case: {:?} ({total_cells} cellules)", t0.elapsed());
if let Some((slot, col, d, s)) = first_mismatch {
eprintln!("PREMIER MISMATCH: slot={slot} col={col} dense={d} sparse={s}");
println!("MISMATCHES: {mismatched_cells} cellules differentes sur {total_cells}");
return Err("dense et sparse divergent".into());
}
println!("OK: toutes les {total_cells} cellules sont identiques entre dense et sparse.");
Ok(())
}
/// Content stats computed straight from the dense matrix (source of truth
/// for what's actually stored) — sparsity, singleton-row fraction, and the
/// value-magnitude buckets this crate's overflow encoding is tuned around
/// (< 127, < 255, >= 255).
fn content_stats(dense: &PersistentCompactIntMatrix) {
let n = dense.n();
let n_cols = dense.n_cols();
let total_cells = n as u64 * n_cols as u64;
let mut nonzero_cells = 0u64;
let mut singleton_rows = 0u64;
let mut sum: u128 = 0;
let mut max_value = 0u32;
let mut under_127 = 0u64;
let mut under_255 = 0u64;
let mut overflow = 0u64;
let mut row = vec![0u32; n_cols];
for slot in 0..n {
dense.fill_row(slot, &mut row);
let mut row_nonzero = 0u32;
for &v in &row {
if v == 0 {
continue;
}
row_nonzero += 1;
nonzero_cells += 1;
sum += v as u128;
max_value = max_value.max(v);
if v < 127 {
under_127 += 1;
}
if v < 255 {
under_255 += 1;
} else {
overflow += 1;
}
}
if row_nonzero == 1 {
singleton_rows += 1;
}
}
println!("\n--- stats de contenu (source: dense) ---");
println!("cellules totales: {total_cells}");
println!(
"cellules non-nulles: {nonzero_cells} ({:.3}% du total)",
100.0 * nonzero_cells as f64 / total_cells as f64
);
println!(
"lignes singleton: {singleton_rows} / {n} ({:.3}%)",
100.0 * singleton_rows as f64 / n as f64
);
if nonzero_cells > 0 {
println!("valeur moyenne (non-nulles): {:.2}", sum as f64 / nonzero_cells as f64);
println!("valeur max: {max_value}");
println!(
"valeurs < 127: {under_127} ({:.3}% des non-nulles)",
100.0 * under_127 as f64 / nonzero_cells as f64
);
println!(
"valeurs < 255: {under_255} ({:.3}% des non-nulles)",
100.0 * under_255 as f64 / nonzero_cells as f64
);
println!(
"valeurs >= 255 (overflow): {overflow} ({:.3}% des non-nulles)",
100.0 * overflow as f64 / nonzero_cells as f64
);
}
}
/// Directory size in bytes — sums every regular file, one level deep
/// (matches both the dense `counts/` layout and the sparse builder's flat
/// output directory, neither of which nests further).
fn dir_size(dir: &Path) -> std::io::Result<u64> {
let mut total = 0u64;
for entry in std::fs::read_dir(dir)? {
let entry = entry?;
if entry.file_type()?.is_file() {
total += entry.metadata()?.len();
}
}
Ok(total)
}
fn compaction_stats(
dense_layer_dir: &Path,
sparse_dir: &Path,
dense: &PersistentCompactIntMatrix,
_sparse: &PersistentSparseCompactIntMatrix,
) -> Result<(), Box<dyn Error>> {
let dense_size = dir_size(&dense_layer_dir.join("counts"))?;
let sparse_size = dir_size(sparse_dir)?;
let n = dense.n();
let n_cols = dense.n_cols();
let nonzero_cells: u64 = {
let mut row = vec![0u32; n_cols];
let mut count = 0u64;
for slot in 0..n {
dense.fill_row(slot, &mut row);
count += row.iter().filter(|&&v| v != 0).count() as u64;
}
count
};
println!("\n--- compaction ---");
println!("taille dense (counts/): {dense_size} octets");
println!("taille sparse: {sparse_size} octets");
println!(
"ratio sparse/dense: {:.4} ({:.2}% de la taille dense)",
sparse_size as f64 / dense_size as f64,
100.0 * sparse_size as f64 / dense_size as f64
);
if nonzero_cells > 0 {
println!(
"bits/valeur (dense, sur cellules non-nulles): {:.3}",
dense_size as f64 * 8.0 / nonzero_cells as f64
);
println!(
"bits/valeur (sparse, sur cellules non-nulles): {:.3}",
sparse_size as f64 * 8.0 / nonzero_cells as f64
);
}
let total_cells = n as u64 * n_cols as u64;
println!(
"bits/valeur (dense, sur toutes les cellules): {:.3}",
dense_size as f64 * 8.0 / total_cells as f64
);
println!(
"bits/valeur (sparse, sur toutes les cellules): {:.3}",
sparse_size as f64 * 8.0 / total_cells as f64
);
Ok(())
}
@@ -0,0 +1,169 @@
//! Diagnostic: build a `PersistentSparseBitMatrix` from a real dense
//! `PersistentBitMatrix` (batched `build_from_dense`), verify the result is
//! bit-for-bit identical on an actual presence matrix, and report the
//! on-disk compaction ratio, content stats, and effective bits/set-bit.
//! Same shape as `compare_sparse_dense_count.rs`, for the presence (bit)
//! side instead of counts.
//!
//! Usage: cargo run --release --example compare_sparse_dense_presence -p obicompactvec -- <layer_dir>
//! (layer_dir is the directory containing `presence/matrix.pbmx`, e.g.
//! `benchmark/tmp/bacteria/partitions/part_00000/index/layer_1`)
use std::error::Error;
use std::path::Path;
use std::time::Instant;
use obicompactvec::{PersistentBitMatrix, PersistentSparseBitMatrix, PersistentSparseBitMatrixBuilder};
fn main() -> Result<(), Box<dyn Error>> {
let layer_dir = std::env::args().nth(1).expect("usage: compare_sparse_dense_presence <layer_dir>");
let layer_dir = Path::new(&layer_dir);
let t0 = Instant::now();
let dense = PersistentBitMatrix::open(layer_dir)?;
println!("dense ouverte en {:?} ({} lignes x {} colonnes)", t0.elapsed(), dense.n(), dense.n_cols());
let out_dir = tempfile::tempdir()?;
let t0 = Instant::now();
let sparse = PersistentSparseBitMatrixBuilder::build_from_dense(&dense, out_dir.path())?.finish()?;
println!("build_from_dense: {:?}", t0.elapsed());
compare(&dense, &sparse)?;
content_stats(&dense);
compaction_stats(layer_dir, out_dir.path(), &dense)?;
Ok(())
}
fn compare(dense: &PersistentBitMatrix, sparse: &PersistentSparseBitMatrix) -> Result<(), Box<dyn Error>> {
let n = dense.n();
let n_cols = dense.n_cols();
assert_eq!(n, sparse.n(), "n mismatch");
assert_eq!(n_cols, sparse.n_cols(), "n_cols mismatch");
let mut dense_row = vec![0u32; n_cols];
let mut sparse_row = vec![0u32; n_cols];
let mut total_cells = 0usize;
let mut mismatched_cells = 0usize;
let mut first_mismatch = None;
let t0 = Instant::now();
for slot in 0..n {
dense.fill_row(slot, &mut dense_row);
sparse.fill_row(slot, &mut sparse_row);
for c in 0..n_cols {
total_cells += 1;
if dense_row[c] != sparse_row[c] {
mismatched_cells += 1;
if first_mismatch.is_none() {
first_mismatch = Some((slot, c, dense_row[c], sparse_row[c]));
}
}
}
}
println!("comparaison case-a-case: {:?} ({total_cells} cellules)", t0.elapsed());
if let Some((slot, col, d, s)) = first_mismatch {
eprintln!("PREMIER MISMATCH: slot={slot} col={col} dense={d} sparse={s}");
println!("MISMATCHES: {mismatched_cells} cellules differentes sur {total_cells}");
return Err("dense et sparse divergent".into());
}
println!("OK: toutes les {total_cells} cellules sont identiques entre dense et sparse.");
Ok(())
}
/// Content stats computed straight from the dense matrix — sparsity and the
/// singleton-row fraction the dictionary-dedup design targets (see
/// `compare_sparse_dense_count.rs`'s own `content_stats`).
fn content_stats(dense: &PersistentBitMatrix) {
let n = dense.n();
let n_cols = dense.n_cols();
let total_cells = n as u64 * n_cols as u64;
let mut set_cells = 0u64;
let mut singleton_rows = 0u64;
let mut row = vec![0u32; n_cols];
for slot in 0..n {
dense.fill_row(slot, &mut row);
let row_set = row.iter().filter(|&&v| v != 0).count() as u32;
set_cells += row_set as u64;
if row_set == 1 {
singleton_rows += 1;
}
}
println!("\n--- stats de contenu (source: dense) ---");
println!("cellules totales: {total_cells}");
println!(
"cellules non-nulles: {set_cells} ({:.3}% du total)",
100.0 * set_cells as f64 / total_cells as f64
);
println!(
"lignes singleton: {singleton_rows} / {n} ({:.3}%)",
100.0 * singleton_rows as f64 / n as f64
);
}
/// Directory size in bytes — sums every regular file, one level deep
/// (matches both the dense `presence/` layout and the sparse builder's
/// flat output directory).
fn dir_size(dir: &Path) -> std::io::Result<u64> {
let mut total = 0u64;
for entry in std::fs::read_dir(dir)? {
let entry = entry?;
if entry.file_type()?.is_file() {
total += entry.metadata()?.len();
}
}
Ok(total)
}
fn compaction_stats(dense_layer_dir: &Path, sparse_dir: &Path, dense: &PersistentBitMatrix) -> Result<(), Box<dyn Error>> {
let dense_size = dir_size(&dense_layer_dir.join("presence"))?;
let sparse_size = dir_size(sparse_dir)?;
let n = dense.n();
let n_cols = dense.n_cols();
let set_cells: u64 = {
let mut row = vec![0u32; n_cols];
let mut count = 0u64;
for slot in 0..n {
dense.fill_row(slot, &mut row);
count += row.iter().filter(|&&v| v != 0).count() as u64;
}
count
};
println!("\n--- compaction ---");
println!("taille dense (presence/): {dense_size} octets");
println!("taille sparse: {sparse_size} octets");
println!(
"ratio sparse/dense: {:.4} ({:.2}% de la taille dense)",
sparse_size as f64 / dense_size as f64,
100.0 * sparse_size as f64 / dense_size as f64
);
if set_cells > 0 {
println!(
"bits/valeur (dense, sur cellules non-nulles): {:.3}",
dense_size as f64 * 8.0 / set_cells as f64
);
println!(
"bits/valeur (sparse, sur cellules non-nulles): {:.3}",
sparse_size as f64 * 8.0 / set_cells as f64
);
}
let total_cells = n as u64 * n_cols as u64;
println!(
"bits/valeur (dense, sur toutes les cellules): {:.3}",
dense_size as f64 * 8.0 / total_cells as f64
);
println!(
"bits/valeur (sparse, sur toutes les cellules): {:.3}",
sparse_size as f64 * 8.0 / total_cells as f64
);
Ok(())
}
+1
View File
@@ -23,6 +23,7 @@ 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, pack_sparse_bit_matrix}; pub use sparse::{PersistentSparseBitMatrix, PersistentSparseBitMatrixBuilder, pack_sparse_bit_matrix};
pub(crate) use sparse::RowRank;
pub(crate) use pairwise::{fill_row_generic, par_col_reduce, pairwise2_matrix, pairwise_matrix, row_generic}; pub(crate) use pairwise::{fill_row_generic, par_col_reduce, pairwise2_matrix, pairwise_matrix, row_generic};
+22 -1
View File
@@ -111,6 +111,16 @@ fn read_varint(data: &[u8], pos: &mut usize) -> u32 {
// ── PersistentSparseBitMatrix ─────────────────────────────────────────────── // ── PersistentSparseBitMatrix ───────────────────────────────────────────────
/// A row's position within the singleton/multi split — the same
/// partition `for_each_genome_in_row` resolves internally. Exposed so
/// value-layer companions (e.g. `PersistentSparseCompactIntMatrix`) can
/// index their own per-row value streams with the identical rank,
/// without recomputing the split.
pub(crate) enum RowRank {
Singleton(usize),
Multi(usize),
}
/// Row-sparse `n × n_cols` binary matrix: each of the `n` rows (k-mer /// Row-sparse `n × n_cols` binary matrix: each of the `n` rows (k-mer
/// slots) stores only the small set of `n_cols` (genome) columns it has set /// slots) stores only the small set of `n_cols` (genome) columns it has set
/// — a singleton genome index, or a shared, deduplicated multi-genome set, /// — a singleton genome index, or a shared, deduplicated multi-genome set,
@@ -200,11 +210,22 @@ impl PersistentSparseBitMatrix {
out.into_boxed_slice() out.into_boxed_slice()
} }
/// This row's rank within whichever of the singleton/multi arrays holds
/// it — see [`RowRank`].
#[inline]
pub(crate) fn row_rank(&self, slot: usize) -> RowRank {
if self.is_multi.get(slot) {
RowRank::Multi(self.is_multi.rank1(slot) as usize)
} else {
RowRank::Singleton(self.is_multi.rank0(slot) as usize)
}
}
/// Calls `f` once per genome index present at `slot` — the shared /// Calls `f` once per genome index present at `slot` — the shared
/// decode branch (singleton vs. varint-encoded multi set) behind /// decode branch (singleton vs. varint-encoded multi set) behind
/// `fill_row`, `fill_row_bool`, and `fill_sub_matrix`. /// `fill_row`, `fill_row_bool`, and `fill_sub_matrix`.
#[inline] #[inline]
fn for_each_genome_in_row(&self, slot: usize, mut f: impl FnMut(usize)) { pub(crate) fn for_each_genome_in_row(&self, slot: usize, mut f: impl FnMut(usize)) {
if self.is_multi.get(slot) { if self.is_multi.get(slot) {
let pos = self.is_multi.rank1(slot) as usize; let pos = self.is_multi.rank1(slot) as usize;
let dict_id = self.multi.get(pos) as usize; let dict_id = self.multi.get(pos) as usize;
+1 -1
View File
@@ -16,7 +16,7 @@ use crate::tempbitvec::{TempBitVec, TempBitVecBuilder};
use crate::tempintvec::{TempCompactIntVec, TempCompactIntVecBuilder}; use crate::tempintvec::{TempCompactIntVec, TempCompactIntVecBuilder};
use crate::views::{nonzero_triples, IntSliceView}; use crate::views::{nonzero_triples, IntSliceView};
fn col_path(dir: &Path, col: usize) -> PathBuf { pub(crate) fn col_path(dir: &Path, col: usize) -> PathBuf {
dir.join(format!("col_{col:06}.pciv")) dir.join(format!("col_{col:06}.pciv"))
} }
+2
View File
@@ -11,6 +11,7 @@ mod layer_meta;
mod meta; mod meta;
mod mmap_file; mod mmap_file;
mod reader; mod reader;
mod sparse_intmatrix;
mod storage_kind; mod storage_kind;
mod tempbitvec; mod tempbitvec;
mod tempintvec; mod tempintvec;
@@ -27,6 +28,7 @@ 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};
pub use layer_meta::LayerMeta; pub use layer_meta::LayerMeta;
pub use reader::{PersistentCompactIntVec, Iter as CompactIntVecIter}; pub use reader::{PersistentCompactIntVec, Iter as CompactIntVecIter};
pub use sparse_intmatrix::{PersistentSparseCompactIntMatrix, PersistentSparseCompactIntMatrixBuilder, pack_sparse_compact_int_matrix};
pub use storage_kind::StorageKind; pub use storage_kind::StorageKind;
pub use tempbitvec::{TempBitVec, TempBitVecBuilder}; pub use tempbitvec::{TempBitVec, TempBitVecBuilder};
pub use tempintvec::{TempCompactIntVec, TempCompactIntVecBuilder}; pub use tempintvec::{TempCompactIntVec, TempCompactIntVecBuilder};
+337
View File
@@ -0,0 +1,337 @@
//! `PersistentSparseCompactIntMatrix` — row-major sparse integer matrix,
//! composed on top of [`PersistentSparseBitMatrix`] rather than
//! reimplementing it: the "which columns are non-zero" support is exactly
//! a sparse bit matrix (dictionary-deduplicated — presence patterns repeat
//! across rows, same measured benefit as the pure bit case), and only the
//! values layer is new.
//!
//! Unlike the support, values are **not** deduplicated: two rows can share
//! the same non-zero columns (same `dict_id`) while carrying different
//! counts, so a value run belongs to a row, never to a dictionary entry.
//! Two companion streams:
//! - `singleton_values.pciv`: one value per singleton row, dense over
//! `n_singleton`, indexed by the same `rank0(slot)` the support already
//! uses for its own `singleton` array (see
//! [`PersistentSparseBitMatrix::row_rank`]).
//! - `multi_values.pciv` + `multi_offsets` (`.efl`/`.efh`): one flat,
//! non-deduplicated value per `(multi row, occurrence)` pair, addressed
//! by a per-*row* Elias-Fano offset (indexed by `rank1(slot)`) —
//! deliberately distinct from the support's own `dict_offsets`, which is
//! keyed by `dict_id` and therefore locates a *shared* column range, not
//! a row's private value range.
//!
//! Both value streams reuse [`PersistentCompactIntVec`]'s 1-byte +
//! overflow encoding — tuned for exactly this crate's measured value
//! distribution (the large majority under 127, virtually all under 255).
use std::fs;
use std::io;
use std::path::{Path, PathBuf};
use ndarray::Array1;
use crate::bitmatrix::RowRank;
use crate::builder::PersistentCompactIntVecBuilder;
use crate::eliasfano::{EliasFano, EliasFanoBuilder};
use crate::reader::PersistentCompactIntVec;
use crate::traits::ColumnWeights;
use crate::{PersistentCompactIntMatrix, PersistentSparseBitMatrix, PersistentSparseBitMatrixBuilder};
fn singleton_values_path(dir: &Path) -> PathBuf { dir.join("singleton_values.pciv") }
fn multi_values_path(dir: &Path) -> PathBuf { dir.join("multi_values.pciv") }
fn multi_offsets_base(dir: &Path) -> PathBuf { dir.join("multi_offsets") }
/// Row batch size for [`PersistentSparseCompactIntMatrixBuilder::build_from_dense`].
const BUILD_FROM_DENSE_BATCH: usize = 4096;
// ── PersistentSparseCompactIntMatrix ────────────────────────────────────────
pub struct PersistentSparseCompactIntMatrix {
support: PersistentSparseBitMatrix,
singleton_values: PersistentCompactIntVec,
multi_values: PersistentCompactIntVec,
multi_offsets: EliasFano,
}
impl PersistentSparseCompactIntMatrix {
pub fn open(dir: &Path) -> io::Result<Self> {
Ok(Self {
support: PersistentSparseBitMatrix::open(dir)?,
singleton_values: PersistentCompactIntVec::open(&singleton_values_path(dir))?,
multi_values: PersistentCompactIntVec::open(&multi_values_path(dir))?,
multi_offsets: EliasFano::open(&multi_offsets_base(dir))?,
})
}
#[inline]
pub fn n(&self) -> usize { self.support.n() }
#[inline]
pub fn n_cols(&self) -> usize { self.support.n_cols() }
/// Calls `f(col, value)` once per non-zero column at `slot`, in
/// ascending column order — the shared decode branch behind every
/// accessor below. Zips the support's own column decode
/// ([`PersistentSparseBitMatrix::for_each_genome_in_row`]) with this
/// row's private slice of whichever value stream it belongs to.
#[inline]
fn for_each_cell_in_row(&self, slot: usize, mut f: impl FnMut(usize, u32)) {
match self.support.row_rank(slot) {
RowRank::Singleton(rank) => {
let v = self.singleton_values.get(rank);
self.support.for_each_genome_in_row(slot, |c| f(c, v));
}
RowRank::Multi(rank) => {
let start = self.multi_offsets.get(rank) as usize;
let mut i = 0usize;
self.support.for_each_genome_in_row(slot, |c| {
f(c, self.multi_values.get(start + i));
i += 1;
});
}
}
}
/// Column-major point lookup: value at column `c`, slot `slot` (0 if
/// absent).
pub fn get(&self, c: usize, slot: usize) -> u32 {
let mut found = 0u32;
self.for_each_cell_in_row(slot, |g, v| if g == c { found = v; });
found
}
pub fn row(&self, slot: usize) -> Box<[u32]> {
let mut out = vec![0u32; self.n_cols()];
self.fill_row(slot, &mut out);
out.into_boxed_slice()
}
/// Fill `buf[i]` with column `i`'s value at `slot` (0 if absent) —
/// mirrors [`PersistentSparseBitMatrix::fill_row`]'s signature.
pub fn fill_row(&self, slot: usize, buf: &mut [u32]) {
buf[..self.n_cols()].fill(0);
self.for_each_cell_in_row(slot, |g, v| buf[g] = v);
}
/// Like [`PersistentSparseBitMatrix::fill_sub_matrix`]: `out` has one
/// entry per column, each filled with that column's values at `slots`,
/// in `slots` order.
pub fn fill_sub_matrix(&self, slots: &[usize], out: &mut [Vec<u32>]) {
assert_eq!(out.len(), self.n_cols());
for col in out.iter_mut() {
col.clear();
col.resize(slots.len(), 0);
}
for (i, &slot) in slots.iter().enumerate() {
self.for_each_cell_in_row(slot, |g, v| out[g][i] = v);
}
}
/// Yields every `(idx into slots, col, value)` triple, row by row, in
/// `slots` order — mirrors
/// [`PersistentSparseBitMatrix::nonzero_iter`].
pub fn nonzero_iter<'a>(&'a self, slots: &'a [usize]) -> impl Iterator<Item = (usize, usize, u32)> + 'a {
let mut next_slot = 0usize;
let mut cur_idx = 0usize;
let mut buf: Vec<(usize, u32)> = Vec::new();
let mut buf_pos = 0usize;
std::iter::from_fn(move || loop {
if buf_pos < buf.len() {
let (g, v) = buf[buf_pos];
buf_pos += 1;
return Some((cur_idx, g, v));
}
if next_slot >= slots.len() {
return None;
}
cur_idx = next_slot;
let slot = slots[next_slot];
next_slot += 1;
buf.clear();
self.for_each_cell_in_row(slot, |g, v| buf.push((g, v)));
buf_pos = 0;
})
}
/// Per-column sum, via a naive row-first accumulation. Column-pairwise
/// distance matrices (`CountPartials`) are explicitly deferred here,
/// same as `BitPartials` was initially for
/// `PersistentSparseBitMatrix`: unlike the support, values aren't
/// deduplicated across rows, so the dict-driven co-occurrence shortcut
/// (`col_weights_and_pair_counts`) doesn't carry over — a row-major
/// pairwise rewrite for counts is future work, not part of this type
/// yet.
pub fn sum(&self) -> Array1<u64> {
let mut sums = vec![0u64; self.n_cols()];
for slot in 0..self.n() {
self.for_each_cell_in_row(slot, |g, v| sums[g] += v as u64);
}
Array1::from(sums)
}
pub fn count_nonzero(&self) -> Array1<u64> {
let mut counts = vec![0u64; self.n_cols()];
for slot in 0..self.n() {
self.for_each_cell_in_row(slot, |g, _| counts[g] += 1);
}
Array1::from(counts)
}
}
impl ColumnWeights for PersistentSparseCompactIntMatrix {
#[inline]
fn col_weights(&self) -> Array1<u64> { self.sum() }
#[inline]
fn partial_kmer_counts(&self) -> Array1<u64> { self.count_nonzero() }
}
// ── PersistentSparseCompactIntMatrixBuilder ─────────────────────────────────
pub struct PersistentSparseCompactIntMatrixBuilder {
support: PersistentSparseBitMatrixBuilder,
dir: PathBuf,
singleton_values: Vec<u32>,
multi_values: Vec<u32>,
multi_row_offsets: Vec<u64>,
}
impl PersistentSparseCompactIntMatrixBuilder {
pub fn new(n: usize, n_cols: usize, dir: &Path) -> io::Result<Self> {
fs::create_dir_all(dir)?;
Ok(Self {
support: PersistentSparseBitMatrixBuilder::new(n, n_cols, dir)?,
dir: dir.to_path_buf(),
singleton_values: Vec::new(),
multi_values: Vec::new(),
multi_row_offsets: Vec::new(),
})
}
/// Appends one row: `cols` sorted ascending, each `< n_cols`;
/// `values[i]` is the value at column `cols[i]`. Rows must be pushed
/// in slot order — same contract as
/// [`PersistentSparseBitMatrixBuilder::push_row`], which this
/// delegates the column side to.
///
/// Cardinality 0 or ≥2 both route through the multi/dict path (0 as a
/// genuine empty entry, same reasoning as the bit builder); only
/// cardinality exactly 1 takes the singleton shortcut.
pub fn push_row(&mut self, cols: &[u32], values: &[u32]) {
assert_eq!(cols.len(), values.len(), "cols/values length mismatch");
match cols.len() {
1 => self.singleton_values.push(values[0]),
_ => {
self.multi_row_offsets.push(self.multi_values.len() as u64);
self.multi_values.extend_from_slice(values);
}
}
self.support.push_row(cols);
}
/// Builds a sparse matrix from an already-built dense
/// [`PersistentCompactIntMatrix`] (`Columnar` or `Packed`) — same
/// batched-transpose shape as
/// [`PersistentSparseBitMatrixBuilder::build_from_dense`]: processes
/// rows in batches of [`BUILD_FROM_DENSE_BATCH`], driven by
/// `dense.nonzero_iter`, which walks each column once per batch
/// instead of one `get()` per cell in row-major order.
///
/// `nonzero_iter`'s triples arrive grouped by ascending column (outer
/// loop over columns) — so within a batch, each row's `(col, value)`
/// pairs are appended to its own buffer in strictly ascending column
/// order "for free", matching `push_row`'s ascending-`cols` contract
/// without an explicit sort.
pub fn build_from_dense(dense: &PersistentCompactIntMatrix, dir: &Path) -> io::Result<Self> {
let n = dense.n();
let n_cols = dense.n_cols();
let mut builder = Self::new(n, n_cols, dir)?;
let mut row_cols: Vec<Vec<u32>> = vec![Vec::new(); BUILD_FROM_DENSE_BATCH];
let mut row_vals: Vec<Vec<u32>> = vec![Vec::new(); BUILD_FROM_DENSE_BATCH];
let mut start = 0;
while start < n {
let end = (start + BUILD_FROM_DENSE_BATCH).min(n);
let slots: Vec<usize> = (start..end).collect();
for (row, c, v) in dense.nonzero_iter(&slots) {
row_cols[row].push(c as u32);
row_vals[row].push(v);
}
for i in 0..(end - start) {
builder.push_row(&row_cols[i], &row_vals[i]);
row_cols[i].clear();
row_vals[i].clear();
}
start = end;
}
Ok(builder)
}
pub fn close(self) -> io::Result<()> {
let Self { support, dir, singleton_values, multi_values, multi_row_offsets } = self;
support.close()?;
let mut sv = PersistentCompactIntVecBuilder::new(singleton_values.len(), &singleton_values_path(&dir))?;
for (i, &v) in singleton_values.iter().enumerate() {
sv.set(i, v);
}
sv.close()?;
let mut mv = PersistentCompactIntVecBuilder::new(multi_values.len(), &multi_values_path(&dir))?;
for (i, &v) in multi_values.iter().enumerate() {
mv.set(i, v);
}
mv.close()?;
let universe = multi_values.len() as u64 + 1;
let mut mo = EliasFanoBuilder::new(multi_row_offsets.len(), universe, &multi_offsets_base(&dir))?;
for &o in &multi_row_offsets {
mo.push(o);
}
mo.close()?;
Ok(())
}
pub fn finish(self) -> io::Result<PersistentSparseCompactIntMatrix> {
let dir = self.dir.clone();
self.close()?;
PersistentSparseCompactIntMatrix::open(&dir)
}
}
// ── pack_sparse_compact_int_matrix ──────────────────────────────────────────
/// Converts whichever dense count-matrix format is on disk at `dir`
/// (`Packed`, i.e. `matrix.pcmx`, or `Columnar`, i.e. `col_*.pciv` +
/// `meta.json`) into the sparse on-disk format in place — generic over the
/// dense source the same way [`super::bitmatrix::pack_sparse_bit_matrix`]
/// is: detects the format, transposes straight from it, and only removes
/// the old dense files once the new format is fully written (a crash
/// mid-conversion leaves the previous, still-valid format rather than a
/// half-written one). Idempotent — a no-op if `singleton_values.pciv`
/// already exists.
pub fn pack_sparse_compact_int_matrix(dir: &Path) -> io::Result<()> {
use crate::intmatrix::{col_path, ColumnarCompactIntMatrix, PackedCompactIntMatrix};
if singleton_values_path(dir).exists() {
return Ok(());
}
let packed_path = dir.join("matrix.pcmx");
if packed_path.exists() {
let dense = PersistentCompactIntMatrix::Packed(PackedCompactIntMatrix::open(&packed_path)?);
PersistentSparseCompactIntMatrixBuilder::build_from_dense(&dense, dir)?.close()?;
drop(dense);
fs::remove_file(&packed_path)?;
} else {
let dense = PersistentCompactIntMatrix::Columnar(ColumnarCompactIntMatrix::open(dir)?);
let n_cols = dense.n_cols();
PersistentSparseCompactIntMatrixBuilder::build_from_dense(&dense, dir)?.close()?;
drop(dense);
for c in 0..n_cols {
let _ = fs::remove_file(col_path(dir, c));
}
let _ = fs::remove_file(dir.join("meta.json"));
}
Ok(())
}
+1
View File
@@ -6,6 +6,7 @@ mod fixedintvec;
mod intmatrix; mod intmatrix;
mod rankselect; mod rankselect;
mod sparse; mod sparse;
mod sparse_intmatrix;
use tempfile::tempdir; use tempfile::tempdir;
@@ -0,0 +1,306 @@
use tempfile::tempdir;
use crate::{
pack_compact_int_matrix, pack_sparse_compact_int_matrix, ColumnWeights, PersistentCompactIntMatrix,
PersistentCompactIntMatrixBuilder, PersistentSparseCompactIntMatrix, PersistentSparseCompactIntMatrixBuilder,
};
/// Builds a dense `PersistentCompactIntMatrix` from column-major `u32` data
/// — mirrors `tests/intmatrix.rs`'s own `make_matrix` helper.
fn make_dense(cols: &[&[u32]]) -> (tempfile::TempDir, PersistentCompactIntMatrix) {
let n = cols.first().map_or(0, |c| c.len());
let dir = tempdir().unwrap();
let counts_dir = dir.path().join("counts");
let mut b = PersistentCompactIntMatrixBuilder::new(n, &counts_dir).unwrap();
for &col in cols {
let mut cb = b.add_col().unwrap();
for (slot, &v) in col.iter().enumerate() {
cb.set(slot, v);
}
cb.close().unwrap();
}
b.close().unwrap();
let m = PersistentCompactIntMatrix::open(dir.path()).unwrap();
(dir, m)
}
/// Builds a sparse int matrix directly from row-major `u32` data (one
/// slice per row, `n_cols` values each) — mirrors `tests/sparse.rs`'s
/// `make_sparse` helper.
fn make_sparse(rows: &[&[u32]], n_cols: usize) -> (tempfile::TempDir, PersistentSparseCompactIntMatrix) {
let dir = tempdir().unwrap();
let sparse_dir = dir.path().join("sparse");
let mut b = PersistentSparseCompactIntMatrixBuilder::new(rows.len(), n_cols, &sparse_dir).unwrap();
let mut cols = Vec::new();
let mut values = Vec::new();
for row in rows {
cols.clear();
values.clear();
for (c, &v) in row.iter().enumerate() {
if v != 0 {
cols.push(c as u32);
values.push(v);
}
}
b.push_row(&cols, &values);
}
let m = b.finish().unwrap();
(dir, m)
}
#[test]
fn basic_roundtrip_singletons_and_multi() {
// 5 rows, 4 genomes: mix of singleton, multi-value, and one row that
// shares its *support* (columns {0,1}) with another but has different
// values — the whole point of not deduplicating values.
let rows: Vec<&[u32]> = vec![
&[5, 0, 0, 0], // singleton: genome 0, value 5
&[0, 0, 3, 0], // singleton: genome 2, value 3
&[7, 2, 0, 0], // multi: {0:7, 1:2}
&[0, 0, 4, 9], // multi: {2:4, 3:9}
&[1, 6, 0, 0], // multi: same support {0,1} as row 2, different values
];
let (_dir, m) = make_sparse(&rows, 4);
assert_eq!(m.n(), 5);
assert_eq!(m.n_cols(), 4);
for (slot, &expected) in rows.iter().enumerate() {
assert_eq!(&*m.row(slot), expected, "row {slot}");
}
}
#[test]
fn get_matches_row() {
let rows: Vec<&[u32]> = vec![
&[3, 0, 5],
&[0, 8, 0],
&[1, 2, 4],
];
let (_dir, m) = make_sparse(&rows, 3);
for (slot, row) in rows.iter().enumerate() {
for (c, &expected) in row.iter().enumerate() {
assert_eq!(m.get(c, slot), expected, "slot {slot}, col {c}");
}
}
}
#[test]
fn fill_row_matches_row() {
let rows: Vec<&[u32]> = vec![
&[9, 0, 2],
&[0, 4, 0],
&[1, 1, 1],
];
let (_dir, m) = make_sparse(&rows, 3);
let mut buf = vec![0u32; 3];
for slot in 0..3 {
m.fill_row(slot, &mut buf);
assert_eq!(&buf, rows[slot], "slot {slot}");
}
}
#[test]
fn fill_sub_matrix_matches_row() {
let rows: Vec<&[u32]> = vec![
&[9, 0, 2, 0],
&[0, 4, 0, 7],
&[1, 1, 1, 1],
&[0, 0, 0, 3],
&[5, 0, 0, 0],
];
let (_dir, m) = make_sparse(&rows, 4);
let slots = [4usize, 0, 2];
let mut sub: Vec<Vec<u32>> = vec![Vec::new(); 4];
m.fill_sub_matrix(&slots, &mut sub);
for (c, col) in sub.iter().enumerate() {
for (i, &slot) in slots.iter().enumerate() {
assert_eq!(col[i], rows[slot][c], "col {c}, slot {slot}");
}
}
}
#[test]
fn nonzero_iter_matches_row() {
let rows: Vec<&[u32]> = vec![
&[9, 0, 2],
&[0, 4, 0],
&[1, 1, 1],
];
let (_dir, m) = make_sparse(&rows, 3);
let slots = [2usize, 0, 1];
let mut expected: Vec<(usize, usize, u32)> = Vec::new();
for (idx, &slot) in slots.iter().enumerate() {
for (c, &v) in rows[slot].iter().enumerate() {
if v != 0 {
expected.push((idx, c, v));
}
}
}
let mut got: Vec<(usize, usize, u32)> = m.nonzero_iter(&slots).collect();
got.sort();
expected.sort();
assert_eq!(got, expected);
}
#[test]
fn sum_and_count_nonzero_match_naive() {
let rows: Vec<&[u32]> = vec![
&[9, 0, 2],
&[0, 4, 0],
&[1, 1, 1],
&[0, 0, 6],
];
let (_dir, m) = make_sparse(&rows, 3);
let mut expected_sum = [0u64; 3];
let mut expected_count = [0u64; 3];
for row in &rows {
for (c, &v) in row.iter().enumerate() {
expected_sum[c] += v as u64;
if v != 0 {
expected_count[c] += 1;
}
}
}
assert_eq!(m.sum().to_vec(), expected_sum.to_vec());
assert_eq!(m.count_nonzero().to_vec(), expected_count.to_vec());
assert_eq!(m.col_weights().to_vec(), expected_sum.to_vec());
assert_eq!(m.partial_kmer_counts().to_vec(), expected_count.to_vec());
}
#[test]
fn single_genome_matrix() {
// n_cols=1 — every row is necessarily a singleton or empty.
let rows: Vec<&[u32]> = vec![&[5], &[0], &[3]];
let (_dir, m) = make_sparse(&rows, 1);
for (slot, &expected) in rows.iter().enumerate() {
assert_eq!(&*m.row(slot), expected, "row {slot}");
}
}
#[test]
fn empty_row_roundtrips() {
// Cardinality 0 (no genome present) routes through the multi/dict
// path as a genuine empty entry — must not be confused with a
// singleton or crash.
let rows: Vec<&[u32]> = vec![&[0, 0, 0], &[7, 0, 0], &[0, 0, 0]];
let (_dir, m) = make_sparse(&rows, 3);
for (slot, &expected) in rows.iter().enumerate() {
assert_eq!(&*m.row(slot), expected, "row {slot}");
}
}
#[test]
fn reopen_after_close_matches_original() {
let rows: Vec<Vec<u32>> = (0..500)
.map(|i| {
(0..37)
.map(|c| if (i * 7 + c * 3) % 11 == 0 { ((i + c) % 250 + 1) as u32 } else { 0 })
.collect()
})
.collect();
let dir = tempdir().unwrap();
let sparse_dir = dir.path().join("sparse");
{
let mut b = PersistentSparseCompactIntMatrixBuilder::new(rows.len(), 37, &sparse_dir).unwrap();
let mut cols = Vec::new();
let mut values = Vec::new();
for row in &rows {
cols.clear();
values.clear();
for (c, &v) in row.iter().enumerate() {
if v != 0 {
cols.push(c as u32);
values.push(v);
}
}
b.push_row(&cols, &values);
}
b.close().unwrap();
} // builder + every mmap fully dropped here
let m = PersistentSparseCompactIntMatrix::open(&sparse_dir).unwrap();
assert_eq!(m.n(), 500);
assert_eq!(m.n_cols(), 37);
for (slot, expected) in rows.iter().enumerate() {
assert_eq!(&*m.row(slot), expected.as_slice(), "row {slot}");
}
}
#[test]
fn overflow_values_roundtrip() {
// Values >= 255 exercise the PersistentCompactIntVec overflow path in
// both the singleton and multi value streams.
let rows: Vec<&[u32]> = vec![
&[1000, 0, 0],
&[0, 50_000, 300],
];
let (_dir, m) = make_sparse(&rows, 3);
for (slot, &expected) in rows.iter().enumerate() {
assert_eq!(&*m.row(slot), expected, "row {slot}");
}
}
/// Common case: `pack_sparse_compact_int_matrix` runs right after a build,
/// while `counts/` is still `Columnar` (`matrix.pcmx` doesn't exist yet).
/// Must transpose straight from the per-column `.pciv` files — never write
/// a full packed copy just to read it back and delete it. Mirrors
/// `tests::sparse::pack_sparse_bit_matrix_transposes_columnar_directly`.
#[test]
fn pack_sparse_compact_int_matrix_transposes_columnar_directly() {
let cols: &[&[u32]] = &[
&[9, 0, 2, 0],
&[0, 4, 0, 7],
&[1, 1, 300, 1],
];
let (dir, dense) = make_dense(cols);
let counts_dir = dir.path().join("counts");
assert!(!counts_dir.join("matrix.pcmx").exists(), "still Columnar before packing");
let expected_rows: Vec<Box<[u32]>> = (0..dense.n()).map(|s| dense.row(s)).collect();
let n_cols = dense.n_cols();
drop(dense);
pack_sparse_compact_int_matrix(&counts_dir).unwrap();
assert!(counts_dir.join("singleton_values.pciv").exists());
assert!(!counts_dir.join("meta.json").exists(), "columnar meta.json cleaned up");
assert!(!counts_dir.join("col_000000.pciv").exists(), "columnar column files cleaned up");
assert!(!counts_dir.join("matrix.pcmx").exists(), "never materialised");
let sparse = PersistentSparseCompactIntMatrix::open(&counts_dir).unwrap();
assert_eq!(sparse.n(), expected_rows.len());
assert_eq!(sparse.n_cols(), n_cols);
for (slot, expected) in expected_rows.iter().enumerate() {
assert_eq!(&*sparse.row(slot), &**expected, "slot {slot}");
}
// Idempotent — a second call is a no-op, doesn't error.
pack_sparse_compact_int_matrix(&counts_dir).unwrap();
}
/// Regression for the pre-existing path: if `matrix.pcmx` is already there
/// (e.g. `pack_compact_int_matrix` ran earlier for a non-sparse consumer),
/// it's still used as the transpose source and cleaned up afterward.
#[test]
fn pack_sparse_compact_int_matrix_from_already_packed_matrix() {
let cols: &[&[u32]] = &[
&[9, 0, 2],
&[0, 4, 300],
];
let (dir, dense) = make_dense(cols);
let counts_dir = dir.path().join("counts");
let expected_rows: Vec<Box<[u32]>> = (0..dense.n()).map(|s| dense.row(s)).collect();
let n_cols = dense.n_cols();
drop(dense);
pack_compact_int_matrix(&counts_dir).unwrap();
assert!(counts_dir.join("matrix.pcmx").exists());
pack_sparse_compact_int_matrix(&counts_dir).unwrap();
assert!(!counts_dir.join("matrix.pcmx").exists(), "packed intermediate cleaned up");
let sparse = PersistentSparseCompactIntMatrix::open(&counts_dir).unwrap();
assert_eq!(sparse.n_cols(), n_cols);
for (slot, expected) in expected_rows.iter().enumerate() {
assert_eq!(&*sparse.row(slot), &**expected, "slot {slot}");
}
}