Move sibling annex implementation to obikphylo siblings module

Relocates `SiblingAnnex`, `FamilyMask`, and `SiblingAnnexBuilder` from `obicompactvec` to the local `obikphylo::siblings` module. Replaces the previous implementation with a memory-mapped version using `memmap2`, featuring a 2-byte-per-slot layout, explicit bitfield manipulation, and support for concurrent atomic writes. Updates all sibling module imports to local paths and adds the `memmap2` dependency to `obikphylo`.
This commit is contained in:
Eric Coissac
2026-08-16 21:04:14 +02:00
parent c990087ef3
commit fecfe84ea6
13 changed files with 404 additions and 293 deletions
+1
View File
@@ -1846,6 +1846,7 @@ dependencies = [
name = "obikphylo" name = "obikphylo"
version = "0.1.0" version = "0.1.0"
dependencies = [ dependencies = [
"memmap2",
"ndarray", "ndarray",
"obicompactvec", "obicompactvec",
"obikindex", "obikindex",
-2
View File
@@ -7,7 +7,6 @@ mod intmatrix;
mod layer_meta; mod layer_meta;
mod meta; mod meta;
mod reader; mod reader;
mod siblingannex;
mod tempbitvec; mod tempbitvec;
mod tempintvec; mod tempintvec;
mod views; mod views;
@@ -19,7 +18,6 @@ pub use builder::PersistentCompactIntVecBuilder;
pub use colgroup::{ColGroup, FilterMask, MatrixGroupOps, eval_filter_mask}; pub use colgroup::{ColGroup, FilterMask, MatrixGroupOps, eval_filter_mask};
pub use intmatrix::{PersistentCompactIntMatrix, PersistentCompactIntMatrixBuilder, pack_compact_int_matrix}; pub use intmatrix::{PersistentCompactIntMatrix, PersistentCompactIntMatrixBuilder, pack_compact_int_matrix};
pub use layer_meta::LayerMeta; pub use layer_meta::LayerMeta;
pub use siblingannex::{FamilyMask, SiblingAnnex, SiblingAnnexBuilder};
pub use reader::{PersistentCompactIntVec, Iter as CompactIntVecIter}; pub use reader::{PersistentCompactIntVec, Iter as CompactIntVecIter};
pub use tempbitvec::{TempBitVec, TempBitVecBuilder}; pub use tempbitvec::{TempBitVec, TempBitVecBuilder};
pub use tempintvec::{TempCompactIntVec, TempCompactIntVecBuilder}; pub use tempintvec::{TempCompactIntVec, TempCompactIntVecBuilder};
-280
View File
@@ -1,280 +0,0 @@
//! Family presence-mask annex: a compact, read-only-after-build, per-slot
//! derived value used by the central-position SNP distance estimator (see
//! `docmd/theory/evolutionary_distances.md`, "Step 2b" and "Definitions:
//! family, and the canonical form of a family").
//!
//! One byte is stored per MPHF slot of a partition/layer, its low 4 bits
//! encoding a **presence mask** for the slot's k-mer's "family" (the up to 4
//! k-mers sharing the same flanks, differing only at the central base):
//! bit `b` (`b` = 0..3, in the fixed A/C/G/T = 0/1/2/3 encoding already used
//! for a single nucleotide) is set iff the family member whose *own* central
//! base — in its own canonical orientation — is `b`, is observed anywhere in
//! the current multi-genome index. This is a property of the whole index,
//! not of any one genome.
//!
//! Both facts the earlier (superseded) 3-bit design stored explicitly are
//! derived from the mask instead, not stored:
//! - sibling count = `popcount(mask) - 1`;
//! - minorant = regenerate the family's 4 canonical forms from the slot's
//! own k-mer (`CanonicalKmerOf::central_canonical_neighbors`, cheap, no
//! lookup), compare the raw encodings of whichever are set in the mask,
//! take the smallest — see `obikindex::siblings`.
//!
//! Mask value 0 is logically unreachable as a real result (a slot's own base
//! is always present in its own family) and is reused as the "not yet
//! computed" sentinel: annex files are pre-initialised to all-zero, and a
//! real value is only ever written once, by the computation pass.
//!
//! Deliberately simpler than a true 4-bit pack (1 byte/slot instead of 4
//! bits/slot): correctness and simplicity first, for a first implementation.
//! Packing to 4 bits/slot is a pure storage-density follow-up, not a
//! behavioural change, left for later.
use std::fs::{File, OpenOptions};
use std::io;
use std::path::{Path, PathBuf};
use memmap2::{Mmap, MmapMut};
const MAGIC: [u8; 4] = *b"PSIB";
// Header: magic(4) + _pad(4) + n(8) = 16 bytes. Data (1 byte/slot) follows.
const HEADER_SIZE: usize = 16;
/// A family presence mask: bit `b` (0..3) set iff the member whose own
/// canonical central base is `b` (0=A, 1=C, 2=G, 3=T) is observed in the
/// index. Bit 4 is a second, independent fact piggy-backed onto the same
/// byte: whether *this slot's own k-mer* is its family's minorant (the
/// smallest raw encoding among the family's observed members) — computed
/// once, when the whole family's mask is already final (see
/// `obikindex::siblings::build_sibling_annex`), and read back by every
/// consumer that would otherwise have to reconstruct this slot's k-mer from
/// `unitigs.bin` and hash it through the MPHF again just to ask the same
/// question. Bits 5-7 unused.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub struct FamilyMask(u8);
const PRESENCE_BITS: u8 = 0b0000_1111;
const MINORANT_BIT: u8 = 0b0001_0000;
impl FamilyMask {
/// The empty mask — never a valid *computed* result (a slot's own base
/// is always present in its own family) — used only to build up a mask
/// via repeated [`with`](Self::with) calls before storing it.
pub const EMPTY: FamilyMask = FamilyMask(0);
/// Set bit `base` (0=A, 1=C, 2=G, 3=T).
#[inline]
pub fn with(self, base: u8) -> Self {
debug_assert!(base < 4, "base out of range: {base}");
FamilyMask(self.0 | (1 << base))
}
/// Is the member with central base `base` (0..3) present?
#[inline]
pub fn has(self, base: u8) -> bool {
debug_assert!(base < 4, "base out of range: {base}");
self.0 & (1 << base) != 0
}
/// Number of family members observed anywhere in the index (1..=4).
#[inline]
pub fn family_size(self) -> u32 {
(self.0 & PRESENCE_BITS).count_ones()
}
/// Number of *other* members observed (0..=3) — `family_size() - 1`.
#[inline]
pub fn siblings(self) -> u32 {
self.family_size() - 1
}
/// Set or clear the minorant flag (see the struct docs).
#[inline]
pub fn with_minorant(self, is_minorant: bool) -> Self {
if is_minorant {
FamilyMask(self.0 | MINORANT_BIT)
} else {
FamilyMask(self.0 & !MINORANT_BIT)
}
}
/// Is this slot's own k-mer its family's minorant? Only meaningful once
/// `with_minorant` has been called with the family's *final* mask (i.e.
/// after construction) — see the struct docs.
#[inline]
pub fn is_minorant(self) -> bool {
self.0 & MINORANT_BIT != 0
}
/// Raw presence bitmask (bit `b` = base `b` present, low 4 bits only —
/// never includes the minorant flag) — for callers that build up a mask
/// via their own bit operations (e.g. concurrently, via an `AtomicU8`)
/// and only need the `FamilyMask` wrapper at the end.
#[inline]
pub fn bits(self) -> u8 {
self.0 & PRESENCE_BITS
}
/// Construct from a raw presence bitmask (only the low 4 bits are kept —
/// the minorant flag is not part of this, use [`with_minorant`](Self::with_minorant)
/// separately).
#[inline]
pub fn from_bits(bits: u8) -> Self {
FamilyMask(bits & PRESENCE_BITS)
}
#[inline]
fn encode(self) -> u8 {
self.0
}
#[inline]
fn decode(byte: u8) -> Option<Self> {
if byte == 0 {
// Unreachable for a real result — reserved as the "not yet
// computed" sentinel. Still safe as a sentinel with the
// minorant bit added: a real entry always has at least one
// presence bit set (bits 0-3), so byte == 0 still means
// "nothing written yet" and never a genuine minorant-only value.
return None;
}
Some(FamilyMask(byte))
}
}
// ── SiblingAnnex (reader) ───────────────────────────────────────────────────
pub struct SiblingAnnex {
mmap: Mmap,
n: usize,
path: PathBuf,
}
impl SiblingAnnex {
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, "PSIB file too short"));
}
if mmap[0..4] != MAGIC {
return Err(io::Error::new(io::ErrorKind::InvalidData, "bad PSIB magic"));
}
let n = u64::from_le_bytes(mmap[8..16].try_into().unwrap()) as usize;
if mmap.len() < HEADER_SIZE + n {
return Err(io::Error::new(io::ErrorKind::InvalidData, "PSIB file truncated"));
}
Ok(Self { mmap, n, path: path.to_path_buf() })
}
pub fn path(&self) -> &Path { &self.path }
pub fn len(&self) -> usize { self.n }
pub fn is_empty(&self) -> bool { self.n == 0 }
/// `None` means the slot has not (yet) been computed — see module docs.
pub fn get(&self, slot: usize) -> Option<FamilyMask> {
FamilyMask::decode(self.mmap[HEADER_SIZE + slot])
}
}
// ── SiblingAnnexBuilder (writer) ────────────────────────────────────────────
pub struct SiblingAnnexBuilder {
mmap: MmapMut,
n: usize,
path: PathBuf,
}
impl SiblingAnnexBuilder {
/// Create a new annex of `n` slots at `path`, pre-initialised to the
/// "not yet computed" sentinel (all-zero).
pub fn new(n: usize, path: &Path) -> io::Result<Self> {
let file_size = HEADER_SIZE + n;
let file = OpenOptions::new()
.read(true).write(true).create(true).truncate(true)
.open(path)?;
file.set_len(file_size as u64)?;
let mut mmap = unsafe { MmapMut::map_mut(&file)? };
mmap[0..4].copy_from_slice(&MAGIC);
mmap[4..8].copy_from_slice(&[0u8; 4]);
mmap[8..16].copy_from_slice(&(n as u64).to_le_bytes());
// Data region left at 0 by `set_len`/mmap — the sentinel value.
Ok(Self { mmap, n, path: path.to_path_buf() })
}
pub fn len(&self) -> usize { self.n }
pub fn is_empty(&self) -> bool { self.n == 0 }
pub fn get(&self, slot: usize) -> Option<FamilyMask> {
FamilyMask::decode(self.mmap[HEADER_SIZE + slot])
}
pub fn set(&mut self, slot: usize, mask: FamilyMask) {
// Redundant concurrent writes from independent recomputation paths
// converge to the same encoded byte for a given slot, so a plain
// store here is safe even without external synchronisation, as long
// as the byte write itself is atomic (true for a single aligned
// byte on every platform this project targets).
self.mmap[HEADER_SIZE + slot] = mask.encode();
}
pub fn close(self) -> io::Result<()> { self.mmap.flush() }
pub fn finish(self) -> io::Result<SiblingAnnex> {
let path = self.path.clone();
self.close()?;
SiblingAnnex::open(&path)
}
}
#[cfg(test)]
mod tests {
use super::*;
use tempfile::tempdir;
#[test]
fn sentinel_is_zero_and_unset_slots_read_as_uncomputed() {
let dir = tempdir().unwrap();
let path = dir.path().join("test.psib");
let builder = SiblingAnnexBuilder::new(4, &path).unwrap();
for slot in 0..4 {
assert_eq!(builder.get(slot), None);
}
builder.close().unwrap();
}
#[test]
fn roundtrip_all_valid_masks() {
let dir = tempdir().unwrap();
let path = dir.path().join("test.psib");
let mut builder = SiblingAnnexBuilder::new(4, &path).unwrap();
let masks = [
FamilyMask::EMPTY.with(0), // just A: family size 1
FamilyMask::EMPTY.with(0).with(3), // A + T: size 2
FamilyMask::EMPTY.with(1).with(2).with(3), // C+G+T: size 3
FamilyMask::EMPTY.with(0).with(1).with(2).with(3), // all 4
];
for (slot, mask) in masks.iter().enumerate() {
builder.set(slot, *mask);
}
let annex = builder.finish().unwrap();
for (slot, mask) in masks.iter().enumerate() {
assert_eq!(annex.get(slot), Some(*mask));
}
assert_eq!(annex.get(0).unwrap().siblings(), 0);
assert_eq!(annex.get(1).unwrap().siblings(), 1);
assert_eq!(annex.get(2).unwrap().siblings(), 2);
assert_eq!(annex.get(3).unwrap().siblings(), 3);
assert_eq!(annex.get(3).unwrap().family_size(), 4);
}
#[test]
fn has_reflects_individual_bits() {
let mask = FamilyMask::EMPTY.with(0).with(2);
assert!(mask.has(0));
assert!(!mask.has(1));
assert!(mask.has(2));
assert!(!mask.has(3));
}
}
+1
View File
@@ -13,6 +13,7 @@ obicompactvec = { path = "../obicompactvec" }
obilayeredmap = { path = "../obilayeredmap" } obilayeredmap = { path = "../obilayeredmap" }
obiskbuilder = { path = "../obiskbuilder" } obiskbuilder = { path = "../obiskbuilder" }
obipipeline = { path = "../obipipeline" } obipipeline = { path = "../obipipeline" }
memmap2 = "0.9"
ndarray = "0.16" ndarray = "0.16"
rayon = "1" rayon = "1"
tracing = "0.1.44" tracing = "0.1.44"
+1 -2
View File
@@ -4,7 +4,6 @@ use std::sync::Arc;
use rayon::prelude::*; use rayon::prelude::*;
use obicompactvec::{FamilyMask, SiblingAnnexBuilder};
use obikpartitionner::KmerPartition; use obikpartitionner::KmerPartition;
use obipipeline::ThrottleGuard; use obipipeline::ThrottleGuard;
use obikseq::CanonicalKmer; use obikseq::CanonicalKmer;
@@ -17,7 +16,7 @@ use obikindex::KmerIndex;
use super::cache::PartitionCache; use super::cache::PartitionCache;
use super::helpers::{central_base, is_minorant}; use super::helpers::{central_base, is_minorant};
use super::{olm_to_ok, ANNEX_FILE_NAME, INDEX_SUBDIR}; use super::{olm_to_ok, FamilyMask, SiblingAnnexBuilder, ANNEX_FILE_NAME, INDEX_SUBDIR};
// ── obipipeline data types ───────────────────────────────────────────────── // ── obipipeline data types ─────────────────────────────────────────────────
+2 -2
View File
@@ -10,7 +10,7 @@ use obisys::progress_bar;
use obikindex::OKIResult; use obikindex::OKIResult;
use super::iter::SiblingLayerExt; use super::iter::SiblingLayerExt;
use super::{olm_to_ok, INDEX_SUBDIR}; use super::{olm_to_ok, SiblingAnnex, INDEX_SUBDIR};
/// Every partition's already-open layers, built **once** for the whole /// Every partition's already-open layers, built **once** for the whole
/// `build_sibling_annex` run and shared (read-only) across every lookup, in /// `build_sibling_annex` run and shared (read-only) across every lookup, in
@@ -68,7 +68,7 @@ impl Mat {
/// doesn't depend on `D`), so no boxing is needed. /// doesn't depend on `D`), so no boxing is needed.
pub(super) fn iter_minorants_batch( pub(super) fn iter_minorants_batch(
&self, &self,
annex: std::sync::Arc<obicompactvec::SiblingAnnex>, annex: std::sync::Arc<SiblingAnnex>,
batch_size: usize, batch_size: usize,
) -> super::iter::MinorantBatchIter { ) -> super::iter::MinorantBatchIter {
match self { match self {
+2 -2
View File
@@ -53,7 +53,7 @@ use std::sync::atomic::{AtomicU8, Ordering};
use rayon::prelude::*; use rayon::prelude::*;
use obicompactvec::{FamilyMask, PersistentBitMatrix, PersistentCompactIntMatrix, SiblingAnnex}; use obicompactvec::{PersistentBitMatrix, PersistentCompactIntMatrix};
use obikseq::CanonicalKmer; use obikseq::CanonicalKmer;
use obilayeredmap::Layer; use obilayeredmap::Layer;
use obilayeredmap::meta::PartitionMeta; use obilayeredmap::meta::PartitionMeta;
@@ -65,7 +65,7 @@ use obikindex::KmerIndex;
use super::cache::{Mat, PartitionCache}; use super::cache::{Mat, PartitionCache};
use super::helpers::central_base; use super::helpers::central_base;
use super::iter::SiblingEntry; use super::iter::SiblingEntry;
use super::{olm_to_ok, ANNEX_FILE_NAME, INDEX_SUBDIR}; use super::{olm_to_ok, FamilyMask, SiblingAnnex, ANNEX_FILE_NAME, INDEX_SUBDIR};
/// Families per batch — see the module docs for the memory-vs-per-partition- /// Families per batch — see the module docs for the memory-vs-per-partition-
/// density trade-off this picks a point on. At ~90 genomes and a few /// density trade-off this picks a point on. At ~90 genomes and a few
+1 -1
View File
@@ -1,6 +1,6 @@
use obikseq::CanonicalKmer; use obikseq::CanonicalKmer;
use obicompactvec::FamilyMask; use super::FamilyMask;
/// Central-position base of a canonical k-mer, in the fixed 0=A/1=C/2=G/3=T /// Central-position base of a canonical k-mer, in the fixed 0=A/1=C/2=G/3=T
/// encoding — the mask's bit index. `k` must be odd (project invariant). /// encoding — the mask's bit index. `k` must be odd (project invariant).
+2 -1
View File
@@ -28,10 +28,11 @@
use std::sync::Arc; use std::sync::Arc;
use obicompactvec::{FamilyMask, SiblingAnnex};
use obikseq::CanonicalKmer; use obikseq::CanonicalKmer;
use obilayeredmap::{KmerIter, Layer, LayerData}; use obilayeredmap::{KmerIter, Layer, LayerData};
use super::{FamilyMask, SiblingAnnex};
/// One layer entry: a k-mer's position in the layer's iteration order (the /// One layer entry: a k-mer's position in the layer's iteration order (the
/// same index the sibling annex is keyed on — not an MPHF slot), the k-mer /// same index the sibling annex is keyed on — not an MPHF slot), the k-mer
/// itself, and its family mask. /// itself, and its family mask.
+2
View File
@@ -50,6 +50,7 @@ mod distance;
mod family_scan; mod family_scan;
mod helpers; mod helpers;
mod iter; mod iter;
mod siblingannex;
mod stats; mod stats;
#[cfg(test)] #[cfg(test)]
@@ -60,6 +61,7 @@ pub use build::SiblingAnnexBuildExt;
pub use cardinality::{CardinalityExt, CardinalityTally}; pub use cardinality::{CardinalityExt, CardinalityTally};
pub use distance::{BasePairTally, DistanceExt, RawSnpDistanceOutput}; pub use distance::{BasePairTally, DistanceExt, RawSnpDistanceOutput};
pub use iter::{MinorantBatchIter, MinorantIter, SiblingBatchIter, SiblingEntry, SiblingIter, SiblingLayerExt}; pub use iter::{MinorantBatchIter, MinorantIter, SiblingBatchIter, SiblingEntry, SiblingIter, SiblingLayerExt};
pub(crate) use siblingannex::{FamilyMask, SiblingAnnex, SiblingAnnexBuilder};
pub use stats::{SiblingAnnexStats, SiblingStatsExt}; pub use stats::{SiblingAnnexStats, SiblingStatsExt};
use obilayeredmap::OLMError; use obilayeredmap::OLMError;
+390
View File
@@ -0,0 +1,390 @@
//! Family presence-mask annex: a compact, read-only-after-build, per-slot
//! derived value used by the central-position SNP distance estimator (see
//! `docmd/theory/evolutionary_distances.md`, "Step 2b" and "Definitions:
//! family, and the canonical form of a family").
//!
//! Two bytes are stored per MPHF slot of a partition/layer, packed as four
//! 3-bit fields (one per central base, in the fixed A/C/G/T = 0/1/2/3
//! encoding) plus a minorant flag:
//!
//! - Each 3-bit field: `0` = the family member whose *own* central base —
//! in its own canonical orientation — is that base is absent from the
//! whole multi-genome index (a property of the whole index, not of any
//! one genome); a non-zero value `v` (`1..=7`) means it's present, and
//! *may* additionally encode `layer = v - 1` (the destination layer
//! within its partition) once a caller populates it via
//! [`with_layer`](FamilyMask::with_layer) — see that method's docs.
//! `FamilyMask` itself is policy-free about whether a given non-zero
//! value is a trustworthy layer index or just "present, no layer
//! recorded": that decision belongs to callers, who already have
//! `PartitionMeta::n_layers` in scope.
//! - The minorant bit: whether *this slot's own k-mer* is its family's
//! minorant (the smallest raw encoding among the family's observed
//! members) — computed once, when the whole family's mask is already
//! final (see `build::build_layer_sibling_annex`), and read back by every
//! consumer that would otherwise have to reconstruct this slot's k-mer
//! from `unitigs.bin` and hash it through the MPHF again just to ask the
//! same question.
//!
//! Both facts an even earlier design stored explicitly are still derived,
//! not stored:
//! - sibling count = `family_size() - 1`;
//! - which of the (up to 4) family members are present, and their canonical
//! form — [`FamilyMask::family_members`], regenerated from the slot's own
//! k-mer (`CanonicalKmer::central_canonical_neighbors`, cheap, no
//! lookup), not stored.
//!
//! Mask value 0 is logically unreachable as a real result (a slot's own
//! base is always present in its own family) and is reused as the "not yet
//! computed" sentinel: annex files are pre-initialised to all-zero, and a
//! real value is only ever written once, by the computation pass.
//!
//! Widened from the original 1-byte/slot design (4 presence bits + 1
//! minorant bit) to carry a per-member layer number without an extra
//! lookup pass — see the project discussion this implements. An old
//! (1-byte/slot) `.psib` file is *not* silently misread by the new reader:
//! its length no longer matches `HEADER_SIZE + n * 2`, so `SiblingAnnex::open`
//! fails loudly ("PSIB file truncated") rather than producing garbage.
//! Annexes built before this change must be rebuilt (`--sibling-annex`).
use std::fs::{File, OpenOptions};
use std::io;
use std::path::{Path, PathBuf};
use memmap2::{Mmap, MmapMut};
use obikseq::CanonicalKmer;
use super::helpers::central_base;
const MAGIC: [u8; 4] = *b"PSIB";
// Header: magic(4) + _pad(4) + n(8) = 16 bytes. Data (2 bytes/slot) follows.
const HEADER_SIZE: usize = 16;
/// Width, in bits, of one base's field.
const FIELD_BITS: u32 = 3;
/// Mask for one base's field once shifted into position.
const FIELD_MASK: u16 = 0b111;
/// Highest raw field value a base can carry (`1..=MAX_FIELD_VALUE`, `0` = absent).
const MAX_FIELD_VALUE: u8 = 7;
const MINORANT_BIT: u16 = 1 << 12;
/// A family presence mask — see the module docs for the full bit layout.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub struct FamilyMask(u16);
impl FamilyMask {
/// The empty mask — never a valid *computed* result (a slot's own base
/// is always present in its own family) — used only to build up a mask
/// via repeated [`with`](Self::with)/[`with_layer`](Self::with_layer)
/// calls before storing it.
pub const EMPTY: FamilyMask = FamilyMask(0);
#[inline]
fn field_shift(base: u8) -> u32 {
debug_assert!(base < 4, "base out of range: {base}");
base as u32 * FIELD_BITS
}
#[inline]
fn field(self, base: u8) -> u8 {
((self.0 >> Self::field_shift(base)) & FIELD_MASK) as u8
}
#[inline]
fn with_field(self, base: u8, value: u8) -> Self {
debug_assert!(value <= MAX_FIELD_VALUE, "field value out of range: {value}");
let shift = Self::field_shift(base);
let cleared = self.0 & !(FIELD_MASK << shift);
FamilyMask(cleared | ((value as u16) << shift))
}
/// Mark the member with central base `base` (0=A, 1=C, 2=G, 3=T) as
/// present, without recording a layer (compat path: same observable
/// effect as the original 1-byte design's `with`). Use
/// [`with_layer`](Self::with_layer) instead when the destination layer
/// is already known.
#[inline]
pub fn with(self, base: u8) -> Self {
self.with_field(base, 1)
}
/// Mark the member with central base `base` as present *and* record its
/// destination layer. `layer` must be `< 7` (`debug_assert`ed) — a
/// caller facing more layers than that has no compact field to record
/// them in and must fall back to [`with`](Self::with) instead (still
/// correct, just without the fast-path payoff); see the module docs on
/// why `FamilyMask` doesn't decide that threshold itself.
#[inline]
pub fn with_layer(self, base: u8, layer: usize) -> Self {
debug_assert!(layer < MAX_FIELD_VALUE as usize, "layer out of range: {layer}");
self.with_field(base, layer as u8 + 1)
}
/// Is the member with central base `base` (0..3) present?
#[inline]
pub fn has(self, base: u8) -> bool {
self.field(base) != 0
}
/// Raw stored value for `base`'s field, if present: `Some(v - 1)` where
/// `v` is the non-zero field value. This is the *raw* stored value, not
/// a validated layer index — a caller must cross-check it against
/// `PartitionMeta::n_layers` (`<= 7`) before trusting it as a real
/// layer, since a mask written via [`with`](Self::with) (no layer
/// known) also reads back as `Some(0)` here.
#[inline]
pub fn layer_value(self, base: u8) -> Option<u8> {
let v = self.field(base);
if v == 0 { None } else { Some(v - 1) }
}
/// Number of family members observed anywhere in the index (1..=4).
#[inline]
pub fn family_size(self) -> u32 {
(0..4).filter(|&b| self.has(b)).count() as u32
}
/// Number of *other* members observed (0..=3) — `family_size() - 1`.
#[inline]
pub fn siblings(self) -> u32 {
self.family_size() - 1
}
/// Set or clear the minorant flag (see the module docs).
#[inline]
pub fn with_minorant(self, is_minorant: bool) -> Self {
if is_minorant {
FamilyMask(self.0 | MINORANT_BIT)
} else {
FamilyMask(self.0 & !MINORANT_BIT)
}
}
/// Is this slot's own k-mer its family's minorant? Only meaningful once
/// `with_minorant` has been called with the family's *final* mask (i.e.
/// after construction) — see the module docs.
#[inline]
pub fn is_minorant(self) -> bool {
self.0 & MINORANT_BIT != 0
}
/// Raw presence bitmask (bit `b` = base `b` present, low 4 bits only —
/// never includes the minorant flag or any layer information) — for
/// callers that only need the old 1-bit-per-base view, e.g. to compare
/// against a previously-computed value built via their own bit
/// operations.
#[inline]
pub fn bits(self) -> u8 {
(0..4).fold(0u8, |acc, b| if self.has(b) { acc | (1 << b) } else { acc })
}
/// Construct from a raw presence bitmask (only the low 4 bits are kept
/// — the minorant flag is not part of this, use
/// [`with_minorant`](Self::with_minorant) separately). Each set bit is
/// recorded as "present, no layer known" — see [`with`](Self::with).
#[inline]
pub fn from_bits(bits: u8) -> Self {
(0..4).fold(FamilyMask::EMPTY, |mask, b| {
if bits & (1 << b) != 0 { mask.with(b) } else { mask }
})
}
/// This family's present members, in canonical form, paired with their
/// raw stored field value (see [`layer_value`](Self::layer_value) for
/// why it's not directly a validated layer index). `owner` is this
/// slot's own k-mer (any member works — the 4 canonical forms are
/// invariant regardless of which member you start from).
pub fn family_members(self, owner: CanonicalKmer, k: usize) -> impl Iterator<Item = (CanonicalKmer, Option<u8>)> {
owner
.central_canonical_neighbors()
.into_iter()
.filter_map(move |member| {
let base = central_base(member, k);
self.has(base).then(|| (member, self.layer_value(base)))
})
}
#[inline]
fn encode(self) -> u16 {
self.0
}
#[inline]
fn decode(word: u16) -> Option<Self> {
if word == 0 {
// Unreachable for a real result — reserved as the "not yet
// computed" sentinel. A real entry always has at least one
// non-zero base field, so word == 0 always means "nothing
// written yet", never a genuine minorant-only value.
return None;
}
Some(FamilyMask(word))
}
}
// ── SiblingAnnex (reader) ───────────────────────────────────────────────────
pub struct SiblingAnnex {
mmap: Mmap,
n: usize,
path: PathBuf,
}
impl SiblingAnnex {
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, "PSIB file too short"));
}
if mmap[0..4] != MAGIC {
return Err(io::Error::new(io::ErrorKind::InvalidData, "bad PSIB magic"));
}
let n = u64::from_le_bytes(mmap[8..16].try_into().unwrap()) as usize;
if mmap.len() < HEADER_SIZE + n * 2 {
return Err(io::Error::new(io::ErrorKind::InvalidData, "PSIB file truncated"));
}
Ok(Self { mmap, n, path: path.to_path_buf() })
}
pub fn path(&self) -> &Path { &self.path }
pub fn len(&self) -> usize { self.n }
pub fn is_empty(&self) -> bool { self.n == 0 }
/// `None` means the slot has not (yet) been computed — see module docs.
pub fn get(&self, slot: usize) -> Option<FamilyMask> {
let off = HEADER_SIZE + slot * 2;
FamilyMask::decode(u16::from_le_bytes(self.mmap[off..off + 2].try_into().unwrap()))
}
}
// ── SiblingAnnexBuilder (writer) ────────────────────────────────────────────
pub struct SiblingAnnexBuilder {
mmap: MmapMut,
}
impl SiblingAnnexBuilder {
/// Create a new annex of `n` slots at `path`, pre-initialised to the
/// "not yet computed" sentinel (all-zero).
pub fn new(n: usize, path: &Path) -> io::Result<Self> {
let file_size = HEADER_SIZE + n * 2;
let file = OpenOptions::new()
.read(true).write(true).create(true).truncate(true)
.open(path)?;
file.set_len(file_size as u64)?;
let mut mmap = unsafe { MmapMut::map_mut(&file)? };
mmap[0..4].copy_from_slice(&MAGIC);
mmap[4..8].copy_from_slice(&[0u8; 4]);
mmap[8..16].copy_from_slice(&(n as u64).to_le_bytes());
// Data region left at 0 by `set_len`/mmap — the sentinel value.
Ok(Self { mmap })
}
pub fn set(&mut self, slot: usize, mask: FamilyMask) {
// Redundant concurrent writes from independent recomputation paths
// converge to the same encoded value for a given slot, so a plain
// store here is safe even without external synchronisation, as long
// as the write itself doesn't tear — true for a 2-byte-aligned
// `u16` store on every platform this project targets.
let off = HEADER_SIZE + slot * 2;
self.mmap[off..off + 2].copy_from_slice(&mask.encode().to_le_bytes());
}
pub fn close(self) -> io::Result<()> { self.mmap.flush() }
}
#[cfg(test)]
mod tests {
use super::*;
use tempfile::tempdir;
#[test]
fn sentinel_is_zero_and_unset_slots_read_as_uncomputed() {
let dir = tempdir().unwrap();
let path = dir.path().join("test.psib");
let builder = SiblingAnnexBuilder::new(4, &path).unwrap();
builder.close().unwrap();
let annex = SiblingAnnex::open(&path).unwrap();
for slot in 0..4 {
assert_eq!(annex.get(slot), None);
}
}
#[test]
fn roundtrip_all_valid_masks() {
let dir = tempdir().unwrap();
let path = dir.path().join("test.psib");
let mut builder = SiblingAnnexBuilder::new(4, &path).unwrap();
let masks = [
FamilyMask::EMPTY.with(0), // just A: family size 1
FamilyMask::EMPTY.with(0).with(3), // A + T: size 2
FamilyMask::EMPTY.with(1).with(2).with(3), // C+G+T: size 3
FamilyMask::EMPTY.with(0).with(1).with(2).with(3), // all 4
];
for (slot, mask) in masks.iter().enumerate() {
builder.set(slot, *mask);
}
builder.close().unwrap();
let annex = SiblingAnnex::open(&path).unwrap();
for (slot, mask) in masks.iter().enumerate() {
assert_eq!(annex.get(slot), Some(*mask));
}
assert_eq!(annex.get(0).unwrap().siblings(), 0);
assert_eq!(annex.get(1).unwrap().siblings(), 1);
assert_eq!(annex.get(2).unwrap().siblings(), 2);
assert_eq!(annex.get(3).unwrap().siblings(), 3);
assert_eq!(annex.get(3).unwrap().family_size(), 4);
}
#[test]
fn has_reflects_individual_bits() {
let mask = FamilyMask::EMPTY.with(0).with(2);
assert!(mask.has(0));
assert!(!mask.has(1));
assert!(mask.has(2));
assert!(!mask.has(3));
}
#[test]
fn with_layer_roundtrips_and_leaves_absent_bases_none() {
let mask = FamilyMask::EMPTY.with_layer(0, 0).with_layer(2, 5);
assert_eq!(mask.layer_value(0), Some(0));
assert_eq!(mask.layer_value(2), Some(5));
assert_eq!(mask.layer_value(1), None);
assert_eq!(mask.layer_value(3), None);
assert!(mask.has(0) && mask.has(2));
assert!(!mask.has(1) && !mask.has(3));
assert_eq!(mask.family_size(), 2);
}
#[test]
fn with_marks_present_without_a_layer() {
// Compat path: `with` still means "present", now reading back as
// `layer_value == Some(0)` (raw field value 1, i.e. "no layer
// recorded") — distinguishable from a real layer 0 only by a
// caller that already knows whether layers were ever wired in.
let mask = FamilyMask::EMPTY.with(1);
assert!(mask.has(1));
assert_eq!(mask.layer_value(1), Some(0));
}
#[test]
fn family_members_reconstructs_present_canonical_forms() {
use obikseq::{Kmer, Sequence};
const K: usize = 11;
obikseq::params::set_k(K);
// Centre (index 5) is 'C' (base 1) in this k-mer's own orientation.
let owner = Kmer::from_ascii(b"AACCGCTTAAG").unwrap().canonical();
let mask = FamilyMask::EMPTY.with_layer(1, 0).with_layer(2, 3); // C (own) + G present
let members: Vec<_> = mask.family_members(owner, K).collect();
assert_eq!(members.len(), 2, "only the 2 present members should be yielded");
assert!(members.iter().any(|(k, layer)| *k == owner && *layer == Some(0)));
assert!(members.iter().any(|(k, layer)| *k != owner && *layer == Some(3)));
}
}
+1 -1
View File
@@ -2,7 +2,6 @@ use std::sync::Arc;
use rayon::prelude::*; use rayon::prelude::*;
use obicompactvec::SiblingAnnex;
use obikpartitionner::KmerPartition; use obikpartitionner::KmerPartition;
use obisys::progress_bar; use obisys::progress_bar;
@@ -10,6 +9,7 @@ use obikindex::{OKIError, OKIResult};
use obikindex::KmerIndex; use obikindex::KmerIndex;
use super::ANNEX_FILE_NAME; use super::ANNEX_FILE_NAME;
use super::SiblingAnnex;
use super::cache::PartitionCache; use super::cache::PartitionCache;
use super::family_scan::scan_layer_families; use super::family_scan::scan_layer_families;
+1 -2
View File
@@ -1,7 +1,6 @@
use std::io::Write; use std::io::Write;
use std::path::Path; use std::path::Path;
use obicompactvec::{FamilyMask, SiblingAnnex};
use obikseq::{CanonicalKmer, Kmer, Sequence}; use obikseq::{CanonicalKmer, Kmer, Sequence};
use obilayeredmap::MphfLayer; use obilayeredmap::MphfLayer;
use obilayeredmap::meta::PartitionMeta; use obilayeredmap::meta::PartitionMeta;
@@ -16,7 +15,7 @@ use super::cardinality::CardinalityExt;
use super::distance::DistanceExt; use super::distance::DistanceExt;
use super::helpers::is_minorant; use super::helpers::is_minorant;
use super::stats::SiblingStatsExt; use super::stats::SiblingStatsExt;
use super::{ANNEX_FILE_NAME, INDEX_SUBDIR}; use super::{FamilyMask, SiblingAnnex, ANNEX_FILE_NAME, INDEX_SUBDIR};
// k must be >= 11 (project constraint, "k ∈ [11,31]"); k=11, level_max=1, // k must be >= 11 (project constraint, "k ∈ [11,31]"); k=11, level_max=1,
// theta=0.0 mirror `obiskbuilder`'s own tests (smaller k/level_max // theta=0.0 mirror `obiskbuilder`'s own tests (smaller k/level_max