OBIKmers/obikmer
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 27 Wiki Activity

feat(obicompactvec): introduce sparse bit matrix with supporting primitives

Browse Source 
Implements a compact, row-major sparse bit matrix backed by memory-mapped components, introducing EliasFano, PersistentFixedIntVec, and PersistentRankSelectBitVec primitives for efficient storage and decoding. Adds a BinaryMatrix trait to unify row-level operations across dense and sparse implementations. Corrects edge-case behaviors for zero-width bit storage and cardinality-0 rows. Delivers reduced on-disk size and faster random row access, with column reads remaining dense-only. Test suites and benchmarks are included but currently marked as ignored.
This commit is contained in:
Eric Coissac
2026-08-16 21:43:09 +02:00
parent 45b19503a1
commit 50f4820cb9
17 changed files with 2024 additions and 2 deletions
Download Patch File Download Diff File Expand all files Collapse all files
docmd/architecture/siblings.md
+93
View File
@@ -403,3 +403,96 @@ the total base-occurrence count, not the genome count (the two coincide
only when no genome carries more than one base). Kept side by side
specifically to measure, on real data, how much the two diverge — not yet
analyzed.
## `PersistentSparseBitMatrix` — implemented and measured (2026-08-15)
A row-major (k-mer-major), deduplicated sparse alternative to
`obicompactvec::PersistentBitMatrix`, motivated by the same sparsity that
drove `--subsample`/`--shannon` above, but pursued as a foundational
storage-layer change rather than an index-level workaround. Full design
history, rationale, and rejected alternatives (external Elias-Fano crates,
`cacheline-ef`, a single unsplit `dict_id` array) are in the dedicated
implementation plan (`vivid-mapping-tiger.md` at the time of writing — the
content below is the durable summary, not a pointer to a session-scoped
file). Also directly informed by Alanko, Bille, Gørtz, Navarro, Puglisi,
"Compact Data Structures for Collections of Sets" (2025,
`biblio/Alanko et al. - Compact Data Structures for Collections of
Sets.pdf`) — this design implements only their exact-duplicate special
case (a plain dedup dictionary), not their full subset-containment
hierarchy.
**Design**: four on-disk components, each mmap-backed, built once per
layer (matching how the rest of the build pipeline already works — never
the whole multi-billion-row index at once): an `is_multi` rank-capable
flag per row (singleton vs. multi-genome), a fixed-bit-width array for
singleton rows (genome index directly, `ceil(log2(n_cols))` bits), a
separate fixed-bit-width array for multi-genome rows (`dict_id`,
`ceil(log2(n_distinct_multi_sets))` bits — kept apart from the singleton
array specifically because `n_distinct_multi_sets` can be large in
absolute terms even when multi-genome rows are a small *fraction* of all
rows, and a single shared array would force every row, singletons
included, to pay the wider width), and a deduplicated dictionary of
distinct multi-genome sets (Elias-Fano-encoded byte offsets + a
varint-encoded values blob). New low-level primitives added to
`obicompactvec` to build this: `PersistentFixedIntVec` (arbitrary,
runtime-parameterized bit width, width 0 included — needed once a real
bug surfaced, see below), `PersistentRankSelectBitVec` (rank1/rank0/select1
on top of the crate's existing `count_ones`, using
`common_traits::SelectInWord`), `EliasFano` (composes the two). A new
`BinaryMatrix` trait (`n`, `n_cols`, `row`/`fill_row`, `fill_sub_matrix`,
`count_ones`) unifies dense and sparse at the one call site that needs
both interchangeably (`obikphylo::siblings::cache::Mat`) — column-oriented
methods (`col`, `col_view`, the `partial_*_dist_matrix` family) stay
dense-only.
**Two real bugs caught by tests, not by inspection**: (1) `EliasFano::open`
re-derived its low-bits width from the persisted low-vector file's own
width byte; the zero-width case was built with a dummy 1-bit placeholder
(the builder rejected true width 0), so every reopened value silently
doubled. Fixed by making `PersistentFixedIntVec` genuinely support width 0
(no storage, `get` always 0) instead of working around the limitation in
`EliasFano`. (2) An empty row (cardinality 0 — not expected on a real
built index, but not guarded against either) was recorded as a singleton
at genome 0, indistinguishable on read-back from a *real* singleton at
genome 0. Fixed by routing cardinality-0 rows through the dictionary path
(a genuine empty entry) instead of the singleton shortcut. Both caught by
`obicompactvec`'s test suite (142 tests, including disk-reopen round-trips
that drop every builder/mmap before reopening fresh), not by manual
review — worth remembering next time a "this edge case can't happen in
practice" shortcut is tempting.
**Measured on real data** (`layer_1` of `phyloskims_sal_vac`'s
`part_00018`, 30,246,774 rows, 91 genomes — `#[ignore]`d benchmarks in
`obikphylo/src/siblings/tests.rs`):
| | dense | sparse | ratio |
|---|---|---|---|
| on-disk size | 328.1MB | 43.7MB | **7.5x** smaller |
| build time / peak RSS | — | 4.26s / 628MB | (per-layer, in-memory construction — comfortable) |
| row access, sequential (2M reads) | 43ns/row | 32ns/row | sparse **faster** (smaller structure, better cache fit) |
| row access, random (2M reads) | 409ns/row | 85ns/row | sparse **~4.8x faster** (the real `--shannon`/family-lookup shape) |
| column access, one full column (30.2M rows) | 11.5ms | 993ms | sparse **86x slower** (no native column method — every read decodes a full row to keep one bit) |
The row-access wins (both directions) weren't the design's stated goal —
compactness was — but turn out real: dense's genome-major layout scatters
a single row read across a much bigger file, which costs more than
sparse's rank/select/varint decode once the file is this much smaller.
The column-access cost is the flip side of the same layout choice, and is
exactly what the next item below exists to fix.
**Next, not yet planned**: rewrite `partial_jaccard_dist_matrix`/
`partial_hamming_dist_matrix`/etc. (`obicompactvec/src/bitmatrix/pairwise.rs`)
as a row-major co-occurrence accumulation (`O(Σ_rows k²)`, per-row
increments into an `NxN` genome-pair counter — the known alternative to
today's column-fold, plausibly cheaper on data this sparse, not just a
fallback) so `obikindex`'s `--metric`/distance-matrix path can use the
sparse type without the measured 86x column-access penalty. Needs its own
design pass (in particular how it plugs into the `BitPartials`/
`ColumnWeights` traits so both matrix types keep serving `--metric`)
before implementation — not just "port the loop", a genuinely different
algorithm.
Also still deferred, unchanged from the implementation plan: full Alanko
et al. subset-hierarchy compression (only the exact-duplicate special case
is built), a sparse `PersistentCompactIntMatrix` (count matrices), and
BRWT-style column-correlation exploitation.
src/Cargo.lock
Generated
+1
View File
@@ -1721,6 +1721,7 @@ checksum = "830b246a0e5f20af87141b25c173cd1b609bd7779a4617d6ec582abaf90870f3"
name = "obicompactvec"
version = "0.1.0"
dependencies = [
"common_traits",
"memmap2",
"ndarray",
"rayon",
src/obicompactvec/Cargo.toml
+1
View File
@@ -4,6 +4,7 @@ version = "0.1.0"
edition = "2024"
[dependencies]
common_traits = "0.11"
memmap2 = "0.9"
ndarray = "0.16"
rayon = "1"
src/obicompactvec/src/bitmatrix/mod.rs
+2
View File
@@ -17,10 +17,12 @@ mod group_ops;
mod packed;
mod pairwise;
mod persistent;
mod sparse;
pub use builder::PersistentBitMatrixBuilder;
pub use packed::pack_bit_matrix;
pub use persistent::PersistentBitMatrix;
pub use sparse::{PersistentSparseBitMatrix, PersistentSparseBitMatrixBuilder};
pub(crate) use pairwise::{pairwise_matrix, pairwise2_matrix};
src/obicompactvec/src/bitmatrix/persistent.rs
+15
View File
@@ -241,3 +241,18 @@ impl BitPartials for PersistentBitMatrix {
self.partial_hamming_dist_matrix()
}
}
impl crate::traits::BinaryMatrix for PersistentBitMatrix {
#[inline]
fn n(&self) -> usize { self.n() }
#[inline]
fn n_cols(&self) -> usize { self.n_cols() }
#[inline]
fn row(&self, slot: usize) -> Box<[bool]> { self.row(slot) }
#[inline]
fn fill_row(&self, slot: usize, buf: &mut [u32]) { self.fill_row(slot, buf) }
#[inline]
fn fill_sub_matrix(&self, slots: &[usize], out: &mut [Vec<bool>]) { self.fill_sub_matrix(slots, out) }
#[inline]
fn count_ones(&self) -> Array1<u64> { self.count_ones() }
}
src/obicompactvec/src/bitmatrix/sparse.rs
+372
View File
@@ -0,0 +1,372 @@
//! `PersistentSparseBitMatrix` — row-major (k-mer-major), deduplicated
//! sparse alternative to [`super::PersistentBitMatrix`]. See
//! `docmd/architecture/siblings.md` and the sparse-matrix design plan for
//! the full rationale (measured sparsity/duplication on real data). Not
//! used by any production code path yet — a new type, not a replacement.
//!
//! On-disk layout (a directory, mirroring [`super::PersistentBitMatrix`]'s
//! own directory-of-files convention): `meta.json` plus four components —
//! `is_multi.prsb` (rank-capable flag: singleton row vs. multi-genome
//! row), `singleton.pfiv` (genome index, one entry per singleton row,
//! `ceil(log2(n_cols))` bits each), `multi.pfiv` (`dict_id`, one entry per
//! multi-genome row, `ceil(log2(n_distinct_multi))` bits each),
//! `dict_offsets` (Elias-Fano — `.efl`/`.efh` — one entry per *distinct*
//! multi-genome set, byte offset into `dict_values.bin`), `dict_values.bin`
//! (varint-encoded sorted genome-index list per distinct set).
//!
//! A row's genome set is read by: check `is_multi[slot]`; if singleton,
//! `singleton[rank0(is_multi, slot)]` is the genome index directly; if
//! multi, `multi[rank1(is_multi, slot)]` is a `dict_id`, and
//! `dict_values[dict_offsets[dict_id]..dict_offsets[dict_id+1]]` is its
//! varint-encoded genome list.
use std::collections::HashMap;
use std::fs;
use std::io;
use std::path::{Path, PathBuf};
use ndarray::Array1;
use crate::eliasfano::{EliasFano, EliasFanoBuilder};
use crate::fixedintvec::{PersistentFixedIntVec, PersistentFixedIntVecBuilder, bit_width_for_range};
use crate::meta::field;
use crate::rankselect::{PersistentRankSelectBitVec, PersistentRankSelectBitVecBuilder};
use super::PersistentBitMatrix;
fn is_multi_path(dir: &Path) -> PathBuf { dir.join("is_multi.prsb") }
fn singleton_path(dir: &Path) -> PathBuf { dir.join("singleton.pfiv") }
fn multi_path(dir: &Path) -> PathBuf { dir.join("multi.pfiv") }
fn dict_offsets_base(dir: &Path) -> PathBuf { dir.join("dict_offsets") }
fn dict_values_path(dir: &Path) -> PathBuf { dir.join("dict_values.bin") }
fn meta_path(dir: &Path) -> PathBuf { dir.join("meta.json") }
struct SparseMeta {
n: usize,
n_cols: usize,
n_singleton: usize,
n_multi: usize,
n_distinct_multi: usize,
}
impl SparseMeta {
fn load(dir: &Path) -> io::Result<Self> {
let s = fs::read_to_string(meta_path(dir))?;
let get = |name: &str| {
field(&s, name).ok_or_else(|| io::Error::new(io::ErrorKind::InvalidData, format!("bad meta.json: missing {name}")))
};
Ok(Self {
n: get("n")?,
n_cols: get("n_cols")?,
n_singleton: get("n_singleton")?,
n_multi: get("n_multi")?,
n_distinct_multi: get("n_distinct_multi")?,
})
}
fn save(&self, dir: &Path) -> io::Result<()> {
fs::write(
meta_path(dir),
format!(
"{{\"n\":{},\"n_cols\":{},\"n_singleton\":{},\"n_multi\":{},\"n_distinct_multi\":{}}}\n",
self.n, self.n_cols, self.n_singleton, self.n_multi, self.n_distinct_multi,
),
)
}
}
// ── varint (LEB128-style) ────────────────────────────────────────────────────
fn write_varint(buf: &mut Vec<u8>, mut v: u32) {
loop {
let byte = (v & 0x7F) as u8;
v >>= 7;
if v == 0 {
buf.push(byte);
break;
}
buf.push(byte | 0x80);
}
}
fn read_varint(data: &[u8], pos: &mut usize) -> u32 {
let mut result = 0u32;
let mut shift = 0u32;
loop {
let byte = data[*pos];
*pos += 1;
result |= ((byte & 0x7F) as u32) << shift;
if byte & 0x80 == 0 {
break;
}
shift += 7;
}
result
}
// ── PersistentSparseBitMatrix ───────────────────────────────────────────────
pub struct PersistentSparseBitMatrix {
is_multi: PersistentRankSelectBitVec,
singleton: PersistentFixedIntVec,
multi: PersistentFixedIntVec,
dict_offsets: EliasFano,
dict_values: Vec<u8>,
n: usize,
n_cols: usize,
n_distinct_multi: usize,
}
impl PersistentSparseBitMatrix {
pub fn open(dir: &Path) -> io::Result<Self> {
let meta = SparseMeta::load(dir)?;
let is_multi = PersistentRankSelectBitVec::open(&is_multi_path(dir))?;
let singleton = PersistentFixedIntVec::open(&singleton_path(dir))?;
let multi = PersistentFixedIntVec::open(&multi_path(dir))?;
let dict_offsets = EliasFano::open(&dict_offsets_base(dir))?;
let dict_values = fs::read(dict_values_path(dir))?;
Ok(Self {
is_multi, singleton, multi, dict_offsets, dict_values,
n: meta.n, n_cols: meta.n_cols, n_distinct_multi: meta.n_distinct_multi,
})
}
#[inline]
pub fn n(&self) -> usize { self.n }
#[inline]
pub fn n_cols(&self) -> usize { self.n_cols }
/// Byte range of `dict_id`'s varint-encoded genome list within
/// `dict_values`.
#[inline]
fn dict_entry_range(&self, dict_id: usize) -> (usize, usize) {
let start = self.dict_offsets.get(dict_id) as usize;
let end = if dict_id + 1 < self.n_distinct_multi {
self.dict_offsets.get(dict_id + 1) as usize
} else {
self.dict_values.len()
};
(start, end)
}
pub fn row(&self, slot: usize) -> Box<[bool]> {
let mut out = vec![false; self.n_cols];
self.fill_row_bool(slot, &mut out);
out.into_boxed_slice()
}
/// Fill `buf[i]` with `1` iff genome `i` is present at `slot`, else `0`
/// — mirrors [`super::PersistentBitMatrix::fill_row`]'s signature.
pub fn fill_row(&self, slot: usize, buf: &mut [u32]) {
buf[..self.n_cols].fill(0);
if self.is_multi.get(slot) {
let pos = self.is_multi.rank1(slot) as usize;
let dict_id = self.multi.get(pos) as usize;
let (start, end) = self.dict_entry_range(dict_id);
let mut p = start;
while p < end {
buf[read_varint(&self.dict_values, &mut p) as usize] = 1;
}
} else {
let pos = self.is_multi.rank0(slot) as usize;
buf[self.singleton.get(pos) as usize] = 1;
}
}
fn fill_row_bool(&self, slot: usize, buf: &mut [bool]) {
buf.fill(false);
if self.is_multi.get(slot) {
let pos = self.is_multi.rank1(slot) as usize;
let dict_id = self.multi.get(pos) as usize;
let (start, end) = self.dict_entry_range(dict_id);
let mut p = start;
while p < end {
buf[read_varint(&self.dict_values, &mut p) as usize] = true;
}
} else {
let pos = self.is_multi.rank0(slot) as usize;
buf[self.singleton.get(pos) as usize] = true;
}
}
/// Column-oriented per-genome k-mer totals — a naive row-by-row scan,
/// not the O(1)-per-column reduction the dense matrix's `count_ones`
/// is. Deliberately not optimised: see `docmd/architecture/siblings.md`
/// and the sparse-matrix plan's "Explicitly deferred" — column-side
/// access stays correct but slow on this type for now.
pub fn count_ones(&self) -> Array1<u64> {
let mut counts = vec![0u64; self.n_cols];
let mut buf = vec![false; self.n_cols];
for slot in 0..self.n {
self.fill_row_bool(slot, &mut buf);
for (c, &present) in buf.iter().enumerate() {
if present {
counts[c] += 1;
}
}
}
Array1::from(counts)
}
/// Like [`super::PersistentBitMatrix::fill_sub_matrix`]: `out` has one
/// entry per genome column, each filled with that column's values at
/// `slots`, in `slots` order. Naive (row-major decode, scattered into
/// column buffers) — see [`count_ones`](Self::count_ones)'s doc comment.
pub fn fill_sub_matrix(&self, slots: &[usize], out: &mut [Vec<bool>]) {
assert_eq!(out.len(), self.n_cols);
for col in out.iter_mut() {
col.clear();
col.resize(slots.len(), false);
}
let mut buf = vec![false; self.n_cols];
for (i, &slot) in slots.iter().enumerate() {
self.fill_row_bool(slot, &mut buf);
for (c, &present) in buf.iter().enumerate() {
out[c][i] = present;
}
}
}
}
impl crate::traits::BinaryMatrix for PersistentSparseBitMatrix {
#[inline]
fn n(&self) -> usize { self.n() }
#[inline]
fn n_cols(&self) -> usize { self.n_cols() }
#[inline]
fn row(&self, slot: usize) -> Box<[bool]> { self.row(slot) }
#[inline]
fn fill_row(&self, slot: usize, buf: &mut [u32]) { self.fill_row(slot, buf) }
#[inline]
fn fill_sub_matrix(&self, slots: &[usize], out: &mut [Vec<bool>]) { self.fill_sub_matrix(slots, out) }
#[inline]
fn count_ones(&self) -> Array1<u64> { self.count_ones() }
}
// ── PersistentSparseBitMatrixBuilder ────────────────────────────────────────
pub struct PersistentSparseBitMatrixBuilder {
dir: PathBuf,
n_cols: usize,
is_multi: Vec<bool>,
singleton_values: Vec<u32>,
multi_dict_ids: Vec<u32>,
dedup: HashMap<Vec<u32>, u32>,
dict_sets_in_order: Vec<Vec<u32>>,
}
impl PersistentSparseBitMatrixBuilder {
pub fn new(n: usize, n_cols: usize, dir: &Path) -> io::Result<Self> {
fs::create_dir_all(dir)?;
Ok(Self {
dir: dir.to_path_buf(),
n_cols,
is_multi: Vec::with_capacity(n),
singleton_values: Vec::new(),
multi_dict_ids: Vec::new(),
dedup: HashMap::new(),
dict_sets_in_order: Vec::new(),
})
}
/// Appends one row's genome set — `genomes` sorted ascending, each
/// `< n_cols`. Rows must be pushed in slot order (0, 1, 2, ...), one
/// call per k-mer slot of the layer being converted/built.
///
/// Cardinality 0 (no genome present — not expected on a real built
/// index, every slot has at least one, but handled correctly rather
/// than assumed away) routes through the dictionary path as a genuine
/// empty entry, *not* the singleton shortcut: a placeholder singleton
/// value would be indistinguishable from a real singleton at that same
/// genome index on read-back.
pub fn push_row(&mut self, genomes: &[u32]) {
match genomes.len() {
1 => {
self.is_multi.push(false);
self.singleton_values.push(genomes[0]);
}
_ => {
self.is_multi.push(true);
let id = if let Some(&id) = self.dedup.get(genomes) {
id
} else {
let id = self.dict_sets_in_order.len() as u32;
self.dict_sets_in_order.push(genomes.to_vec());
self.dedup.insert(genomes.to_vec(), id);
id
};
self.multi_dict_ids.push(id);
}
}
}
pub fn close(self) -> io::Result<()> {
let n = self.is_multi.len();
let n_singleton = self.singleton_values.len();
let n_multi = self.multi_dict_ids.len();
let n_distinct_multi = self.dict_sets_in_order.len();
let mut is_multi_b = PersistentRankSelectBitVecBuilder::new(n, &is_multi_path(&self.dir))?;
for (i, &v) in self.is_multi.iter().enumerate() {
is_multi_b.set(i, v);
}
is_multi_b.close()?;
let singleton_width = bit_width_for_range(self.n_cols as u64);
let mut singleton_b = PersistentFixedIntVecBuilder::new(n_singleton, singleton_width, &singleton_path(&self.dir))?;
for (i, &v) in self.singleton_values.iter().enumerate() {
singleton_b.set(i, v as u64);
}
singleton_b.close()?;
let multi_width = bit_width_for_range(n_distinct_multi as u64);
let mut multi_b = PersistentFixedIntVecBuilder::new(n_multi, multi_width, &multi_path(&self.dir))?;
for (i, &v) in self.multi_dict_ids.iter().enumerate() {
multi_b.set(i, v as u64);
}
multi_b.close()?;
let mut dict_values = Vec::new();
let mut offsets: Vec<u64> = Vec::with_capacity(n_distinct_multi);
for set in &self.dict_sets_in_order {
offsets.push(dict_values.len() as u64);
for &g in set {
write_varint(&mut dict_values, g);
}
}
let universe = dict_values.len() as u64 + 1;
let mut offsets_b = EliasFanoBuilder::new(offsets.len(), universe, &dict_offsets_base(&self.dir))?;
for &o in &offsets {
offsets_b.push(o);
}
offsets_b.close()?;
fs::write(dict_values_path(&self.dir), &dict_values)?;
SparseMeta { n, n_cols: self.n_cols, n_singleton, n_multi, n_distinct_multi }.save(&self.dir)?;
Ok(())
}
pub fn finish(self) -> io::Result<PersistentSparseBitMatrix> {
let dir = self.dir.clone();
self.close()?;
PersistentSparseBitMatrix::open(&dir)
}
/// Builds a sparse matrix from an already-built dense
/// [`PersistentBitMatrix`] — row-by-row transpose via
/// [`PersistentBitMatrix::fill_row`], for migrating an existing index.
pub fn build_from_dense(dense: &PersistentBitMatrix, 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 buf = vec![0u32; n_cols];
let mut genomes: Vec<u32> = Vec::new();
for slot in 0..n {
dense.fill_row(slot, &mut buf);
genomes.clear();
genomes.extend((0..n_cols).filter(|&c| buf[c] != 0).map(|c| c as u32));
builder.push_row(&genomes);
}
Ok(builder)
}
}
src/obicompactvec/src/eliasfano.rs
+133
View File
@@ -0,0 +1,133 @@
//! Elias-Fano encoding of a monotone (non-decreasing) `u64` sequence —
//! used only for the sparse presence matrix's dictionary offsets (see
//! `docmd/architecture/siblings.md` and the sparse-matrix plan): a
//! monotone, unbounded-magnitude sequence, the one component in that
//! design that genuinely needs this.
//!
//! Classic two-array construction: each value's low `l` bits are packed
//! at fixed width ([`crate::fixedintvec::PersistentFixedIntVec`]), its
//! high bits (`value >> l`) are unary-coded into a bitvector
//! ([`crate::rankselect::PersistentRankSelectBitVec`]) — a 1 bit at
//! position `(value >> l) + i` for the `i`-th value, which is itself
//! monotone non-decreasing whenever the input is, so it can be built by a
//! single forward pass with no lookahead. Decoding value `i`:
//! `((select1(i) - i) << l) | low[i]`.
//!
//! `l` is chosen as `floor(log2(universe / n))` (0 if `n >= universe`),
//! the standard choice that keeps the high bitvector's length — `n +
//! (universe >> l) + 1` — within a small constant factor of `n`.
use std::io;
use std::path::{Path, PathBuf};
use crate::fixedintvec::{PersistentFixedIntVec, PersistentFixedIntVecBuilder};
use crate::rankselect::{PersistentRankSelectBitVec, PersistentRankSelectBitVecBuilder};
fn low_bits_width(n: usize, universe: u64) -> u32 {
if n == 0 || universe <= n as u64 {
0
} else {
// floor(log2(universe / n))
(universe / n as u64).max(1).ilog2()
}
}
fn low_path(base: &Path) -> PathBuf {
let mut p = base.as_os_str().to_owned();
p.push(".efl");
PathBuf::from(p)
}
fn high_path(base: &Path) -> PathBuf {
let mut p = base.as_os_str().to_owned();
p.push(".efh");
PathBuf::from(p)
}
// ── EliasFano ────────────────────────────────────────────────────────────────
pub struct EliasFano {
low: PersistentFixedIntVec,
high: PersistentRankSelectBitVec,
n: usize,
low_width: u32,
}
impl EliasFano {
/// Opens a structure previously built at `base` (i.e. `{base}.efl` and
/// `{base}.efh`).
pub fn open(base: &Path) -> io::Result<Self> {
let low = PersistentFixedIntVec::open(&low_path(base))?;
let high = PersistentRankSelectBitVec::open(&high_path(base))?;
let n = low.len();
let low_width = low.width();
Ok(Self { low, high, n, low_width })
}
#[inline]
pub fn len(&self) -> usize { self.n }
#[inline]
pub fn is_empty(&self) -> bool { self.n == 0 }
/// The `i`-th value of the encoded sequence.
pub fn get(&self, i: usize) -> u64 {
debug_assert!(i < self.n);
let pos = self.high.select1(i as u64);
let high_part = pos as u64 - i as u64;
(high_part << self.low_width) | self.low.get(i)
}
}
// ── EliasFanoBuilder ─────────────────────────────────────────────────────────
pub struct EliasFanoBuilder {
low: PersistentFixedIntVecBuilder,
high: PersistentRankSelectBitVecBuilder,
low_width: u32,
n: usize,
next: usize,
last_value: u64,
}
impl EliasFanoBuilder {
/// `n` values will be [`push`](Self::push)ed, in non-decreasing order,
/// each `< universe`. Writes `{base}.efl` and `{base}.efh`.
pub fn new(n: usize, universe: u64, base: &Path) -> io::Result<Self> {
let low_width = low_bits_width(n, universe);
let high_len = n + (universe >> low_width) as usize + 1;
// `PersistentFixedIntVec` genuinely supports width 0 (every value
// 0, no storage) — the persisted width byte is then the source of
// truth on reopen, no separate bookkeeping of `low_width` needed.
let low = PersistentFixedIntVecBuilder::new(n, low_width, &low_path(base))?;
let high = PersistentRankSelectBitVecBuilder::new(high_len, &high_path(base))?;
Ok(Self { low, high, low_width, n, next: 0, last_value: 0 })
}
/// Appends the next value — must be `>= ` every previously pushed
/// value (monotone non-decreasing), and `< universe` as given to
/// [`new`](Self::new).
pub fn push(&mut self, value: u64) {
assert!(self.next < self.n, "push() called more than n={} times", self.n);
assert!(
self.next == 0 || value >= self.last_value,
"push({value}) breaks monotonicity: last value was {}", self.last_value
);
let low_mask = if self.low_width >= 64 { u64::MAX } else { (1u64 << self.low_width) - 1 };
self.low.set(self.next, value & low_mask);
let high_part = value >> self.low_width;
self.high.set(high_part as usize + self.next, true);
self.last_value = value;
self.next += 1;
}
pub fn close(self) -> io::Result<()> {
assert_eq!(self.next, self.n, "push() called {} times, expected n={}", self.next, self.n);
self.low.close()?;
self.high.close()
}
pub fn finish(self, base: &Path) -> io::Result<EliasFano> {
self.close()?;
EliasFano::open(base)
}
}
src/obicompactvec/src/fixedintvec.rs
+202
View File
@@ -0,0 +1,202 @@
//! Fixed-bit-width packed integer vector — `n` values, each exactly `width`
//! bits (1..=64, chosen once at construction from the actual data being
//! stored), no overflow scheme. Mirrors [`crate::bitvec::PersistentBitVec`]'s
//! mmap-backed reader/builder split, generalised from 1 bit/value to an
//! arbitrary width.
//!
//! Deliberately distinct from [`crate::reader::PersistentCompactIntVec`]
//! (`PCIV`): that format is "1 primary byte + overflow map for values ≥
//! 255", tuned for mostly-small-value distributions. This format is for
//! values roughly uniform over a known range (e.g. genome indices, or
//! `dict_id`s) — every value costs the same, `width` bits, by construction.
use std::fs::{File, OpenOptions};
use std::io::{self, Seek, SeekFrom, Write as _};
use std::path::{Path, PathBuf};
use memmap2::{Mmap, MmapMut};
const MAGIC: [u8; 4] = *b"PFIV";
// Header: magic(4) + width(1) + _pad(3) + n(8) = 16 bytes.
// Data starts at offset 16, u64-aligned (mmap base is page-aligned, 16 % 8 == 0).
const HEADER_SIZE: usize = 16;
/// Smallest bit width that can hold every value in `0..range` (`range`
/// itself excluded, i.e. the number of distinct values needed) — the
/// standard `ceil(log2(range))` sizing used throughout the sparse matrix
/// design, with the same `.max(1)` floor (a single-value range still needs
/// 1 bit, not 0, since `0..1` is a real distinct value).
pub fn bit_width_for_range(range: u64) -> u32 {
if range <= 1 {
1
} else {
(u64::BITS - (range - 1).leading_zeros()).max(1)
}
}
#[inline]
fn n_words_for(n: usize, width: u32) -> usize {
// +1 guard word: a value straddling the last word boundary reads/writes
// one bit into a word past the nominal bit count otherwise.
(n as u64 * width as u64).div_ceil(64) as usize + 1
}
#[inline]
fn mask(width: u32) -> u64 {
if width >= 64 { u64::MAX } else { (1u64 << width) - 1 }
}
// ── PersistentFixedIntVec ───────────────────────────────────────────────────
pub struct PersistentFixedIntVec {
mmap: Mmap,
n: usize,
width: u32,
path: PathBuf,
}
impl PersistentFixedIntVec {
pub fn open(path: &Path) -> io::Result<Self> {
let mmap = unsafe { Mmap::map(&File::open(path)?)? };
if mmap.len() < HEADER_SIZE {
return Err(io::Error::new(io::ErrorKind::InvalidData, "PFIV file too short"));
}
if mmap[0..4] != MAGIC {
return Err(io::Error::new(io::ErrorKind::InvalidData, "bad PFIV magic"));
}
let width = mmap[4] as u32;
let n = u64::from_le_bytes(mmap[8..16].try_into().unwrap()) as usize;
Ok(Self { mmap, n, width, path: path.to_path_buf() })
}
#[inline]
pub fn path(&self) -> &Path { &self.path }
#[inline]
pub fn len(&self) -> usize { self.n }
#[inline]
pub fn is_empty(&self) -> bool { self.n == 0 }
#[inline]
pub fn width(&self) -> u32 { self.width }
// SAFETY: mmap is page-aligned, HEADER_SIZE=16 divisible by 8 → u64-aligned.
#[inline]
fn data_words(&self) -> &[u64] {
let nw = n_words_for(self.n, self.width);
let ptr = self.mmap[HEADER_SIZE..].as_ptr() as *const u64;
unsafe { std::slice::from_raw_parts(ptr, nw) }
}
#[inline]
pub fn get(&self, slot: usize) -> u64 {
debug_assert!(slot < self.n);
get_packed(self.data_words(), slot, self.width)
}
}
#[inline]
fn get_packed(words: &[u64], slot: usize, width: u32) -> u64 {
let bit_offset = slot as u64 * width as u64;
let word_idx = (bit_offset / 64) as usize;
let bit_in_word = (bit_offset % 64) as u32;
let m = mask(width);
let lo = words[word_idx] >> bit_in_word;
if bit_in_word + width <= 64 {
lo & m
} else {
let hi = words[word_idx + 1] << (64 - bit_in_word);
(lo | hi) & m
}
}
#[inline]
fn set_packed(words: &mut [u64], slot: usize, width: u32, value: u64) {
let m = mask(width);
debug_assert!(value & !m == 0, "value {value} does not fit in {width} bits");
let bit_offset = slot as u64 * width as u64;
let word_idx = (bit_offset / 64) as usize;
let bit_in_word = (bit_offset % 64) as u32;
words[word_idx] &= !(m << bit_in_word);
words[word_idx] |= (value & m) << bit_in_word;
if bit_in_word + width > 64 {
let hi_bits = bit_in_word + width - 64;
let hi_mask = (1u64 << hi_bits) - 1;
words[word_idx + 1] &= !hi_mask;
words[word_idx + 1] |= value >> (64 - bit_in_word);
}
}
// ── PersistentFixedIntVecBuilder ────────────────────────────────────────────
pub struct PersistentFixedIntVecBuilder {
mmap: MmapMut,
n: usize,
width: u32,
path: PathBuf,
}
impl PersistentFixedIntVecBuilder {
/// `width` (1..=64) is chosen by the caller from the actual data being
/// stored — see [`bit_width_for_range`] — never hardcoded. `width == 0`
/// is allowed (every value is 0, e.g. a whole `EliasFano` sequence
/// packed entirely into its high bits) — `get` always returns 0,
/// `set` is a no-op, no storage beyond the header.
pub fn new(n: usize, width: u32, path: &Path) -> io::Result<Self> {
assert!(width <= 64, "width must be 0..=64, got {width}");
let file_size = HEADER_SIZE + n_words_for(n, width) * 8;
let mut file = OpenOptions::new()
.read(true).write(true).create(true).truncate(true)
.open(path)?;
file.write_all(&MAGIC)?;
file.write_all(&[width as u8, 0, 0, 0])?;
file.write_all(&(n as u64).to_le_bytes())?;
file.seek(SeekFrom::Start(0))?;
file.set_len(file_size as u64)?;
let mmap = unsafe { MmapMut::map_mut(&file)? };
Ok(Self { mmap, n, width, path: path.to_path_buf() })
}
#[inline]
pub fn len(&self) -> usize { self.n }
#[inline]
pub fn is_empty(&self) -> bool { self.n == 0 }
#[inline]
pub fn width(&self) -> u32 { self.width }
#[inline]
fn data_words_mut(&mut self) -> &mut [u64] {
let nw = n_words_for(self.n, self.width);
let ptr = self.mmap[HEADER_SIZE..].as_mut_ptr() as *mut u64;
unsafe { std::slice::from_raw_parts_mut(ptr, nw) }
}
#[inline]
fn data_words(&self) -> &[u64] {
let nw = n_words_for(self.n, self.width);
let ptr = self.mmap[HEADER_SIZE..].as_ptr() as *const u64;
unsafe { std::slice::from_raw_parts(ptr, nw) }
}
#[inline]
pub fn get(&self, slot: usize) -> u64 {
debug_assert!(slot < self.n);
get_packed(self.data_words(), slot, self.width)
}
#[inline]
pub fn set(&mut self, slot: usize, value: u64) {
debug_assert!(slot < self.n);
let width = self.width;
set_packed(self.data_words_mut(), slot, width, value);
}
pub fn close(self) -> io::Result<()> {
self.mmap.flush()
}
pub fn finish(self) -> io::Result<PersistentFixedIntVec> {
let path = self.path.clone();
self.close()?;
PersistentFixedIntVec::open(&path)
}
}
src/obicompactvec/src/lib.rs
+8 -2
View File
@@ -2,7 +2,10 @@ mod bitvec;
mod bitmatrix;
mod builder;
mod colgroup;
mod eliasfano;
mod fixedintvec;
mod format;
mod rankselect;
mod intmatrix;
mod layer_meta;
mod meta;
@@ -13,7 +16,10 @@ mod views;
pub mod traits;
pub use bitvec::{BitIter, PersistentBitVec, PersistentBitVecBuilder};
pub use bitmatrix::{PersistentBitMatrix, PersistentBitMatrixBuilder, pack_bit_matrix};
pub use fixedintvec::{PersistentFixedIntVec, PersistentFixedIntVecBuilder, bit_width_for_range};
pub use rankselect::{PersistentRankSelectBitVec, PersistentRankSelectBitVecBuilder};
pub use eliasfano::{EliasFano, EliasFanoBuilder};
pub use bitmatrix::{PersistentBitMatrix, PersistentBitMatrixBuilder, PersistentSparseBitMatrix, PersistentSparseBitMatrixBuilder, pack_bit_matrix};
pub use builder::PersistentCompactIntVecBuilder;
pub use colgroup::{ColGroup, FilterMask, MatrixGroupOps, eval_filter_mask};
pub use intmatrix::{PersistentCompactIntMatrix, PersistentCompactIntMatrixBuilder, pack_compact_int_matrix};
@@ -21,7 +27,7 @@ pub use layer_meta::LayerMeta;
pub use reader::{PersistentCompactIntVec, Iter as CompactIntVecIter};
pub use tempbitvec::{TempBitVec, TempBitVecBuilder};
pub use tempintvec::{TempCompactIntVec, TempCompactIntVecBuilder};
pub use traits::{BitPartials, ColumnWeights, CountPartials};
pub use traits::{BinaryMatrix, BitPartials, ColumnWeights, CountPartials};
pub use views::{BitSliceView, BitSliceIter, IntSliceView, IntSliceViewIter};
#[cfg(test)]
src/obicompactvec/src/rankselect.rs
+265
View File
@@ -0,0 +1,265 @@
//! Rank/select-capable bitvector — extends the crate's existing bit-count
//! (`count_ones`, a global reduction) with `rank1`/`rank0` (count of 1s/0s
//! in a prefix) and `select1` (position of the k-th 1 bit). Needed by two
//! consumers in the sparse presence-matrix design (see
//! `docmd/architecture/siblings.md` and the sparse-matrix plan): the
//! `is_multi` row-kind flag (needs rank, to locate a row's position within
//! whichever of the two split arrays it belongs to) and the Elias-Fano high
//! bits (needs select).
//!
//! Two-level structure, the standard succinct-bitvector approach: a
//! cumulative rank sampled every [`BLOCK_WORDS`] words, plus a linear scan
//! within the block (at most [`BLOCK_WORDS`] popcounts) for the remainder.
//! `select1` binary-searches the block samples, then
//! `common_traits::SelectInWord` locates the exact bit within the winning
//! word.
use std::fs::{File, OpenOptions};
use std::io::{self, Seek, SeekFrom, Write as _};
use std::path::{Path, PathBuf};
use common_traits::SelectInWord;
use memmap2::{Mmap, MmapMut};
const MAGIC: [u8; 4] = *b"PRSB";
/// Words per rank-sample block — 8 words = 512 bits, one cache line's
/// worth of linear-scan work for the within-block remainder.
const BLOCK_WORDS: usize = 8;
// Header: magic(4) + _pad(4) + n(8) + total_ones(8) = 24 bytes (8-aligned).
const HEADER_SIZE: usize = 24;
#[inline]
fn n_words(n: usize) -> usize {
n.div_ceil(64)
}
#[inline]
fn n_blocks(n: usize) -> usize {
n_words(n).div_ceil(BLOCK_WORDS)
}
#[inline]
fn bits_bytes(n: usize) -> usize {
n_words(n) * 8
}
#[inline]
fn blocks_bytes(n: usize) -> usize {
n_blocks(n) * 8
}
// ── PersistentRankSelectBitVec ──────────────────────────────────────────────
pub struct PersistentRankSelectBitVec {
mmap: Mmap,
n: usize,
total_ones: u64,
path: PathBuf,
}
impl PersistentRankSelectBitVec {
pub fn open(path: &Path) -> io::Result<Self> {
let mmap = unsafe { Mmap::map(&File::open(path)?)? };
if mmap.len() < HEADER_SIZE {
return Err(io::Error::new(io::ErrorKind::InvalidData, "PRSB file too short"));
}
if mmap[0..4] != MAGIC {
return Err(io::Error::new(io::ErrorKind::InvalidData, "bad PRSB magic"));
}
let n = u64::from_le_bytes(mmap[8..16].try_into().unwrap()) as usize;
let total_ones = u64::from_le_bytes(mmap[16..24].try_into().unwrap());
Ok(Self { mmap, n, total_ones, path: path.to_path_buf() })
}
#[inline]
pub fn path(&self) -> &Path { &self.path }
#[inline]
pub fn len(&self) -> usize { self.n }
#[inline]
pub fn is_empty(&self) -> bool { self.n == 0 }
#[inline]
pub fn count_ones(&self) -> u64 { self.total_ones }
#[inline]
pub fn count_zeros(&self) -> u64 { self.n as u64 - self.total_ones }
// SAFETY: mmap is page-aligned, HEADER_SIZE=24 divisible by 8 → u64-aligned.
#[inline]
fn bit_words(&self) -> &[u64] {
let nw = n_words(self.n);
let ptr = self.mmap[HEADER_SIZE..].as_ptr() as *const u64;
unsafe { std::slice::from_raw_parts(ptr, nw) }
}
#[inline]
fn block_ranks(&self) -> &[u64] {
let nb = n_blocks(self.n);
let off = HEADER_SIZE + bits_bytes(self.n);
let ptr = self.mmap[off..].as_ptr() as *const u64;
unsafe { std::slice::from_raw_parts(ptr, nb) }
}
#[inline]
pub fn get(&self, pos: usize) -> bool {
debug_assert!(pos < self.n);
(self.bit_words()[pos >> 6] >> (pos & 63)) & 1 != 0
}
/// Number of 1 bits in `[0, pos)`.
pub fn rank1(&self, pos: usize) -> u64 {
debug_assert!(pos <= self.n);
if pos == 0 {
return 0;
}
let words = self.bit_words();
let block_ranks = self.block_ranks();
let word_idx = (pos - 1) / 64;
let block_idx = word_idx / BLOCK_WORDS;
let mut rank = block_ranks[block_idx];
let block_start_word = block_idx * BLOCK_WORDS;
for w in &words[block_start_word..word_idx] {
rank += w.count_ones() as u64;
}
let bit_in_word = pos - word_idx * 64;
let last = words[word_idx];
let masked = if bit_in_word >= 64 { last } else { last & ((1u64 << bit_in_word) - 1) };
rank += masked.count_ones() as u64;
rank
}
/// Number of 0 bits in `[0, pos)`.
#[inline]
pub fn rank0(&self, pos: usize) -> u64 {
pos as u64 - self.rank1(pos)
}
/// Position of the `k`-th (0-indexed) 1 bit. Panics if fewer than
/// `k + 1` ones exist.
pub fn select1(&self, k: u64) -> usize {
assert!(k < self.total_ones, "select1({k}) out of range: only {} ones", self.total_ones);
let block_ranks = self.block_ranks();
let words = self.bit_words();
// Binary search: last block whose cumulative rank is <= k.
let mut lo = 0usize;
let mut hi = block_ranks.len();
while lo + 1 < hi {
let mid = lo + (hi - lo) / 2;
if block_ranks[mid] <= k { lo = mid; } else { hi = mid; }
}
let block_idx = lo;
let mut remaining = k - block_ranks[block_idx];
let start_word = block_idx * BLOCK_WORDS;
let end_word = (start_word + BLOCK_WORDS).min(words.len());
for (i, &w) in words[start_word..end_word].iter().enumerate() {
let c = w.count_ones() as u64;
if remaining < c {
return (start_word + i) * 64 + w.select_in_word(remaining as usize);
}
remaining -= c;
}
unreachable!("select1({k}): ran off the end of its own block — rank/select index inconsistent");
}
}
// ── PersistentRankSelectBitVecBuilder ───────────────────────────────────────
pub struct PersistentRankSelectBitVecBuilder {
mmap: MmapMut,
n: usize,
path: PathBuf,
}
impl PersistentRankSelectBitVecBuilder {
pub fn new(n: usize, path: &Path) -> io::Result<Self> {
let file_size = HEADER_SIZE + bits_bytes(n) + blocks_bytes(n);
let mut file = OpenOptions::new()
.read(true).write(true).create(true).truncate(true)
.open(path)?;
file.write_all(&MAGIC)?;
file.write_all(&[0u8; 4])?;
file.write_all(&(n as u64).to_le_bytes())?;
file.write_all(&0u64.to_le_bytes())?; // total_ones, patched in close()
file.seek(SeekFrom::Start(0))?;
file.set_len(file_size as u64)?;
let mmap = unsafe { MmapMut::map_mut(&file)? };
Ok(Self { mmap, n, path: path.to_path_buf() })
}
#[inline]
pub fn len(&self) -> usize { self.n }
#[inline]
pub fn is_empty(&self) -> bool { self.n == 0 }
#[inline]
fn bit_words_mut(&mut self) -> &mut [u64] {
let nw = n_words(self.n);
let ptr = self.mmap[HEADER_SIZE..].as_mut_ptr() as *mut u64;
unsafe { std::slice::from_raw_parts_mut(ptr, nw) }
}
#[inline]
pub fn get(&self, pos: usize) -> bool {
debug_assert!(pos < self.n);
let nw = n_words(self.n);
let ptr = self.mmap[HEADER_SIZE..].as_ptr() as *const u64;
let words = unsafe { std::slice::from_raw_parts(ptr, nw) };
(words[pos >> 6] >> (pos & 63)) & 1 != 0
}
#[inline]
pub fn set(&mut self, pos: usize, value: bool) {
debug_assert!(pos < self.n);
let bit = 1u64 << (pos & 63);
let words = self.bit_words_mut();
if value {
words[pos >> 6] |= bit;
} else {
words[pos >> 6] &= !bit;
}
}
/// Computes the block-rank index and total-ones count from the bits
/// written so far, then flushes. Must run after every `set()` call —
/// the index is a function of the final bit pattern, not maintainable
/// incrementally through arbitrary overwrites.
pub fn close(mut self) -> io::Result<()> {
let nb = n_blocks(self.n);
let nw = n_words(self.n);
let mut block_ranks = vec![0u64; nb];
let mut running = 0u64;
{
let words = self.bit_words_mut();
for (block_idx, block_rank) in block_ranks.iter_mut().enumerate() {
*block_rank = running;
let start = block_idx * BLOCK_WORDS;
let end = (start + BLOCK_WORDS).min(nw);
for &w in &words[start..end] {
running += w.count_ones() as u64;
}
}
}
self.mmap[16..24].copy_from_slice(&running.to_le_bytes());
let block_off = HEADER_SIZE + bits_bytes(self.n);
let block_bytes = u64_slice_to_le_bytes(&block_ranks);
self.mmap[block_off..block_off + block_bytes.len()].copy_from_slice(&block_bytes);
self.mmap.flush()
}
pub fn finish(self) -> io::Result<PersistentRankSelectBitVec> {
let path = self.path.clone();
self.close()?;
PersistentRankSelectBitVec::open(&path)
}
}
/// Minimal, local `u64` slice -> byte vec conversion (little-endian,
/// matching every other on-disk format in this crate).
fn u64_slice_to_le_bytes(values: &[u64]) -> Vec<u8> {
let mut out = Vec::with_capacity(values.len() * 8);
for &v in values {
out.extend_from_slice(&v.to_le_bytes());
}
out
}
src/obicompactvec/src/tests/eliasfano.rs
+114
View File
@@ -0,0 +1,114 @@
use tempfile::tempdir;
use crate::{EliasFano, EliasFanoBuilder};
fn build(values: &[u64], universe: u64) -> (tempfile::TempDir, EliasFano) {
let dir = tempdir().unwrap();
let base = dir.path().join("test.ef");
let mut b = EliasFanoBuilder::new(values.len(), universe, &base).unwrap();
for &v in values {
b.push(v);
}
let ef = b.finish(&base).unwrap();
(dir, ef)
}
#[test]
fn empty() {
let (_dir, ef) = build(&[], 1000);
assert_eq!(ef.len(), 0);
assert!(ef.is_empty());
}
#[test]
fn small_monotone_sequence() {
let values = [0u64, 3, 3, 7, 42, 42, 42, 100];
let (_dir, ef) = build(&values, 1000);
assert_eq!(ef.len(), values.len());
for (i, &expected) in values.iter().enumerate() {
assert_eq!(ef.get(i), expected, "index {i}");
}
}
#[test]
fn all_equal_values() {
let values = [5u64; 20];
let (_dir, ef) = build(&values, 100);
for i in 0..values.len() {
assert_eq!(ef.get(i), 5);
}
}
#[test]
fn strictly_increasing() {
let values: Vec<u64> = (0..500).map(|i| i * 3).collect();
let universe = values.last().unwrap() + 1;
let (_dir, ef) = build(&values, universe);
for (i, &expected) in values.iter().enumerate() {
assert_eq!(ef.get(i), expected, "index {i}");
}
}
#[test]
fn realistic_byte_offsets() {
// Simulates the actual use case: cumulative byte offsets into a
// varint blob, one entry per distinct multi-genome set, strictly
// increasing by a small, variable amount each time (1-2 bytes/index,
// several indices per set).
let mut offset = 0u64;
let mut values = Vec::new();
let mut rng_state = 12345u64;
for _ in 0..10_000 {
values.push(offset);
// simple xorshift for a deterministic, dependency-free "random" gap
rng_state ^= rng_state << 13;
rng_state ^= rng_state >> 7;
rng_state ^= rng_state << 17;
offset += 1 + (rng_state % 8); // 1..=8 bytes per set, realistic for n_cols=91
}
let universe = offset + 1;
let (_dir, ef) = build(&values, universe);
for (i, &expected) in values.iter().enumerate() {
assert_eq!(ef.get(i), expected, "index {i}");
}
}
#[test]
fn universe_smaller_than_n_gives_zero_width_low() {
// n=100 values drawn from a universe of only 10 — low_bits_width
// should come out 0 (falls back to pure unary/high-bits encoding).
let values: Vec<u64> = (0..100).map(|i| i / 10).collect(); // 0,0,..,0,1,1,..,9,9
let (_dir, ef) = build(&values, 10);
for (i, &expected) in values.iter().enumerate() {
assert_eq!(ef.get(i), expected, "index {i}");
}
}
#[test]
#[should_panic(expected = "monotonicity")]
fn push_out_of_order_panics() {
let dir = tempdir().unwrap();
let base = dir.path().join("test.ef");
let mut b = EliasFanoBuilder::new(3, 100, &base).unwrap();
b.push(10);
b.push(5); // decreasing — must panic
}
#[test]
fn reopen_after_close_matches_original() {
let values: Vec<u64> = (0..2000).map(|i| i * 5 + (i % 3)).collect();
let universe = *values.last().unwrap() + 1;
let dir = tempdir().unwrap();
let base = dir.path().join("test.ef");
{
let mut b = EliasFanoBuilder::new(values.len(), universe, &base).unwrap();
for &v in &values {
b.push(v);
}
b.close().unwrap();
} // builder + mmaps fully dropped here
let ef = EliasFano::open(&base).unwrap();
for (i, &expected) in values.iter().enumerate() {
assert_eq!(ef.get(i), expected, "index {i}");
}
}
src/obicompactvec/src/tests/fixedintvec.rs
+112
View File
@@ -0,0 +1,112 @@
use tempfile::tempdir;
use crate::{PersistentFixedIntVec, PersistentFixedIntVecBuilder, bit_width_for_range};
#[test]
fn bit_width_for_range_matches_expectations() {
assert_eq!(bit_width_for_range(0), 1);
assert_eq!(bit_width_for_range(1), 1);
assert_eq!(bit_width_for_range(2), 1);
assert_eq!(bit_width_for_range(3), 2);
assert_eq!(bit_width_for_range(4), 2);
assert_eq!(bit_width_for_range(91), 7);
assert_eq!(bit_width_for_range(128), 7);
assert_eq!(bit_width_for_range(129), 8);
assert_eq!(bit_width_for_range(460_591), 19);
}
fn roundtrip(width: u32, values: &[u64]) -> Vec<u64> {
let dir = tempdir().unwrap();
let path = dir.path().join("test.pfiv");
let mut b = PersistentFixedIntVecBuilder::new(values.len(), width, &path).unwrap();
for (i, &v) in values.iter().enumerate() {
b.set(i, v);
}
b.close().unwrap();
let r = PersistentFixedIntVec::open(&path).unwrap();
assert_eq!(r.len(), values.len());
assert_eq!(r.width(), width);
(0..values.len()).map(|s| r.get(s)).collect()
}
#[test]
fn width_1_roundtrip() {
let values = [0u64, 1, 1, 0, 1];
assert_eq!(roundtrip(1, &values), values);
}
#[test]
fn width_7_roundtrip_genome_indices() {
// 91-genome-scale values, width=7 (matches bit_width_for_range(91)).
let values: Vec<u64> = (0..91).collect();
assert_eq!(roundtrip(7, &values), values);
}
#[test]
fn width_19_roundtrip_dict_ids() {
// Values crossing many word boundaries at a non-power-of-two width —
// exactly the case that breaks a naive byte/word-aligned packer.
let values: Vec<u64> = (0..2000).map(|i| (i * 37) % 460_591).collect();
assert_eq!(roundtrip(19, &values), values);
}
#[test]
fn width_64_roundtrip() {
let values = [0u64, u64::MAX, 1, u64::MAX - 1, 1 << 40];
assert_eq!(roundtrip(64, &values), values);
}
#[test]
fn values_spanning_word_boundary() {
// width=19: slot 3's bits start at bit 57, spans into the next word.
let values: Vec<u64> = vec![524_287, 0, 0, 500_000, 1];
assert_eq!(roundtrip(19, &values), values);
}
#[test]
fn mutation_via_set_overwrites() {
let dir = tempdir().unwrap();
let path = dir.path().join("test.pfiv");
let mut b = PersistentFixedIntVecBuilder::new(3, 10, &path).unwrap();
b.set(0, 5);
b.set(1, 1000);
b.set(2, 3);
b.set(1, 42); // overwrite
b.close().unwrap();
let r = PersistentFixedIntVec::open(&path).unwrap();
assert_eq!(r.get(0), 5);
assert_eq!(r.get(1), 42);
assert_eq!(r.get(2), 3);
}
#[test]
fn all_zero_by_default() {
let dir = tempdir().unwrap();
let path = dir.path().join("test.pfiv");
let b = PersistentFixedIntVecBuilder::new(50, 13, &path).unwrap();
b.close().unwrap();
let r = PersistentFixedIntVec::open(&path).unwrap();
for i in 0..50 {
assert_eq!(r.get(i), 0, "slot {i}");
}
}
#[test]
fn reopen_after_close_matches_original() {
// Explicit disk round-trip: drop the builder/mmap entirely, reopen
// from a fresh `Mmap::map` call, not the same in-memory handle.
let values: Vec<u64> = (0..5000).map(|i| (i * 97) % 8192).collect();
let dir = tempdir().unwrap();
let path = dir.path().join("test.pfiv");
{
let mut b = PersistentFixedIntVecBuilder::new(values.len(), 13, &path).unwrap();
for (i, &v) in values.iter().enumerate() {
b.set(i, v);
}
b.close().unwrap();
} // builder + mmap fully dropped here
let r = PersistentFixedIntVec::open(&path).unwrap();
for (i, &expected) in values.iter().enumerate() {
assert_eq!(r.get(i), expected, "slot {i}");
}
}
src/obicompactvec/src/tests/mod.rs
+4
View File
@@ -1,7 +1,11 @@
mod bitmatrix;
mod bitvec;
mod colgroup;
mod eliasfano;
mod fixedintvec;
mod intmatrix;
mod rankselect;
mod sparse;
use tempfile::tempdir;
src/obicompactvec/src/tests/rankselect.rs
+167
View File
@@ -0,0 +1,167 @@
use tempfile::tempdir;
use crate::{PersistentRankSelectBitVec, PersistentRankSelectBitVecBuilder};
fn build(bits: &[bool]) -> (tempfile::TempDir, PersistentRankSelectBitVec) {
let dir = tempdir().unwrap();
let path = dir.path().join("test.prsb");
let mut b = PersistentRankSelectBitVecBuilder::new(bits.len(), &path).unwrap();
for (i, &v) in bits.iter().enumerate() {
b.set(i, v);
}
b.close().unwrap();
let r = PersistentRankSelectBitVec::open(&path).unwrap();
(dir, r)
}
fn naive_rank1(bits: &[bool], pos: usize) -> u64 {
bits[..pos].iter().filter(|&&b| b).count() as u64
}
fn naive_select1(bits: &[bool], k: u64) -> usize {
bits.iter().enumerate().filter(|&(_, &b)| b).nth(k as usize).unwrap().0
}
#[test]
fn get_matches_input() {
let bits = [true, false, true, true, false, false, true];
let (_dir, r) = build(&bits);
for (i, &expected) in bits.iter().enumerate() {
assert_eq!(r.get(i), expected, "bit {i}");
}
}
#[test]
fn count_ones_matches_naive() {
let bits: Vec<bool> = (0..1000).map(|i| i % 3 == 0).collect();
let (_dir, r) = build(&bits);
assert_eq!(r.count_ones(), bits.iter().filter(|&&b| b).count() as u64);
assert_eq!(r.count_zeros(), bits.iter().filter(|&&b| !b).count() as u64);
}
#[test]
fn rank1_matches_naive_small() {
let bits = [true, false, true, true, false, false, true, true, false, true];
let (_dir, r) = build(&bits);
for pos in 0..=bits.len() {
assert_eq!(r.rank1(pos), naive_rank1(&bits, pos), "rank1({pos})");
}
}
#[test]
fn rank1_matches_naive_multi_block() {
// BLOCK_WORDS=8 words=512 bits — exercise several blocks, an odd
// total length, and a non-uniform bit pattern.
let n = 3000;
let bits: Vec<bool> = (0..n).map(|i| (i * 7 + 3) % 11 == 0).collect();
let (_dir, r) = build(&bits);
for pos in (0..=n).step_by(37) {
assert_eq!(r.rank1(pos), naive_rank1(&bits, pos), "rank1({pos})");
}
assert_eq!(r.rank1(n), naive_rank1(&bits, n));
}
#[test]
fn rank0_is_complement_of_rank1() {
let bits: Vec<bool> = (0..777).map(|i| i % 5 < 2).collect();
let (_dir, r) = build(&bits);
for pos in (0..=bits.len()).step_by(13) {
assert_eq!(r.rank0(pos), pos as u64 - r.rank1(pos));
}
}
#[test]
fn select1_matches_naive_small() {
let bits = [true, false, true, true, false, false, true, true, false, true];
let (_dir, r) = build(&bits);
let n_ones = bits.iter().filter(|&&b| b).count() as u64;
for k in 0..n_ones {
assert_eq!(r.select1(k), naive_select1(&bits, k), "select1({k})");
}
}
#[test]
fn select1_matches_naive_multi_block() {
let n = 3000;
let bits: Vec<bool> = (0..n).map(|i| (i * 13 + 5) % 17 == 0).collect();
let (_dir, r) = build(&bits);
let n_ones = bits.iter().filter(|&&b| b).count() as u64;
for k in (0..n_ones).step_by(23) {
assert_eq!(r.select1(k), naive_select1(&bits, k), "select1({k})");
}
}
#[test]
fn rank_select_round_trip() {
// For every one-bit's position p, select1(rank1(p)) == p.
let n = 2000;
let bits: Vec<bool> = (0..n).map(|i| (i * 31 + 1) % 9 == 0).collect();
let (_dir, r) = build(&bits);
for (p, &is_one) in bits.iter().enumerate() {
if is_one {
let k = r.rank1(p);
assert_eq!(r.select1(k), p, "position {p}, rank {k}");
}
}
}
#[test]
#[should_panic]
fn select1_out_of_range_panics() {
let bits = [true, false, false];
let (_dir, r) = build(&bits);
r.select1(1); // only one 1-bit (k=0 valid), k=1 must panic
}
#[test]
fn all_zeros() {
let bits = vec![false; 200];
let (_dir, r) = build(&bits);
assert_eq!(r.count_ones(), 0);
assert_eq!(r.rank1(200), 0);
assert_eq!(r.rank0(200), 200);
}
#[test]
fn all_ones() {
let bits = vec![true; 200];
let (_dir, r) = build(&bits);
assert_eq!(r.count_ones(), 200);
assert_eq!(r.rank1(200), 200);
for k in 0..200 {
assert_eq!(r.select1(k), k as usize);
}
}
#[test]
fn reopen_after_close_matches_original() {
let n = 5000;
let bits: Vec<bool> = (0..n).map(|i| (i * 41 + 7) % 13 == 0).collect();
let dir = tempdir().unwrap();
let path = dir.path().join("test.prsb");
{
let mut b = PersistentRankSelectBitVecBuilder::new(n, &path).unwrap();
for (i, &v) in bits.iter().enumerate() {
b.set(i, v);
}
b.close().unwrap();
} // builder + mmap fully dropped here
let r = PersistentRankSelectBitVec::open(&path).unwrap();
assert_eq!(r.count_ones(), bits.iter().filter(|&&b| b).count() as u64);
for pos in (0..=n).step_by(17) {
assert_eq!(r.rank1(pos), naive_rank1(&bits, pos), "rank1({pos})");
}
}
#[test]
fn select1_run_of_ten_then_gap_then_run() {
// Mirrors the exact bit pattern from the failing EliasFano case:
// positions 0..=9 set, 10 clear, 11..=20 set.
let mut bits = vec![false; 30];
for p in 0..=9 { bits[p] = true; }
for p in 11..=20 { bits[p] = true; }
let (_dir, r) = build(&bits);
assert_eq!(r.select1(9), 9, "select1(9)");
assert_eq!(r.select1(10), 11, "select1(10)");
assert_eq!(r.select1(11), 12, "select1(11)");
}
src/obicompactvec/src/tests/sparse.rs
+251
View File
@@ -0,0 +1,251 @@
use tempfile::tempdir;
use crate::{BinaryMatrix, PersistentBitMatrix, PersistentBitMatrixBuilder, PersistentSparseBitMatrix, PersistentSparseBitMatrixBuilder};
/// Builds a dense `PersistentBitMatrix` from column-major `bool` data —
/// mirrors `tests/bitmatrix.rs`'s own `make_matrix` helper.
fn make_dense(cols: &[&[bool]]) -> (tempfile::TempDir, PersistentBitMatrix) {
let n = cols.first().map_or(0, |c| c.len());
let dir = tempdir().unwrap();
let presence = dir.path().join("presence");
let mut b = PersistentBitMatrixBuilder::new(n, &presence).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 = PersistentBitMatrix::open(dir.path()).unwrap();
(dir, m)
}
/// Builds a sparse matrix directly from row-major `bool` data (one slice
/// per row, `n_cols` bools each) — the natural input shape for this type.
fn make_sparse(rows: &[&[bool]], n_cols: usize) -> (tempfile::TempDir, PersistentSparseBitMatrix) {
let dir = tempdir().unwrap();
let sparse_dir = dir.path().join("sparse");
let mut b = PersistentSparseBitMatrixBuilder::new(rows.len(), n_cols, &sparse_dir).unwrap();
let mut genomes = Vec::new();
for row in rows {
genomes.clear();
genomes.extend((0..n_cols).filter(|&c| row[c]).map(|c| c as u32));
b.push_row(&genomes);
}
let m = b.finish().unwrap();
(dir, m)
}
fn row_as_bool(m: &PersistentSparseBitMatrix, slot: usize) -> Vec<bool> {
m.row(slot).to_vec()
}
#[test]
fn basic_roundtrip_singletons_and_multi() {
// 5 rows, 4 genomes: mix of singleton, multi-genome, and one repeated
// multi-genome set (dedup should collapse it to one dictionary entry).
let rows: Vec<&[bool]> = vec![
&[true, false, false, false], // singleton: genome 0
&[false, false, true, false], // singleton: genome 2
&[true, true, false, false], // multi: {0,1}
&[false, false, true, true], // multi: {2,3}
&[true, true, false, false], // multi: {0,1} again — should dedup
];
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!(row_as_bool(&m, slot), expected, "row {slot}");
}
}
#[test]
fn fill_row_matches_row() {
let rows: Vec<&[bool]> = vec![
&[true, false, true],
&[false, true, false],
&[true, true, true],
];
let (_dir, m) = make_sparse(&rows, 3);
let mut buf = vec![0u32; 3];
for slot in 0..3 {
m.fill_row(slot, &mut buf);
let via_fill: Vec<bool> = buf.iter().map(|&v| v != 0).collect();
assert_eq!(via_fill, row_as_bool(&m, slot), "slot {slot}");
}
}
#[test]
fn dense_to_sparse_matches_dense_on_real_shaped_data() {
// Column-major dense fixture: 4 genomes (columns), 6 k-mer slots (rows),
// deliberately including singletons, a repeated multi-genome set, and
// one all-genomes row.
let col0 = [true, false, true, false, true, true];
let col1 = [false, false, true, false, true, false];
let col2 = [false, true, false, false, false, true];
let col3 = [false, false, false, true, false, true];
let (_dense_dir, dense) = make_dense(&[&col0, &col1, &col2, &col3]);
let sparse_root = tempdir().unwrap();
let sparse_dir = sparse_root.path().join("sparse");
let sparse = PersistentSparseBitMatrixBuilder::build_from_dense(&dense, &sparse_dir)
.unwrap()
.finish()
.unwrap();
assert_eq!(sparse.n(), dense.n());
assert_eq!(sparse.n_cols(), dense.n_cols());
for slot in 0..dense.n() {
assert_eq!(
row_as_bool(&sparse, slot), &*dense.row(slot),
"slot {slot}: dense vs sparse disagree"
);
}
}
#[test]
fn count_ones_matches_dense() {
let col0 = [true, false, true, true];
let col1 = [false, true, true, false];
let (_dense_dir, dense) = make_dense(&[&col0, &col1]);
let sparse_root = tempdir().unwrap();
let sparse_dir = sparse_root.path().join("sparse");
let sparse = PersistentSparseBitMatrixBuilder::build_from_dense(&dense, &sparse_dir)
.unwrap()
.finish()
.unwrap();
let dense_counts = dense.count_ones();
let sparse_counts = sparse.count_ones();
assert_eq!(sparse_counts.to_vec(), dense_counts.to_vec());
}
#[test]
fn fill_sub_matrix_matches_dense() {
let col0 = [true, false, true, false, true];
let col1 = [false, true, true, true, false];
let col2 = [true, true, false, false, true];
let (_dense_dir, dense) = make_dense(&[&col0, &col1, &col2]);
let sparse_root = tempdir().unwrap();
let sparse_dir = sparse_root.path().join("sparse");
let sparse = PersistentSparseBitMatrixBuilder::build_from_dense(&dense, &sparse_dir)
.unwrap()
.finish()
.unwrap();
let slots = [4usize, 0, 2];
let dense_sub = dense.sub_matrix(&slots);
let mut sparse_sub: Vec<Vec<bool>> = vec![Vec::new(); dense.n_cols()];
sparse.fill_sub_matrix(&slots, &mut sparse_sub);
assert_eq!(sparse_sub, dense_sub);
}
#[test]
fn single_genome_matrix() {
// n_cols=1 — every row is necessarily a singleton (cardinality 1 or
// 0), the dictionary/multi-array stay entirely empty.
let rows: Vec<&[bool]> = vec![&[true], &[false], &[true]];
let (_dir, m) = make_sparse(&rows, 1);
for (slot, &expected) in rows.iter().enumerate() {
assert_eq!(row_as_bool(&m, slot), expected, "row {slot}");
}
}
#[test]
fn reopen_after_close_matches_original() {
// Explicit disk round-trip: drop the builder and every mmap it holds
// entirely, reopen fresh from `PersistentSparseBitMatrix::open`.
let rows: Vec<Vec<bool>> = (0..500)
.map(|i| (0..37).map(|c| (i * 7 + c * 3) % 11 == 0).collect())
.collect();
let dir = tempdir().unwrap();
let sparse_dir = dir.path().join("sparse");
{
let mut b = PersistentSparseBitMatrixBuilder::new(rows.len(), 37, &sparse_dir).unwrap();
let mut genomes = Vec::new();
for row in &rows {
genomes.clear();
genomes.extend((0..37).filter(|&c| row[c]).map(|c| c as u32));
b.push_row(&genomes);
}
b.close().unwrap();
} // builder + every mmap fully dropped here
let m = PersistentSparseBitMatrix::open(&sparse_dir).unwrap();
assert_eq!(m.n(), 500);
assert_eq!(m.n_cols(), 37);
for (slot, expected) in rows.iter().enumerate() {
assert_eq!(&row_as_bool(&m, slot), expected, "row {slot}");
}
}
/// Exercises `BinaryMatrix` generically over both concrete types — proves
/// they're actually interchangeable at that call site, not just
/// individually correct.
fn sum_via_trait(m: &dyn BinaryMatrix, slot: usize) -> usize {
m.row(slot).iter().filter(|&&b| b).count()
}
#[test]
fn binary_matrix_trait_interchangeable_dense_and_sparse() {
let col0 = [true, false, true, true];
let col1 = [false, true, true, false];
let col2 = [true, true, false, true];
let (_dense_dir, dense) = make_dense(&[&col0, &col1, &col2]);
let sparse_root = tempdir().unwrap();
let sparse_dir = sparse_root.path().join("sparse");
let sparse = PersistentSparseBitMatrixBuilder::build_from_dense(&dense, &sparse_dir)
.unwrap()
.finish()
.unwrap();
assert_eq!(BinaryMatrix::n(&dense), BinaryMatrix::n(&sparse));
assert_eq!(BinaryMatrix::n_cols(&dense), BinaryMatrix::n_cols(&sparse));
for slot in 0..BinaryMatrix::n(&dense) {
assert_eq!(sum_via_trait(&dense, slot), sum_via_trait(&sparse, slot), "slot {slot}");
}
let dense_counts = BinaryMatrix::count_ones(&dense);
let sparse_counts = BinaryMatrix::count_ones(&sparse);
assert_eq!(dense_counts.to_vec(), sparse_counts.to_vec());
}
#[test]
fn reopen_large_real_shaped_data_with_heavy_dedup() {
// Larger, more realistic case: 91-genome scale, plenty of repeated
// multi-genome sets (dedup exercised for real), reopened from disk.
let n_cols = 91;
let n_rows = 5000;
let rows: Vec<Vec<bool>> = (0..n_rows)
.map(|i| {
let mut row = vec![false; n_cols];
match i % 5 {
0 => { row[i % n_cols] = true; } // singleton, varies per row
1 => { row[3] = true; row[7] = true; } // repeated multi-set A
2 => { row[3] = true; row[7] = true; } // same set A again
3 => { row[10] = true; row[20] = true; row[30] = true; } // repeated multi-set B
_ => { row[(i * 13) % n_cols] = true; row[(i * 29) % n_cols] = true; } // varying pairs
}
row
})
.collect();
let dir = tempdir().unwrap();
let sparse_dir = dir.path().join("sparse");
{
let mut b = PersistentSparseBitMatrixBuilder::new(n_rows, n_cols, &sparse_dir).unwrap();
let mut genomes = Vec::new();
for row in &rows {
genomes.clear();
genomes.extend((0..n_cols).filter(|&c| row[c]).map(|c| c as u32));
b.push_row(&genomes);
}
b.close().unwrap();
}
let m = PersistentSparseBitMatrix::open(&sparse_dir).unwrap();
for (slot, expected) in rows.iter().enumerate() {
assert_eq!(&row_as_bool(&m, slot), expected, "row {slot}");
}
}
src/obicompactvec/src/traits.rs
+28
View File
@@ -1,5 +1,33 @@
use ndarray::{Array1, Array2};
/// Minimal shared surface between `PersistentBitMatrix` (dense) and
/// `PersistentSparseBitMatrix` (row-major, deduplicated) — exactly what
/// real consumers use today (`obikphylo::siblings::cache::Mat`), not the
/// two types' full individual APIs. Column-oriented operations
/// (`col`/`col_view`, the `BitPartials`/`ColumnWeights` distance-matrix
/// traits above) are *not* part of this trait — `PersistentSparseBitMatrix`
/// only offers a naive, row-scanning `count_ones` for now (see
/// `docmd/architecture/siblings.md` and the sparse-matrix design plan,
/// "Explicitly deferred": a row-major co-occurrence rewrite of the
/// pairwise distance matrices is future work, not part of this trait).
pub trait BinaryMatrix {
/// Number of rows (k-mer slots).
fn n(&self) -> usize;
/// Number of columns (genomes).
fn n_cols(&self) -> usize;
/// One row's presence values, one `bool` per genome.
fn row(&self, slot: usize) -> Box<[bool]>;
/// Like [`row`](Self::row), filling a caller-provided `0`/`1` buffer
/// instead of allocating.
fn fill_row(&self, slot: usize, buf: &mut [u32]);
/// Extracts a sub-matrix at `slots`, column-first: `out[c]` holds
/// column `c`'s values at `slots`, in `slots` order. `out.len()` must
/// equal `n_cols()`.
fn fill_sub_matrix(&self, slots: &[usize], out: &mut [Vec<bool>]);
/// Per-genome k-mer totals.
fn count_ones(&self) -> Array1<u64>;
}
/// Convert a Jaccard distance matrix (`1 - J`) into a Mash distance matrix, per
/// https://mash.readthedocs.io/en/latest/distances.html:
/// `D = -1/k * ln(2J / (1+J))`.
src/obikphylo/src/siblings/tests.rs
+256
View File
@@ -535,3 +535,259 @@ fn diag_plant_index_cardinality_distribution() {
);
}
}
#[test]
#[ignore]
fn bench_sparse_matrix_inmemory_construction_ram() {
use std::collections::HashMap;
use std::time::Instant;
let idx = KmerIndex::open("/Users/coissac/travail/obiskim/data/phyloalps/phyloskims_sal_vac")
.expect("open real plant index");
let layer_dirs = super::family_scan::sibling_layer_dirs(&idx).expect("layer dirs");
// A `layer_1` (the larger, merged layer) — pick the one already
// profiled: part_00018/index/layer_1, ~30.2M rows.
let layer_dir = layer_dirs.iter()
.find(|p| p.to_string_lossy().contains("part_00018") && p.to_string_lossy().contains("layer_1"))
.expect("expected layer not found — layer_dirs order may have changed");
let mat = obicompactvec::PersistentBitMatrix::open(layer_dir).expect("open presence matrix");
let n = mat.n();
let n_cols = mat.n_cols();
println!("layer={} n_slots={n} n_cols={n_cols}", layer_dir.display());
let rss_before = obisys::peak_rss_bytes();
let t0 = Instant::now();
// Plan design: is_multi flag (1 bit/row) + two SEPARATE arrays
// (singleton-only, multi-only), each sized to its own value range —
// not one array covering the whole dict_id space. In-memory prototype
// (u32 per entry, not yet bit-packed — the fixed-width primitive
// isn't built yet); this benchmark measures RAM + the split's real
// final size, not the construction-time representation's size.
let mut is_multi: Vec<bool> = Vec::with_capacity(n);
let mut singleton_array: Vec<u32> = Vec::new();
let mut multi_array: Vec<u32> = Vec::new();
let mut dedup: HashMap<Vec<u32>, u32> = HashMap::new();
let mut next_dict_id: u32 = 0;
let mut dict_values_bytes: u64 = 0; // running varint-encoded size estimate
let mut buf = vec![0u32; n_cols];
for slot in 0..n {
mat.fill_row(slot, &mut buf);
let set: Vec<u32> = (0..n_cols).filter(|&c| buf[c] != 0).map(|c| c as u32).collect();
match set.len() {
0 => { is_multi.push(false); singleton_array.push(0); } // shouldn't happen on a real built index; placeholder
1 => {
is_multi.push(false);
singleton_array.push(set[0]);
}
_ => {
is_multi.push(true);
let id = if let Some(&id) = dedup.get(&set) {
id
} else {
let id = next_dict_id;
next_dict_id += 1;
// varint size estimate: 1 byte per index < 128 (always
// true here, n_cols=91), matching the plan's encoding.
dict_values_bytes += set.iter().map(|&v| if v < 128 { 1 } else { 2 }).sum::<u64>();
dedup.insert(set, id);
id
};
multi_array.push(id);
}
};
}
let elapsed = t0.elapsed();
let rss_after = obisys::peak_rss_bytes();
let n_distinct_multi = dedup.len() as u64;
let n_singleton = singleton_array.len() as u64;
let n_multi = multi_array.len() as u64;
let is_multi_bytes = (n as u64).div_ceil(8);
let singleton_width_bits = (32 - (n_cols as u32).max(1).leading_zeros()).max(1);
let multi_width_bits = (32 - (n_distinct_multi as u32).max(1).leading_zeros()).max(1);
let singleton_array_bytes = (n_singleton * singleton_width_bits as u64).div_ceil(8);
let multi_array_bytes = (n_multi * multi_width_bits as u64).div_ceil(8);
let split_total = is_multi_bytes + singleton_array_bytes + multi_array_bytes + dict_values_bytes;
let dense_bytes = (n as u64 * n_cols as u64).div_ceil(8);
println!(
"elapsed={:?} rss_before={} rss_after={} rss_delta={}",
elapsed, fmt_mb(rss_before), fmt_mb(rss_after), fmt_mb(rss_after.saturating_sub(rss_before)),
);
println!(
"n_singleton={n_singleton} n_multi={n_multi} n_distinct_multi={n_distinct_multi} \
singleton_width_bits={singleton_width_bits} multi_width_bits={multi_width_bits}",
);
println!(
"is_multi_bytes={} singleton_array_bytes={} multi_array_bytes={} dict_values_bytes={} \
split_total_bytes={} ({:.1}x vs dense) dense_bytes={}",
fmt_mb(is_multi_bytes), fmt_mb(singleton_array_bytes), fmt_mb(multi_array_bytes),
fmt_mb(dict_values_bytes), fmt_mb(split_total),
dense_bytes as f64 / split_total as f64,
fmt_mb(dense_bytes),
);
}
fn fmt_mb(bytes: u64) -> String {
format!("{:.1}MB", bytes as f64 / (1024.0 * 1024.0))
}
#[test]
#[ignore]
fn bench_real_persistent_sparse_bit_matrix_on_disk_size() {
use std::time::Instant;
let idx = KmerIndex::open("/Users/coissac/travail/obiskim/data/phyloalps/phyloskims_sal_vac")
.expect("open real plant index");
let layer_dirs = super::family_scan::sibling_layer_dirs(&idx).expect("layer dirs");
let layer_dir = layer_dirs.iter()
.find(|p| p.to_string_lossy().contains("part_00018") && p.to_string_lossy().contains("layer_1"))
.expect("expected layer not found — layer_dirs order may have changed");
let dense = obicompactvec::PersistentBitMatrix::open(layer_dir).expect("open presence matrix");
let n = dense.n();
let n_cols = dense.n_cols();
println!("layer={} n_slots={n} n_cols={n_cols}", layer_dir.display());
let out_dir = std::env::temp_dir().join(format!("sparse_bench_{}", std::process::id()));
let _ = std::fs::remove_dir_all(&out_dir);
let rss_before = obisys::peak_rss_bytes();
let t0 = Instant::now();
let sparse = obicompactvec::PersistentSparseBitMatrixBuilder::build_from_dense(&dense, &out_dir)
.expect("build_from_dense")
.finish()
.expect("finish");
let elapsed = t0.elapsed();
let rss_after = obisys::peak_rss_bytes();
// Real on-disk size: sum of every file actually written under out_dir.
let mut sparse_bytes: u64 = 0;
for entry in std::fs::read_dir(&out_dir).unwrap() {
let entry = entry.unwrap();
let size = entry.metadata().unwrap().len();
println!(" {}: {}", entry.file_name().to_string_lossy(), fmt_mb(size));
sparse_bytes += size;
}
let dense_bytes = (n as u64 * n_cols as u64).div_ceil(8);
println!(
"elapsed={:?} rss_before={} rss_after={} rss_delta={}",
elapsed, fmt_mb(rss_before), fmt_mb(rss_after), fmt_mb(rss_after.saturating_sub(rss_before)),
);
println!(
"REAL on-disk: sparse_total={} dense={} ({:.1}x)",
fmt_mb(sparse_bytes), fmt_mb(dense_bytes), dense_bytes as f64 / sparse_bytes as f64,
);
// Sanity: reopen from disk (fresh mmap, not the just-built handle) and
// spot-check a handful of rows against the dense original.
drop(sparse);
let reopened = obicompactvec::PersistentSparseBitMatrix::open(&out_dir).expect("reopen");
let mut dense_buf = vec![0u32; n_cols];
for &slot in &[0usize, 1, 1000, n / 2, n - 1] {
dense.fill_row(slot, &mut dense_buf);
let sparse_row = reopened.row(slot);
for c in 0..n_cols {
assert_eq!(sparse_row[c], dense_buf[c] != 0, "slot {slot}, col {c}");
}
}
println!("spot-check rows: OK");
// ── Access-time comparison ──────────────────────────────────────────
// Both patterns matter in practice: sequential (a full-layer scan, the
// existing `scan_layer_families` shape) and random (a single-family
// cross-partition lookup, the entropy/`--shannon` shape). Same slot
// sequence used against both structures for a fair comparison.
const N_ACCESS: usize = 2_000_000;
// Deterministic xorshift, no extra dependency — fine for a benchmark's
// access pattern, not for anything security- or correctness-sensitive.
let mut rng_state: u64 = 0x9E3779B97F4A7C15;
let mut next_rand = move || {
rng_state ^= rng_state << 13;
rng_state ^= rng_state >> 7;
rng_state ^= rng_state << 17;
rng_state
};
let random_slots: Vec<usize> = (0..N_ACCESS).map(|_| (next_rand() as usize) % n).collect();
let sequential_slots: Vec<usize> = (0..N_ACCESS).map(|i| i % n).collect();
let mut buf = vec![0u32; n_cols];
let t = Instant::now();
for &slot in &sequential_slots {
dense.fill_row(slot, &mut buf);
}
let dense_seq = t.elapsed();
let t = Instant::now();
for &slot in &sequential_slots {
reopened.fill_row(slot, &mut buf);
}
let sparse_seq = t.elapsed();
let t = Instant::now();
for &slot in &random_slots {
dense.fill_row(slot, &mut buf);
}
let dense_rand = t.elapsed();
let t = Instant::now();
for &slot in &random_slots {
reopened.fill_row(slot, &mut buf);
}
let sparse_rand = t.elapsed();
println!(
"ACCESS ({N_ACCESS} rows) sequential: dense={:?} ({:.0}ns/row) sparse={:?} ({:.0}ns/row) — {:.2}x",
dense_seq, dense_seq.as_nanos() as f64 / N_ACCESS as f64,
sparse_seq, sparse_seq.as_nanos() as f64 / N_ACCESS as f64,
sparse_seq.as_secs_f64() / dense_seq.as_secs_f64(),
);
println!(
"ACCESS ({N_ACCESS} rows) random: dense={:?} ({:.0}ns/row) sparse={:?} ({:.0}ns/row) — {:.2}x",
dense_rand, dense_rand.as_nanos() as f64 / N_ACCESS as f64,
sparse_rand, sparse_rand.as_nanos() as f64 / N_ACCESS as f64,
sparse_rand.as_secs_f64() / dense_rand.as_secs_f64(),
);
// ── Column-major access, "just for fun" ─────────────────────────────
// The whole point of `docmd/architecture/siblings.md`'s "Explicitly
// deferred" section: dense is genome-major (native, contiguous column
// access), sparse is k-mer-major (no column method at all — reading
// one column means decoding every row and keeping one bit each time).
// Extract one full column (all n rows) both ways.
let col = n_cols / 2;
let t = Instant::now();
let dense_col_view = dense.col_view(col);
let dense_col_ones: u64 = (0..n).filter(|&s| dense_col_view.get(s)).count() as u64;
let dense_col_time = t.elapsed();
let t = Instant::now();
let mut buf2 = vec![0u32; n_cols];
let sparse_col_ones: u64 = (0..n)
.filter(|&s| { reopened.fill_row(s, &mut buf2); buf2[col] != 0 })
.count() as u64;
let sparse_col_time = t.elapsed();
assert_eq!(dense_col_ones, sparse_col_ones, "column {col} popcount disagrees");
println!(
"COLUMN-MAJOR (col {col}, {n} rows) dense={:?} ({:.0}ns/row) sparse={:?} ({:.0}ns/row) — {:.1}x SLOWER on sparse",
dense_col_time, dense_col_time.as_nanos() as f64 / n as f64,
sparse_col_time, sparse_col_time.as_nanos() as f64 / n as f64,
sparse_col_time.as_secs_f64() / dense_col_time.as_secs_f64(),
);
let _ = std::fs::remove_dir_all(&out_dir);
}