diff --git a/src/Cargo.lock b/src/Cargo.lock index 3efc8866..1cafb680 100644 --- a/src/Cargo.lock +++ b/src/Cargo.lock @@ -1846,6 +1846,7 @@ dependencies = [ name = "obikphylo" version = "0.1.0" dependencies = [ + "memmap2", "ndarray", "obicompactvec", "obikindex", diff --git a/src/obicompactvec/src/lib.rs b/src/obicompactvec/src/lib.rs index 1f85a49e..9041ab78 100644 --- a/src/obicompactvec/src/lib.rs +++ b/src/obicompactvec/src/lib.rs @@ -7,7 +7,6 @@ mod intmatrix; mod layer_meta; mod meta; mod reader; -mod siblingannex; mod tempbitvec; mod tempintvec; mod views; @@ -19,7 +18,6 @@ pub use builder::PersistentCompactIntVecBuilder; pub use colgroup::{ColGroup, FilterMask, MatrixGroupOps, eval_filter_mask}; pub use intmatrix::{PersistentCompactIntMatrix, PersistentCompactIntMatrixBuilder, pack_compact_int_matrix}; pub use layer_meta::LayerMeta; -pub use siblingannex::{FamilyMask, SiblingAnnex, SiblingAnnexBuilder}; pub use reader::{PersistentCompactIntVec, Iter as CompactIntVecIter}; pub use tempbitvec::{TempBitVec, TempBitVecBuilder}; pub use tempintvec::{TempCompactIntVec, TempCompactIntVecBuilder}; diff --git a/src/obicompactvec/src/siblingannex.rs b/src/obicompactvec/src/siblingannex.rs deleted file mode 100644 index 80e72d56..00000000 --- a/src/obicompactvec/src/siblingannex.rs +++ /dev/null @@ -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 { - 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 { - 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::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 { - 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::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 { - 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)); - } -} diff --git a/src/obikphylo/Cargo.toml b/src/obikphylo/Cargo.toml index be66a13c..5dc257ba 100644 --- a/src/obikphylo/Cargo.toml +++ b/src/obikphylo/Cargo.toml @@ -13,6 +13,7 @@ obicompactvec = { path = "../obicompactvec" } obilayeredmap = { path = "../obilayeredmap" } obiskbuilder = { path = "../obiskbuilder" } obipipeline = { path = "../obipipeline" } +memmap2 = "0.9" ndarray = "0.16" rayon = "1" tracing = "0.1.44" diff --git a/src/obikphylo/src/siblings/build.rs b/src/obikphylo/src/siblings/build.rs index b81d160b..5a541e5f 100644 --- a/src/obikphylo/src/siblings/build.rs +++ b/src/obikphylo/src/siblings/build.rs @@ -4,7 +4,6 @@ use std::sync::Arc; use rayon::prelude::*; -use obicompactvec::{FamilyMask, SiblingAnnexBuilder}; use obikpartitionner::KmerPartition; use obipipeline::ThrottleGuard; use obikseq::CanonicalKmer; @@ -17,7 +16,7 @@ use obikindex::KmerIndex; use super::cache::PartitionCache; 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 ───────────────────────────────────────────────── diff --git a/src/obikphylo/src/siblings/cache.rs b/src/obikphylo/src/siblings/cache.rs index 8d706d27..e3047821 100644 --- a/src/obikphylo/src/siblings/cache.rs +++ b/src/obikphylo/src/siblings/cache.rs @@ -10,7 +10,7 @@ use obisys::progress_bar; use obikindex::OKIResult; 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 /// `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. pub(super) fn iter_minorants_batch( &self, - annex: std::sync::Arc, + annex: std::sync::Arc, batch_size: usize, ) -> super::iter::MinorantBatchIter { match self { diff --git a/src/obikphylo/src/siblings/family_scan.rs b/src/obikphylo/src/siblings/family_scan.rs index 56bc5d4c..fdde2249 100644 --- a/src/obikphylo/src/siblings/family_scan.rs +++ b/src/obikphylo/src/siblings/family_scan.rs @@ -53,7 +53,7 @@ use std::sync::atomic::{AtomicU8, Ordering}; use rayon::prelude::*; -use obicompactvec::{FamilyMask, PersistentBitMatrix, PersistentCompactIntMatrix, SiblingAnnex}; +use obicompactvec::{PersistentBitMatrix, PersistentCompactIntMatrix}; use obikseq::CanonicalKmer; use obilayeredmap::Layer; use obilayeredmap::meta::PartitionMeta; @@ -65,7 +65,7 @@ use obikindex::KmerIndex; use super::cache::{Mat, PartitionCache}; use super::helpers::central_base; 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- /// density trade-off this picks a point on. At ~90 genomes and a few diff --git a/src/obikphylo/src/siblings/helpers.rs b/src/obikphylo/src/siblings/helpers.rs index 580162b4..b11c7c52 100644 --- a/src/obikphylo/src/siblings/helpers.rs +++ b/src/obikphylo/src/siblings/helpers.rs @@ -1,6 +1,6 @@ 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 /// encoding — the mask's bit index. `k` must be odd (project invariant). diff --git a/src/obikphylo/src/siblings/iter.rs b/src/obikphylo/src/siblings/iter.rs index 20bea078..932f0eed 100644 --- a/src/obikphylo/src/siblings/iter.rs +++ b/src/obikphylo/src/siblings/iter.rs @@ -28,10 +28,11 @@ use std::sync::Arc; -use obicompactvec::{FamilyMask, SiblingAnnex}; use obikseq::CanonicalKmer; use obilayeredmap::{KmerIter, Layer, LayerData}; +use super::{FamilyMask, SiblingAnnex}; + /// 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 /// itself, and its family mask. diff --git a/src/obikphylo/src/siblings/mod.rs b/src/obikphylo/src/siblings/mod.rs index cd3f5da8..9fa04f5e 100644 --- a/src/obikphylo/src/siblings/mod.rs +++ b/src/obikphylo/src/siblings/mod.rs @@ -50,6 +50,7 @@ mod distance; mod family_scan; mod helpers; mod iter; +mod siblingannex; mod stats; #[cfg(test)] @@ -60,6 +61,7 @@ pub use build::SiblingAnnexBuildExt; pub use cardinality::{CardinalityExt, CardinalityTally}; pub use distance::{BasePairTally, DistanceExt, RawSnpDistanceOutput}; pub use iter::{MinorantBatchIter, MinorantIter, SiblingBatchIter, SiblingEntry, SiblingIter, SiblingLayerExt}; +pub(crate) use siblingannex::{FamilyMask, SiblingAnnex, SiblingAnnexBuilder}; pub use stats::{SiblingAnnexStats, SiblingStatsExt}; use obilayeredmap::OLMError; diff --git a/src/obikphylo/src/siblings/siblingannex.rs b/src/obikphylo/src/siblings/siblingannex.rs new file mode 100644 index 00000000..38f62aee --- /dev/null +++ b/src/obikphylo/src/siblings/siblingannex.rs @@ -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 { + 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)> { + 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 { + 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 { + 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 { + 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 { + 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))); + } +} diff --git a/src/obikphylo/src/siblings/stats.rs b/src/obikphylo/src/siblings/stats.rs index a3edd1f2..e612257f 100644 --- a/src/obikphylo/src/siblings/stats.rs +++ b/src/obikphylo/src/siblings/stats.rs @@ -2,7 +2,6 @@ use std::sync::Arc; use rayon::prelude::*; -use obicompactvec::SiblingAnnex; use obikpartitionner::KmerPartition; use obisys::progress_bar; @@ -10,6 +9,7 @@ use obikindex::{OKIError, OKIResult}; use obikindex::KmerIndex; use super::ANNEX_FILE_NAME; +use super::SiblingAnnex; use super::cache::PartitionCache; use super::family_scan::scan_layer_families; diff --git a/src/obikphylo/src/siblings/tests.rs b/src/obikphylo/src/siblings/tests.rs index d854bcab..702d8365 100644 --- a/src/obikphylo/src/siblings/tests.rs +++ b/src/obikphylo/src/siblings/tests.rs @@ -1,7 +1,6 @@ use std::io::Write; use std::path::Path; -use obicompactvec::{FamilyMask, SiblingAnnex}; use obikseq::{CanonicalKmer, Kmer, Sequence}; use obilayeredmap::MphfLayer; use obilayeredmap::meta::PartitionMeta; @@ -16,7 +15,7 @@ use super::cardinality::CardinalityExt; use super::distance::DistanceExt; use super::helpers::is_minorant; 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, // theta=0.0 mirror `obiskbuilder`'s own tests (smaller k/level_max