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<u8>,
+    /// 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,
 }
 

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<usize> {
+        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<u64> {
+        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<u64> {
+        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<u64> {
+        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 ─────────────────────────────────────────────────────────

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 ───────────────────────────────────────

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<u64> {
+    values.iter().copied().find(|&v| v >= x)
+}
+
+fn naive_successor(values: &[u64], x: u64) -> Option<u64> {
+    values.iter().copied().find(|&v| v > x)
+}
+
+fn naive_predecessor(values: &[u64], x: u64) -> Option<u64> {
+    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<u64> = (0..2000).map(|i| i * 5 + (i % 3)).collect();

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<bool> = (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.