diff --git a/src/obicompactvec/examples/compare_sparse_dense_count.rs b/src/obicompactvec/examples/compare_sparse_dense_count.rs new file mode 100644 index 00000000..36040be7 --- /dev/null +++ b/src/obicompactvec/examples/compare_sparse_dense_count.rs @@ -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 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> { + let layer_dir = std::env::args().nth(1).expect("usage: compare_sparse_dense_count "); + 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> { + 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 { + 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> { + 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(()) +} diff --git a/src/obicompactvec/examples/compare_sparse_dense_presence.rs b/src/obicompactvec/examples/compare_sparse_dense_presence.rs new file mode 100644 index 00000000..04e41c00 --- /dev/null +++ b/src/obicompactvec/examples/compare_sparse_dense_presence.rs @@ -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 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> { + let layer_dir = std::env::args().nth(1).expect("usage: compare_sparse_dense_presence "); + 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> { + 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 { + 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> { + 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(()) +} diff --git a/src/obicompactvec/src/bitmatrix/mod.rs b/src/obicompactvec/src/bitmatrix/mod.rs index 40bb06e9..4aefa62b 100644 --- a/src/obicompactvec/src/bitmatrix/mod.rs +++ b/src/obicompactvec/src/bitmatrix/mod.rs @@ -23,6 +23,7 @@ pub use builder::PersistentBitMatrixBuilder; pub use packed::pack_bit_matrix; pub use persistent::PersistentBitMatrix; 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}; diff --git a/src/obicompactvec/src/bitmatrix/sparse.rs b/src/obicompactvec/src/bitmatrix/sparse.rs index 156071e5..b8967edf 100644 --- a/src/obicompactvec/src/bitmatrix/sparse.rs +++ b/src/obicompactvec/src/bitmatrix/sparse.rs @@ -111,6 +111,16 @@ fn read_varint(data: &[u8], pos: &mut usize) -> u32 { // ── 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 /// 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, @@ -200,11 +210,22 @@ impl PersistentSparseBitMatrix { 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 /// decode branch (singleton vs. varint-encoded multi set) behind /// `fill_row`, `fill_row_bool`, and `fill_sub_matrix`. #[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) { let pos = self.is_multi.rank1(slot) as usize; let dict_id = self.multi.get(pos) as usize; diff --git a/src/obicompactvec/src/intmatrix.rs b/src/obicompactvec/src/intmatrix.rs index 438b328f..9d4776c1 100644 --- a/src/obicompactvec/src/intmatrix.rs +++ b/src/obicompactvec/src/intmatrix.rs @@ -16,7 +16,7 @@ use crate::tempbitvec::{TempBitVec, TempBitVecBuilder}; use crate::tempintvec::{TempCompactIntVec, TempCompactIntVecBuilder}; 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")) } diff --git a/src/obicompactvec/src/lib.rs b/src/obicompactvec/src/lib.rs index 84cf95e8..e1a15a38 100644 --- a/src/obicompactvec/src/lib.rs +++ b/src/obicompactvec/src/lib.rs @@ -11,6 +11,7 @@ mod layer_meta; mod meta; mod mmap_file; mod reader; +mod sparse_intmatrix; mod storage_kind; mod tempbitvec; 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 layer_meta::LayerMeta; 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 tempbitvec::{TempBitVec, TempBitVecBuilder}; pub use tempintvec::{TempCompactIntVec, TempCompactIntVecBuilder}; diff --git a/src/obicompactvec/src/sparse_intmatrix.rs b/src/obicompactvec/src/sparse_intmatrix.rs new file mode 100644 index 00000000..0cd3d2c9 --- /dev/null +++ b/src/obicompactvec/src/sparse_intmatrix.rs @@ -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 { + 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]) { + 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 + '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 { + 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 { + 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 { self.sum() } + #[inline] + fn partial_kmer_counts(&self) -> Array1 { self.count_nonzero() } +} + +// ── PersistentSparseCompactIntMatrixBuilder ───────────────────────────────── + +pub struct PersistentSparseCompactIntMatrixBuilder { + support: PersistentSparseBitMatrixBuilder, + dir: PathBuf, + singleton_values: Vec, + multi_values: Vec, + multi_row_offsets: Vec, +} + +impl PersistentSparseCompactIntMatrixBuilder { + pub fn new(n: usize, n_cols: usize, dir: &Path) -> io::Result { + 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 { + 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![Vec::new(); BUILD_FROM_DENSE_BATCH]; + let mut row_vals: Vec> = 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 = (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 { + 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(()) +} diff --git a/src/obicompactvec/src/tests/mod.rs b/src/obicompactvec/src/tests/mod.rs index ee22d26e..4f6fd9ac 100644 --- a/src/obicompactvec/src/tests/mod.rs +++ b/src/obicompactvec/src/tests/mod.rs @@ -6,6 +6,7 @@ mod fixedintvec; mod intmatrix; mod rankselect; mod sparse; +mod sparse_intmatrix; use tempfile::tempdir; diff --git a/src/obicompactvec/src/tests/sparse_intmatrix.rs b/src/obicompactvec/src/tests/sparse_intmatrix.rs new file mode 100644 index 00000000..c4089da8 --- /dev/null +++ b/src/obicompactvec/src/tests/sparse_intmatrix.rs @@ -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![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> = (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> = (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> = (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}"); + } +}