Push zunrplorkwkt #70

Merged
coissac merged 93 commits from push-zunrplorkwkt into main 2026-08-28 23:15:38 +00:00
10 changed files with 388 additions and 59 deletions
Showing only changes of commit e101f629e6 - Show all commits
+31 -4
View File
@@ -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<String>,
/// 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<f64>,
/// 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
+61 -3
View File
@@ -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<bool> = 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,
)
@@ -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<usize> = 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(())
}
@@ -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<i}(n_genomes-1-r) + (j-i-1)`.
/// See `crate::siblings::helpers::triangle_index`.
fn flat_index(&self, i: usize, j: usize) -> 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]
@@ -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,
@@ -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(())
}
}
@@ -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<EntropyBias>,
gamma_shape: Option<f64>,
) -> OKIResult<Array2<f64>>;
/// 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<FamilyOverlap>;
}
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<Array2<f64>> {
snp_distance(self, kind, n, free_loss, no_ambiguity, excluded, entropy_bias, gamma_shape)
}
fn family_overlap(&self) -> OKIResult<FamilyOverlap> {
self.ensure_entropy_and_overlap_annexes()?;
let path = self.index().dir().join(FAMILY_OVERLAP_FILE_NAME);
FamilyOverlap::open(&path).map_err(OKIError::Io)
}
}
@@ -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<Self> {
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<u64>,
}
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);
}
}
+14
View File
@@ -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<i}(n-1-r) + (j-i-1)`. Shared
/// by every per-genome-pair structure that only stores the upper triangle
/// (`PairwiseTally`, `FamilyOverlap`) rather than each re-deriving the same
/// formula.
#[inline]
pub(super) fn triangle_index(n: usize, i: usize, j: usize) -> 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)
}
+3
View File
@@ -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,
};