diff --git a/src/obikmer2/src/cmd/phylo/args.rs b/src/obikmer2/src/cmd/phylo/args.rs index 5ae930d5..d478463b 100644 --- a/src/obikmer2/src/cmd/phylo/args.rs +++ b/src/obikmer2/src/cmd/phylo/args.rs @@ -82,10 +82,10 @@ impl DistanceArg { /// `--free-loss`, `--no-ambiguity`, `--entropy`/`--entropy-sd`), Sankoff /// cost-matrix calibration (`--sankoff`, `--sankoff-ratio-ceiling`) and its /// TNT/PhyG/IQ-TREE exports (`--tnt`, `--phyg`, `--iqtree`/ -/// `--iqtree-min-freq`, `--sankoff-cost-scale`) — everything else -/// sibling-annex-based (family overlap, ...) stays in `obikmer` until the -/// rest of `obikphylo::siblings` is reconnected (see the project memory on -/// this). +/// `--iqtree-min-freq`, `--sankoff-cost-scale`), and Family Overlap +/// (`--family-overlap`, `--min-shared-family`) — everything else +/// sibling-annex-based stays in `obikmer` until the rest of +/// `obikphylo::siblings` is reconnected (see the project memory on this). #[derive(Args)] pub struct PhyloArgs { /// Index directory @@ -101,6 +101,23 @@ pub struct PhyloArgs { #[arg(long = "exclude-genome", value_name = "LABEL")] pub exclude_genome: Vec, + /// Auto-exclude any genome whose mean shared-variable-family count + /// against every other genome (`FamilyOverlap::mean_row` — the same + /// per-row statistic `--family-overlap`'s own matrix shows) falls below + /// this threshold — same exclusion machinery as `--exclude-genome`, + /// applied on top of it rather than instead of it. Applies to the + /// `snp-*` `--distance`/`--pseudo-alignment`/`--sankoff` computations + /// below (all sibling-annex-based); does *not* affect the whole-index + /// `--distance` metrics (jaccard, hamming, bray-curtis, ...) or their + /// matrix/NJ/UPGMA output — a genome with too little SNP-family + /// coverage to trust is a different concern from one whose plain k-mer + /// profile is simply divergent. Requires the Family Overlap annex + /// (built on demand if missing, same as every other annex here — see + /// `obikphylo::siblings::extensions::SiblingExt::family_overlap`'s own + /// docs). + #[arg(long, value_name = "N")] + pub min_shared_family: Option, + /// Build (or rebuild) the sibling-count/minorant annex — independent of /// the distance metric below, meant to be run routinely, ahead of any /// SNP-family distance computation that will later consume it. @@ -122,6 +139,16 @@ pub struct PhyloArgs { #[arg(long)] pub sibling_hist: bool, + /// Write the Family Overlap matrix (CSV) — number of shared *variable* + /// families per genome pair, index-wide (`obikphylo::siblings::FamilyOverlap`). + /// Built on demand if missing (same as every other annex here); no + /// `--sibling-annex` prerequisite beyond that. Every genome is written, + /// unfiltered by `--exclude-genome`/`--min-shared-family` — a raw + /// coverage diagnostic, not a computation those exclusions are meant to + /// protect. + #[arg(long)] + pub family_overlap: bool, + /// Write a per-family Shannon entropy report (CSV) — requires an /// already-built sibling annex (`--sibling-annex` first, in this /// invocation or an earlier one). Always a full, unsampled scan of diff --git a/src/obikmer2/src/cmd/phylo/mod.rs b/src/obikmer2/src/cmd/phylo/mod.rs index fc2d62f8..e61f91f0 100644 --- a/src/obikmer2/src/cmd/phylo/mod.rs +++ b/src/obikmer2/src/cmd/phylo/mod.rs @@ -151,6 +151,33 @@ pub fn run(args: PhyloArgs) { ); } + // ── Family Overlap matrix (`--family-overlap`) ────────────────────────────── + if args.family_overlap { + let t = Stage::start("family_overlap"); + let overlap = cache.family_overlap().unwrap_or_else(|e| { + eprintln!("error computing family overlap: {e}"); + std::process::exit(1); + }); + rep.push(t.stop()); + + let path = args.output.as_ref() + .map(|p| format!("{}_family_overlap.csv", p.display())) + .unwrap_or_else(|| "family_overlap.csv".into()); + let mut f = BufWriter::new(std::fs::File::create(&path).unwrap_or_else(|e| { + eprintln!("error creating {path}: {e}"); + std::process::exit(1); + })); + write!(f, "genome").unwrap(); + for label in &labels { write!(f, ",{label}").unwrap(); } + writeln!(f).unwrap(); + for (i, label) in labels.iter().enumerate() { + write!(f, "{label}").unwrap(); + for j in 0..n { write!(f, ",{}", overlap.get(i, j)).unwrap(); } + writeln!(f).unwrap(); + } + info!("family-overlap matrix → {path}"); + } + // ── Shannon entropy report (`--shannon`) ──────────────────────────────────── if args.shannon { let path = args.output.as_ref() @@ -166,6 +193,37 @@ pub fn run(args: PhyloArgs) { info!("entropy report → {path}"); } + // ── `--min-shared-family` auto-exclusion ──────────────────────────────────── + // Layered on top of `--exclude-genome`, not instead of it — a separate + // mask (not folded into `exclude_mask` itself) since it must reach + // `--pseudo-alignment`/`--sankoff`/`snp-*` `--distance` only, never the + // whole-index metrics' `kept`/`--shared-kmers` output (see + // `args::PhyloArgs::min_shared_family`'s own docs on why). + let snp_exclude_mask: Vec = match args.min_shared_family { + Some(threshold) => { + let overlap = cache.family_overlap().unwrap_or_else(|e| { + eprintln!("error computing family overlap: {e}"); + std::process::exit(1); + }); + let mut mask = exclude_mask.clone(); + for (g, excluded) in mask.iter_mut().enumerate() { + if *excluded { + continue; + } + let mean = overlap.mean_row(g); + if mean < threshold { + info!( + "auto-excluding {} (--min-shared-family: mean shared-family count {mean:.1} < {threshold})", + labels[g] + ); + *excluded = true; + } + } + mask + } + None => exclude_mask.clone(), + }; + // Shared by `--pseudo-alignment` and `--sankoff` — same activation rule: // either flag given activates entropy-biased sampling, the other // defaults to 1.0/0.5. @@ -188,7 +246,7 @@ pub fn run(args: PhyloArgs) { info!("sampling SNP pseudo-alignment (target {subsample_n} site(s))"); let t = Stage::start("pseudo_alignment"); let alignment = cache - .snp_pseudo_alignment(subsample_n, args.free_loss, args.no_ambiguity, &exclude_mask, entropy_bias) + .snp_pseudo_alignment(subsample_n, args.free_loss, args.no_ambiguity, &snp_exclude_mask, entropy_bias) .unwrap_or_else(|e| { eprintln!("error building pseudo-alignment: {e}"); std::process::exit(1); @@ -223,7 +281,7 @@ pub fn run(args: PhyloArgs) { subsample_n, args.free_loss, args.no_ambiguity, - &exclude_mask, + &snp_exclude_mask, entropy_bias, args.sankoff_ratio_ceiling, ) @@ -310,7 +368,7 @@ pub fn run(args: PhyloArgs) { args.subsample, args.free_loss, args.no_ambiguity, - &exclude_mask, + &snp_exclude_mask, entropy_bias, args.gamma_shape, ) diff --git a/src/obikphylo/src/siblings/algorithms/entropy.rs b/src/obikphylo/src/siblings/algorithms/entropy.rs index c1baad05..175d46ac 100644 --- a/src/obikphylo/src/siblings/algorithms/entropy.rs +++ b/src/obikphylo/src/siblings/algorithms/entropy.rs @@ -16,7 +16,10 @@ use obikindex::{OKIError, OKIResult}; use super::family_scan::{Selection, scan_layer_families}; use super::minorant_selection::{is_non_monomorphic_minorant, non_monomorphic_reindex_layer}; -use crate::siblings::{ANNEX_FILE_NAME, ENTROPY_ANNEX_FILE_NAME, EntropyAnnex, EntropyAnnexBuilder, SiblingAnnex}; +use crate::siblings::{ + ANNEX_FILE_NAME, ENTROPY_ANNEX_FILE_NAME, EntropyAnnex, EntropyAnnexBuilder, + FamilyOverlapAccumulator, SiblingAnnex, +}; /// Shannon entropy (bits) of one family's per-genome states, over the 15 /// non-empty subsets of `{A,C,G,T}` actually observed among the genomes @@ -83,31 +86,45 @@ pub(crate) fn family_entropy_4(genome_mask: &[u8]) -> Option<(f64, usize)> { Some((h, total as usize)) } -/// Builds one layer's entropy annex if it doesn't already exist — a no-op, -/// cheap existence check, once it does. Bounded to non-monomorphic -/// minorants only ([`non_monomorphic_reindex_layer`], a cheap annex-only -/// pass — no cross-partition resolution): monomorphism is already known -/// from the annex bits alone, so there is no reason to pay +/// Builds one layer's entropy annex if it doesn't already exist, and/or +/// feeds `overlap_acc` (`Some` iff the index-wide family-overlap annex is +/// itself missing — see +/// [`crate::siblings::FamilyOverlap`]'s own module docs) — a single shared +/// `scan_layer_families` pass covers both, since both need the exact same +/// per-family `genome_mask` resolution and there is no reason to pay for it +/// twice just because only one of the two happens to be missing. A no-op +/// (no scan at all) when entropy already exists for this layer *and* +/// `overlap_acc` is `None`. Bounded to non-monomorphic minorants only +/// ([`non_monomorphic_reindex_layer`], a cheap annex-only pass — no +/// cross-partition resolution) whenever a scan does happen: monomorphism is +/// already known from the annex bits alone, so there is no reason to pay /// `scan_layer_families`'s expensive per-genome resolution for the ~98% of -/// minorants that are monomorphic only to discard the result. +/// minorants that are monomorphic only to discard the result — true for +/// family-overlap's own "variable families only" scope just as much as for +/// entropy's. pub(crate) fn ensure_layer_entropy_annex( cache: &IndexCache, partition: usize, layer_idx: usize, n_genomes: usize, fast_mode: bool, + mut overlap_acc: Option<&mut FamilyOverlapAccumulator>, ) -> OKIResult<()> { let layer = cache .get_layer(partition, layer_idx) .expect("cache-consistent: caller already bounds layer_idx by n_layer(partition)"); let path = layer.dir().join(ENTROPY_ANNEX_FILE_NAME); - if path.exists() { + let need_entropy = !path.exists(); + if !need_entropy && overlap_acc.is_none() { return Ok(()); } let reindex = non_monomorphic_reindex_layer(layer.dir())?; - let mut builder = EntropyAnnexBuilder::new(reindex.len(), &path).map_err(OKIError::Io)?; let eligible: HashSet = reindex.keys().copied().collect(); + let mut builder = need_entropy + .then(|| EntropyAnnexBuilder::new(reindex.len(), &path)) + .transpose() + .map_err(OKIError::Io)?; scan_layer_families( cache, @@ -121,13 +138,19 @@ pub(crate) fn ensure_layer_entropy_annex( mask.family_size() >= 2, "Selection::Some(eligible) must only visit non-monomorphic minorants" ); - if let Some((h15, _)) = family_entropy(genome_mask) { - let compact = reindex[&family_idx]; - builder.set(compact, h15 as f32); + if let Some(builder) = builder.as_mut() { + if let Some((h15, _)) = family_entropy(genome_mask) { + builder.set(reindex[&family_idx], h15 as f32); + } + } + if let Some(acc) = overlap_acc.as_deref_mut() { + acc.add_family(genome_mask); } }, )?; - builder.close().map_err(OKIError::Io)?; + if let Some(builder) = builder { + builder.close().map_err(OKIError::Io)?; + } Ok(()) } diff --git a/src/obikphylo/src/siblings/algorithms/pairwise.rs b/src/obikphylo/src/siblings/algorithms/pairwise.rs index 07e5cff1..b57f49c1 100644 --- a/src/obikphylo/src/siblings/algorithms/pairwise.rs +++ b/src/obikphylo/src/siblings/algorithms/pairwise.rs @@ -52,15 +52,9 @@ impl PairwiseTally { Self { n_genomes, pairs: vec![PairStats::default(); n_pairs] } } - /// Flat index of unordered pair `(i, j)`, `i != j`, into the upper - /// triangle — standard "row-major over `i < j`" packing: row `i` - /// contributes `n_genomes - 1 - i` entries (columns `i+1..n_genomes`), - /// so pair `(i, j)` sits at `sum_{r usize { - let (i, j) = if i < j { (i, j) } else { (j, i) }; - debug_assert!(j < self.n_genomes && i != j); - let n = self.n_genomes; - (i * (2 * n - i - 1)) / 2 + (j - i - 1) + crate::siblings::helpers::triangle_index(self.n_genomes, i, j) } #[inline] diff --git a/src/obikphylo/src/siblings/algorithms/subsample.rs b/src/obikphylo/src/siblings/algorithms/subsample.rs index 840d08ca..dd52afc7 100644 --- a/src/obikphylo/src/siblings/algorithms/subsample.rs +++ b/src/obikphylo/src/siblings/algorithms/subsample.rs @@ -142,7 +142,7 @@ pub(crate) fn sample_index( // annex (`EntropyWeights::open`, per layer, below) — build any missing // one now, once, rather than mid-walk. if entropy_bias.is_some() { - cache.ensure_entropy_annexes()?; + cache.ensure_entropy_and_overlap_annexes()?; } // Every cached layer's eligibility bitset, built up front — annex-only, diff --git a/src/obikphylo/src/siblings/extensions/sibling_builder.rs b/src/obikphylo/src/siblings/extensions/sibling_builder.rs index b8915941..1464917e 100644 --- a/src/obikphylo/src/siblings/extensions/sibling_builder.rs +++ b/src/obikphylo/src/siblings/extensions/sibling_builder.rs @@ -1,40 +1,53 @@ //! [`SiblingBuilder`] — construction triggered *internally*, as a side //! effect of a [`crate::siblings::extensions::SiblingExt`] method (e.g. -//! `shannon_entropy_csv` ensuring its layers' entropy annexes exist before -//! reading them), never called directly from outside this crate. Unlike -//! `SiblingExt` (explicit, CLI-facing actions), nothing here is meant to be -//! its own CLI flag — same `pub`/`pub(crate)` split `obikindexer` uses for -//! its own `IndexBuilder`/`PrivateBuilder`. +//! `shannon_entropy_csv`/`family_overlap` ensuring their annexes exist +//! before reading them), never called directly from outside this crate. +//! Unlike `SiblingExt` (explicit, CLI-facing actions), nothing here is +//! meant to be its own CLI flag — same `pub`/`pub(crate)` split +//! `obikindexer` uses for its own `IndexBuilder`/`PrivateBuilder`. use obikidxcache::index_cache::IndexCache; -use obikindex::OKIResult; +use obikindex::{OKIError, OKIResult}; use obisys::progress_bar; use crate::siblings::algorithms::{ensure_layer_entropy_annex, is_fast_mode}; +use crate::siblings::{FAMILY_OVERLAP_FILE_NAME, FamilyOverlapAccumulator}; pub(crate) trait SiblingBuilder { - /// Build every cached layer's entropy annex that doesn't already exist - /// — a no-op, cheap existence check per layer, once every layer's - /// annex file exists. See `ensure_layer_entropy_annex`'s own docs for - /// why this is bounded to non-monomorphic minorants only. - fn ensure_entropy_annexes(&self) -> OKIResult<()>; + /// Build every cached layer's entropy annex, and the whole index's + /// family-overlap annex, that don't already exist — one call covers + /// both, whichever of the two a caller actually needs + /// (`shannon_entropy_csv`/entropy-biased sampling on one side, + /// `family_overlap`/`--min-shared-family` on the other): a missing + /// annex always triggers a scan that feeds both, per + /// `ensure_layer_entropy_annex`'s own docs. A no-op, cheap + /// existence-check pass, once both already exist everywhere. + fn ensure_entropy_and_overlap_annexes(&self) -> OKIResult<()>; } impl SiblingBuilder for IndexCache { - fn ensure_entropy_annexes(&self) -> OKIResult<()> { + fn ensure_entropy_and_overlap_annexes(&self) -> OKIResult<()> { let n_genomes = self.meta().genomes().len(); let fast_mode = is_fast_mode(self); - let pb = progress_bar("entropy_annex", self.n_partition() as u64, "partitions"); + let overlap_path = self.index().dir().join(FAMILY_OVERLAP_FILE_NAME); + let mut overlap_acc = + (!overlap_path.exists()).then(|| FamilyOverlapAccumulator::new(n_genomes)); + + let pb = progress_bar("entropy_overlap_annex", self.n_partition() as u64, "partitions"); for part in self.partitions() { let n_layer = self.n_layer(part).unwrap_or(0); for l in 0..n_layer { - ensure_layer_entropy_annex(self, part, l, n_genomes, fast_mode)?; + ensure_layer_entropy_annex(self, part, l, n_genomes, fast_mode, overlap_acc.as_mut())?; } pb.inc(1); } pb.finish_and_clear(); + if let Some(acc) = overlap_acc { + acc.write(&overlap_path).map_err(OKIError::Io)?; + } + Ok(()) } } diff --git a/src/obikphylo/src/siblings/extensions/sibling_ext.rs b/src/obikphylo/src/siblings/extensions/sibling_ext.rs index cc9e1177..69d1c1f7 100644 --- a/src/obikphylo/src/siblings/extensions/sibling_ext.rs +++ b/src/obikphylo/src/siblings/extensions/sibling_ext.rs @@ -21,7 +21,7 @@ use crate::siblings::algorithms::{ snp_distance, snp_pseudo_alignment, }; use crate::siblings::extensions::SiblingBuilder; -use crate::siblings::ENTROPY_ANNEX_FILE_NAME; +use crate::siblings::{ENTROPY_ANNEX_FILE_NAME, FAMILY_OVERLAP_FILE_NAME, FamilyOverlap}; pub trait SiblingExt { /// Build the sibling-count/minorant annex for every layer of every @@ -154,6 +154,18 @@ pub trait SiblingExt { entropy_bias: Option, gamma_shape: Option, ) -> OKIResult>; + + /// The Family Overlap annex — number of shared *variable* families per + /// genome pair, index-wide (see [`FamilyOverlap`]'s own module docs). + /// Built (alongside the entropy annex) if missing, same as every other + /// annex-backed method here. Feeds `--family-overlap`'s own matrix + /// output directly, and `--min-shared-family`'s auto-exclusion + /// (`FamilyOverlap::mean_row(g) < threshold` — same statistic + /// `--family-overlap`'s matrix shows, just averaged per row — + /// triggers adding `g` to the excluded-genome set consumed by + /// `snp_pseudo_alignment`/`sankoff_bundle`/`snp_distance`; this needs + /// no computation of its own beyond that check). + fn family_overlap(&self) -> OKIResult; } impl SiblingExt for IndexCache { @@ -181,10 +193,10 @@ impl SiblingExt for IndexCache { // Invalidate this layer's entropy annex, if any — it was // built against the sibling annex's *previous* // `family_idx` numbering (`non_monomorphic_reindex_layer` - // reads the sibling annex), now stale. `ensure_entropy_annexes` - // below only builds what's *missing*, so the old file must - // be removed here for it to be rebuilt rather than silently - // left stale. + // reads the sibling annex), now stale. + // `ensure_entropy_and_overlap_annexes` below only builds + // what's *missing*, so the old file must be removed here + // for it to be rebuilt rather than silently left stale. let entropy_path = layer.dir().join(ENTROPY_ANNEX_FILE_NAME); if entropy_path.exists() { std::fs::remove_file(&entropy_path).map_err(OKIError::Io)?; @@ -199,24 +211,36 @@ impl SiblingExt for IndexCache { pb.finish_and_clear(); tracing::info!("sibling annex built — {total_slots} kmers across {n_parts} partitions"); - // Precompute the entropy annex too, now that every layer's fresh - // sibling annex is on disk (`ensure_entropy_annexes` reads it via - // `non_monomorphic_reindex_layer`) — same routine-rebuild footing - // as the sibling annex itself, rather than leaving it to be built - // lazily on first `--shannon`/`--entropy-bias` use. - tracing::info!("building entropy annex"); - self.ensure_entropy_annexes()?; + // The family-overlap annex is index-level, not per-layer (see its + // own module docs) — any layer's sibling annex changing makes the + // whole thing stale, so it's invalidated once here rather than + // per-layer above. + let overlap_path = self.index().dir().join(FAMILY_OVERLAP_FILE_NAME); + if overlap_path.exists() { + std::fs::remove_file(&overlap_path).map_err(OKIError::Io)?; + } + + // Precompute the entropy and family-overlap annexes too, now that + // every layer's fresh sibling annex is on disk + // (`non_monomorphic_reindex_layer`, which both depend on, reads + // it) — same routine-rebuild footing as the sibling annex itself, + // rather than leaving them to be built lazily on first + // `--shannon`/`--entropy-bias`/`--min-shared-family`/ + // `--family-overlap` use. + tracing::info!("building entropy and family-overlap annexes"); + self.ensure_entropy_and_overlap_annexes()?; Ok(()) } fn shannon_entropy_csv(&self, path: &Path) -> OKIResult<()> { - // Entropy annexes aren't this method's own output, but building - // them is exactly the same full scan this method is about to run - // anyway — see `SiblingBuilder::ensure_entropy_annexes`'s own docs + // Entropy/family-overlap annexes aren't this method's own output, + // but building them is exactly the same full scan this method is + // about to run anyway — see + // `SiblingBuilder::ensure_entropy_and_overlap_annexes`'s own docs // for why that's a side effect worth taking here rather than a // separate, explicit step. - self.ensure_entropy_annexes()?; + self.ensure_entropy_and_overlap_annexes()?; let n_genomes = self.meta().genomes().len(); let fast_mode = is_fast_mode(self); @@ -311,4 +335,10 @@ impl SiblingExt for IndexCache { ) -> OKIResult> { snp_distance(self, kind, n, free_loss, no_ambiguity, excluded, entropy_bias, gamma_shape) } + + fn family_overlap(&self) -> OKIResult { + self.ensure_entropy_and_overlap_annexes()?; + let path = self.index().dir().join(FAMILY_OVERLAP_FILE_NAME); + FamilyOverlap::open(&path).map_err(OKIError::Io) + } } diff --git a/src/obikphylo/src/siblings/family_overlap.rs b/src/obikphylo/src/siblings/family_overlap.rs new file mode 100644 index 00000000..a1e1ed2e --- /dev/null +++ b/src/obikphylo/src/siblings/family_overlap.rs @@ -0,0 +1,167 @@ +//! Family Overlap annex — the number of shared *variable* families +//! (`FamilyMask::family_size() >= 2`) per genome pair, summed across the +//! **whole index**. See `DevDocMD/theory/evolutionary_distances.md`, +//! "Family Overlap as an index-level annex", for the design discussion: +//! unlike [`crate::siblings::SiblingAnnex`]/[`crate::siblings::EntropyAnnex`] +//! (sharded per layer because they hold one entry per k-mer — potentially +//! billions, impossible to keep in memory), this structure's size depends +//! only on `n_genomes²`, never on index size, so there is no memory-pressure +//! reason to shard it. It lives once, at the index root +//! (`IndexCache::index().dir()`), not inside any layer directory. +//! +//! Built alongside the entropy annex, in the very same +//! `scan_layer_families` sweep — see +//! `algorithms::entropy::ensure_layer_entropy_annex`'s own docs: both need +//! the identical per-family `genome_mask` resolution, so paying for it once +//! and feeding both is (almost) free, whichever of the two triggered the +//! sweep. The only case that still costs a full rescan for just one of them +//! is a pre-existing index where the other annex already happens to be +//! present — accepted as unavoidable, not worth special-casing. + +use std::fs::File; +use std::io; +use std::path::Path; + +use memmap2::Mmap; + +use crate::siblings::helpers::triangle_index; + +const MAGIC: [u8; 4] = *b"POVL"; + +// Header: magic(4) + _pad(4) + n_genomes(8) = 16 bytes. Data (8 bytes/pair, +// upper triangle) follows. +const HEADER_SIZE: usize = 16; + +pub(crate) const FAMILY_OVERLAP_FILE_NAME: &str = "family_overlap.povl"; + +/// Read-only, mmap-backed view of an already-built Family Overlap annex. +pub struct FamilyOverlap { + mmap: Mmap, + n_genomes: usize, +} + +impl FamilyOverlap { + 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, "POVL file too short")); + } + if mmap[0..4] != MAGIC { + return Err(io::Error::new(io::ErrorKind::InvalidData, "bad POVL magic")); + } + let n_genomes = u64::from_le_bytes(mmap[8..16].try_into().unwrap()) as usize; + let n_pairs = n_genomes.saturating_sub(1) * n_genomes / 2; + if mmap.len() < HEADER_SIZE + n_pairs * 8 { + return Err(io::Error::new(io::ErrorKind::InvalidData, "POVL file truncated")); + } + Ok(Self { mmap, n_genomes }) + } + + pub fn n_genomes(&self) -> usize { + self.n_genomes + } + + /// Number of variable families shared by genomes `i` and `j` — `0` for + /// `i == j` (not stored; a genome trivially "shares" every one of its + /// own families with itself, not a meaningful statistic here). + pub fn get(&self, i: usize, j: usize) -> u64 { + if i == j { + return 0; + } + let off = HEADER_SIZE + triangle_index(self.n_genomes, i, j) * 8; + u64::from_le_bytes(self.mmap[off..off + 8].try_into().unwrap()) + } + + /// Mean shared-family count of genome `i` against every other genome — + /// `--min-shared-family`'s own statistic (see + /// `crate::siblings::extensions::SiblingExt::family_overlap`'s docs). + /// `0.0` for a single-genome index (no other genome to average + /// against). + pub fn mean_row(&self, i: usize) -> f64 { + if self.n_genomes <= 1 { + return 0.0; + } + let sum: u64 = (0..self.n_genomes).filter(|&j| j != i).map(|j| self.get(i, j)).sum(); + sum as f64 / (self.n_genomes - 1) as f64 + } +} + +/// In-memory accumulator, filled while sweeping every layer of the index +/// (one `add_family` call per variable family visited, from whichever +/// `scan_layer_families` pass is already running for the entropy annex — +/// see the module docs), written once as a whole file when the sweep +/// completes. Small enough (`O(n_genomes²)`) that there is no reason to +/// build it incrementally on disk the way `SiblingAnnexBuilder`/ +/// `EntropyAnnexBuilder` do for their own, much larger, per-k-mer data. +pub(crate) struct FamilyOverlapAccumulator { + n_genomes: usize, + counts: Vec, +} + +impl FamilyOverlapAccumulator { + pub(crate) fn new(n_genomes: usize) -> Self { + let n_pairs = n_genomes.saturating_sub(1) * n_genomes / 2; + Self { n_genomes, counts: vec![0u64; n_pairs] } + } + + /// One variable family's contribution: every pair of genomes both + /// carrying it (`genome_mask[g] != 0`) gets its count incremented once. + pub(crate) fn add_family(&mut self, genome_mask: &[u8]) { + for i in 0..self.n_genomes { + if genome_mask[i] == 0 { + continue; + } + for j in (i + 1)..self.n_genomes { + if genome_mask[j] != 0 { + self.counts[triangle_index(self.n_genomes, i, j)] += 1; + } + } + } + } + + pub(crate) fn write(&self, path: &Path) -> io::Result<()> { + let mut buf = Vec::with_capacity(HEADER_SIZE + self.counts.len() * 8); + buf.extend_from_slice(&MAGIC); + buf.extend_from_slice(&[0u8; 4]); + buf.extend_from_slice(&(self.n_genomes as u64).to_le_bytes()); + for &c in &self.counts { + buf.extend_from_slice(&c.to_le_bytes()); + } + std::fs::write(path, buf) + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn roundtrip_and_mean_row() { + let dir = tempfile::tempdir().unwrap(); + let path = dir.path().join("test.povl"); + + let mut acc = FamilyOverlapAccumulator::new(3); + acc.add_family(&[1, 1, 0]); // genomes 0,1 share this family + acc.add_family(&[1, 1, 1]); // all three share this one + acc.add_family(&[0, 1, 1]); // genomes 1,2 share this family + acc.write(&path).unwrap(); + + let overlap = FamilyOverlap::open(&path).unwrap(); + assert_eq!(overlap.n_genomes(), 3); + assert_eq!(overlap.get(0, 1), 2); + assert_eq!(overlap.get(0, 2), 1); + assert_eq!(overlap.get(1, 2), 2); + assert_eq!(overlap.get(0, 0), 0); + assert_eq!(overlap.get(1, 0), overlap.get(0, 1)); // symmetric + assert!((overlap.mean_row(0) - 1.5).abs() < 1e-12); // (2+1)/2 + } + + #[test] + fn single_genome_mean_row_is_zero() { + let dir = tempfile::tempdir().unwrap(); + let path = dir.path().join("test.povl"); + FamilyOverlapAccumulator::new(1).write(&path).unwrap(); + let overlap = FamilyOverlap::open(&path).unwrap(); + assert_eq!(overlap.mean_row(0), 0.0); + } +} diff --git a/src/obikphylo/src/siblings/helpers.rs b/src/obikphylo/src/siblings/helpers.rs index d6e9ab8f..b9d92499 100644 --- a/src/obikphylo/src/siblings/helpers.rs +++ b/src/obikphylo/src/siblings/helpers.rs @@ -20,3 +20,17 @@ pub(super) fn is_minorant(kmer: CanonicalKmer, mask: FamilyMask, k: usize) -> bo other == kmer || !mask.has(central_base(other, k)) || kmer.raw() <= other.raw() }) } + +/// Flat index of unordered genome pair `(i, j)`, `i != j`, into an +/// upper-triangle-packed `n`-genome array — standard "row-major over `i < +/// j`" packing: row `i` contributes `n - 1 - i` entries (columns +/// `i+1..n`), so pair `(i, j)` sits at `sum_{r usize { + let (i, j) = if i < j { (i, j) } else { (j, i) }; + debug_assert!(j < n && i != j); + (i * (2 * n - i - 1)) / 2 + (j - i - 1) +} diff --git a/src/obikphylo/src/siblings/mod.rs b/src/obikphylo/src/siblings/mod.rs index 42c1d809..d5067f96 100644 --- a/src/obikphylo/src/siblings/mod.rs +++ b/src/obikphylo/src/siblings/mod.rs @@ -22,12 +22,15 @@ pub mod algorithms; mod entropy_annex; pub mod extensions; +mod family_overlap; mod helpers; pub mod iter; mod siblingannex; pub(crate) use entropy_annex::ENTROPY_ANNEX_FILE_NAME; pub use entropy_annex::{EntropyAnnex, EntropyAnnexBuilder}; +pub(crate) use family_overlap::{FAMILY_OVERLAP_FILE_NAME, FamilyOverlapAccumulator}; +pub use family_overlap::FamilyOverlap; pub use iter::{ MinorantBatchIter, MinorantIter, SiblingBatchIter, SiblingEntry, SiblingIter, SiblingLayerExt, };