diff --git a/src/obicompactvec/src/bitmatrix/sparse.rs b/src/obicompactvec/src/bitmatrix/sparse.rs index c1f1b806..156071e5 100644 --- a/src/obicompactvec/src/bitmatrix/sparse.rs +++ b/src/obicompactvec/src/bitmatrix/sparse.rs @@ -111,14 +111,42 @@ fn read_varint(data: &[u8], pos: &mut usize) -> u32 { // ── PersistentSparseBitMatrix ─────────────────────────────────────────────── +/// 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, +/// rather than an `n_cols`-wide bit row. +/// +/// `is_multi[slot]` is the row-kind flag; it also acts as the split point +/// between `singleton` and `multi` — `rank0(slot)`/`rank1(slot)` convert a +/// slot into its position within whichever of the two arrays holds it +/// (see [`for_each_genome_in_row`](Self::for_each_genome_in_row)): +/// - `is_multi[slot] == false`: `singleton[rank0(slot)]` is the one genome +/// index present at that slot. +/// - `is_multi[slot] == true`: `multi[rank1(slot)]` is a `dict_id` into the +/// deduplicated multi-genome dictionary — many slots sharing the same +/// genome set (common in practice) point at the same entry instead of +/// repeating it. +/// +/// The dictionary itself is `dict_offsets`/`dict_values`: each distinct +/// multi-genome set is varint-encoded (ascending genome indices, delta +/// isn't taken — see [`write_varint`]/[`read_varint`]) as a byte run in +/// `dict_values`, and `dict_offsets` (an [`EliasFano`] over the +/// `n_distinct_multi` sets' *byte* offsets — monotone by construction, one +/// entry per `dict_id`) gives `dict_entry_range(dict_id)`'s `[start, end)` +/// window in +/// `dict_values`. Looked up only by known `dict_id`, never searched by +/// value. pub struct PersistentSparseBitMatrix { is_multi: PersistentRankSelectBitVec, singleton: PersistentFixedIntVec, multi: PersistentFixedIntVec, dict_offsets: EliasFano, dict_values: Vec, + /// Row count (k-mer slots). n: usize, + /// Column count (genomes). n_cols: usize, + /// Number of distinct multi-genome sets in the dictionary. n_distinct_multi: usize, } diff --git a/src/obicompactvec/src/eliasfano.rs b/src/obicompactvec/src/eliasfano.rs index 9449fc20..52758b83 100644 --- a/src/obicompactvec/src/eliasfano.rs +++ b/src/obicompactvec/src/eliasfano.rs @@ -76,6 +76,73 @@ impl EliasFano { let high_part = pos as u64 - i as u64; (high_part << self.low_width) | self.low.get(i) } + + /// Index of the smallest encoded value `>= x`, or `None` if every value + /// is `< x`. Shared by [`next_geq`](Self::next_geq) and + /// [`predecessor`](Self::predecessor) — the latter needs the index, + /// not just the value, to step one position back. + /// + /// The high bitvector is the unary code `1^c0 0 1^c1 0 ... 1^cH 0` + /// (`c_h` = count of values whose high part is `h`) — each `0` marks a + /// bucket boundary. `select0(h - 1) + 1` jumps straight to the start of + /// bucket `h`; `rank1` on that position converts it to the value's + /// index. From there the scan is bounded by bucket `h`'s local density, + /// not by `n`: once the high part exceeds `h` the value already + /// exceeds `x` (`x < (h + 1) << low_width` by construction of `h`). + fn next_geq_index(&self, x: u64) -> Option { + if self.n == 0 { + return None; + } + let h = x >> self.low_width; + let start = if h == 0 { + 0 + } else { + let zeros_before_bucket = h - 1; + if zeros_before_bucket >= self.high.count_zeros() { + return None; // bucket h is beyond the last one present + } + self.high.select0(zeros_before_bucket) + 1 + }; + let mut i = self.high.rank1(start) as usize; + while i < self.n { + if self.get(i) >= x { + return Some(i); + } + i += 1; + } + None + } + + /// Smallest encoded value `>= x`, or `None` if every value is `< x`. + pub fn next_geq(&self, x: u64) -> Option { + self.next_geq_index(x).map(|i| self.get(i)) + } + + /// Whether `x` is present in the sequence. + pub fn contains(&self, x: u64) -> bool { + self.next_geq(x) == Some(x) + } + + /// Smallest encoded value `> x`, or `None` if every value is `<= x`. + pub fn successor(&self, x: u64) -> Option { + let x_plus_1 = x.checked_add(1)?; // x == u64::MAX: nothing can exceed it + self.next_geq(x_plus_1) + } + + /// Largest encoded value `<= x`, or `None` if every value is `> x`. + pub fn predecessor(&self, x: u64) -> Option { + let idx = self.next_geq_index(x).unwrap_or(self.n); + // `idx` is the first index with value `>= x`. If it lands exactly + // on `x`, that's the answer (largest value `<= x` is `x` itself); + // otherwise step back one to the last value `< x`. + if idx < self.n && self.get(idx) == x { + Some(x) + } else if idx == 0 { + None + } else { + Some(self.get(idx - 1)) + } + } } // ── EliasFanoBuilder ───────────────────────────────────────────────────────── diff --git a/src/obicompactvec/src/rankselect.rs b/src/obicompactvec/src/rankselect.rs index a186ab10..704fca54 100644 --- a/src/obicompactvec/src/rankselect.rs +++ b/src/obicompactvec/src/rankselect.rs @@ -1,11 +1,12 @@ //! 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 +//! in a prefix) and `select1`/`select0` (position of the k-th 1/0 bit). +//! Needed by two consumers in the sparse presence-matrix design (see //! `DevDocMD/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). +//! bits (needs select1 to decode a value by index, select0 to jump to a +//! bucket boundary for `next_geq`/`contains`). //! //! Two-level structure, the standard succinct-bitvector approach: a //! cumulative rank sampled every [`BLOCK_WORDS`] words, plus a linear scan @@ -161,6 +162,43 @@ impl PersistentRankSelectBitVec { } unreachable!("select1({k}): ran off the end of its own block — rank/select index inconsistent"); } + + /// Position of the `k`-th (0-indexed) 0 bit. Panics if fewer than + /// `k + 1` zeros exist. + /// + /// Mirrors [`select1`](Self::select1): binary-searches a per-block + /// zero count derived from `block_ranks` (`block_start_bit - + /// block_ranks[block]`, no separate index needed), then locates the + /// exact bit within the winning word via `SelectInWord` on `!w`. + pub fn select0(&self, k: u64) -> usize { + let total_zeros = self.count_zeros(); + assert!(k < total_zeros, "select0({k}) out of range: only {total_zeros} zeros"); + let block_ranks = self.block_ranks(); + let words = self.bit_words(); + + // Binary search: last block whose cumulative zero-count is <= k. + let mut lo = 0usize; + let mut hi = block_ranks.len(); + while lo + 1 < hi { + let mid = lo + (hi - lo) / 2; + let zeros_at_mid = (mid * BLOCK_WORDS * 64) as u64 - block_ranks[mid]; + if zeros_at_mid <= k { lo = mid; } else { hi = mid; } + } + let block_idx = lo; + let zeros_before_block = (block_idx * BLOCK_WORDS * 64) as u64 - block_ranks[block_idx]; + let mut remaining = k - zeros_before_block; + 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 zw = !w; + let c = zw.count_ones() as u64; + if remaining < c { + return (start_word + i) * 64 + zw.select_in_word(remaining as usize); + } + remaining -= c; + } + unreachable!("select0({k}): ran off the end of its own block — rank/select index inconsistent"); + } } // ── PersistentRankSelectBitVecBuilder ─────────────────────────────────────── diff --git a/src/obicompactvec/src/tests/eliasfano.rs b/src/obicompactvec/src/tests/eliasfano.rs index 9254e1ab..891cce7e 100644 --- a/src/obicompactvec/src/tests/eliasfano.rs +++ b/src/obicompactvec/src/tests/eliasfano.rs @@ -94,6 +94,106 @@ fn push_out_of_order_panics() { b.push(5); // decreasing — must panic } +fn naive_next_geq(values: &[u64], x: u64) -> Option { + values.iter().copied().find(|&v| v >= x) +} + +fn naive_successor(values: &[u64], x: u64) -> Option { + values.iter().copied().find(|&v| v > x) +} + +fn naive_predecessor(values: &[u64], x: u64) -> Option { + values.iter().copied().filter(|&v| v <= x).next_back() +} + +#[test] +fn next_geq_and_contains_small() { + let values = [0u64, 3, 3, 7, 42, 42, 42, 100]; + let (_dir, ef) = build(&values, 1000); + for x in 0..=105u64 { + assert_eq!(ef.next_geq(x), naive_next_geq(&values, x), "next_geq({x})"); + assert_eq!(ef.contains(x), values.contains(&x), "contains({x})"); + } +} + +#[test] +fn next_geq_empty_bucket_skips_to_next() { + // Bucket for values in [8, 16) is empty — next_geq(9) must skip it. + let values = [0u64, 1, 20, 21]; + let (_dir, ef) = build(&values, 100); + assert_eq!(ef.next_geq(9), Some(20)); + assert_eq!(ef.next_geq(21), Some(21)); + assert_eq!(ef.next_geq(22), None); + assert!(!ef.contains(9)); + assert!(ef.contains(21)); +} + +#[test] +fn next_geq_empty_sequence() { + let (_dir, ef) = build(&[], 1000); + assert_eq!(ef.next_geq(0), None); + assert!(!ef.contains(0)); +} + +#[test] +fn next_geq_beyond_last_value() { + let values = [5u64, 10, 15]; + let (_dir, ef) = build(&values, 100); + assert_eq!(ef.next_geq(16), None); + assert_eq!(ef.next_geq(999), None); +} + +#[test] +fn successor_and_predecessor_small() { + let values = [0u64, 3, 3, 7, 42, 42, 42, 100]; + let (_dir, ef) = build(&values, 1000); + for x in 0..=105u64 { + assert_eq!(ef.successor(x), naive_successor(&values, x), "successor({x})"); + assert_eq!(ef.predecessor(x), naive_predecessor(&values, x), "predecessor({x})"); + } +} + +#[test] +fn successor_at_u64_max_is_none() { + let values = [1u64, 2, 3]; + let (_dir, ef) = build(&values, 100); + assert_eq!(ef.successor(u64::MAX), None); +} + +#[test] +fn predecessor_before_first_value_is_none() { + let values = [5u64, 10, 15]; + let (_dir, ef) = build(&values, 100); + assert_eq!(ef.predecessor(4), None); + assert_eq!(ef.predecessor(5), Some(5)); +} + +#[test] +fn successor_predecessor_empty_sequence() { + let (_dir, ef) = build(&[], 1000); + assert_eq!(ef.successor(0), None); + assert_eq!(ef.predecessor(0), None); +} + +#[test] +fn next_geq_matches_naive_realistic_offsets() { + let mut offset = 0u64; + let mut values = Vec::new(); + let mut rng_state = 98765u64; + for _ in 0..2000 { + values.push(offset); + rng_state ^= rng_state << 13; + rng_state ^= rng_state >> 7; + rng_state ^= rng_state << 17; + offset += 1 + (rng_state % 8); + } + let universe = offset + 1; + let (_dir, ef) = build(&values, universe); + for x in (0..universe).step_by(7) { + assert_eq!(ef.next_geq(x), naive_next_geq(&values, x), "next_geq({x})"); + } +} + #[test] fn reopen_after_close_matches_original() { let values: Vec = (0..2000).map(|i| i * 5 + (i % 3)).collect(); diff --git a/src/obicompactvec/src/tests/rankselect.rs b/src/obicompactvec/src/tests/rankselect.rs index 92e3e492..6d020ac8 100644 --- a/src/obicompactvec/src/tests/rankselect.rs +++ b/src/obicompactvec/src/tests/rankselect.rs @@ -22,6 +22,10 @@ fn naive_select1(bits: &[bool], k: u64) -> usize { bits.iter().enumerate().filter(|&(_, &b)| b).nth(k as usize).unwrap().0 } +fn naive_select0(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]; @@ -91,6 +95,35 @@ fn select1_matches_naive_multi_block() { } } +#[test] +fn select0_matches_naive_small() { + let bits = [true, false, true, true, false, false, true, true, false, true]; + let (_dir, r) = build(&bits); + let n_zeros = bits.iter().filter(|&&b| !b).count() as u64; + for k in 0..n_zeros { + assert_eq!(r.select0(k), naive_select0(&bits, k), "select0({k})"); + } +} + +#[test] +fn select0_matches_naive_multi_block() { + let n = 3000; + let bits: Vec = (0..n).map(|i| (i * 13 + 5) % 17 == 0).collect(); + let (_dir, r) = build(&bits); + let n_zeros = bits.iter().filter(|&&b| !b).count() as u64; + for k in (0..n_zeros).step_by(23) { + assert_eq!(r.select0(k), naive_select0(&bits, k), "select0({k})"); + } +} + +#[test] +#[should_panic] +fn select0_out_of_range_panics() { + let bits = [true, true, false]; + let (_dir, r) = build(&bits); + r.select0(1); // only one 0-bit (k=0 valid), k=1 must panic +} + #[test] fn rank_select_round_trip() { // For every one-bit's position p, select1(rank1(p)) == p.