From 89629e118db62138921b462045810cb20b35f2be Mon Sep 17 00:00:00 2001 From: Eric Coissac Date: Thu, 27 Aug 2026 15:30:17 +0200 Subject: [PATCH] Add sibling family size distribution statistics and CLI flags Introduces `--sibling-stats` and `--sibling-hist` flags to compute and export family-size distributions from pre-built annexes. The new algorithms module implements parallel and sequential scanning routines, while the CLI layer handles CSV export with deduplication and standardized error handling. --- src/obikmer2/src/cmd/phylo/args.rs | 24 +++- src/obikmer2/src/cmd/phylo/mod.rs | 60 ++++++++ src/obikphylo/src/siblings/algorithms/mod.rs | 3 + .../src/siblings/algorithms/stats.rs | 133 ++++++++++++++++++ .../src/siblings/extensions/sibling_ext.rs | 30 +++- src/obikphylo/src/siblings/mod.rs | 2 +- 6 files changed, 245 insertions(+), 7 deletions(-) create mode 100644 src/obikphylo/src/siblings/algorithms/stats.rs diff --git a/src/obikmer2/src/cmd/phylo/args.rs b/src/obikmer2/src/cmd/phylo/args.rs index dcf9a7d5..85a64a2c 100644 --- a/src/obikmer2/src/cmd/phylo/args.rs +++ b/src/obikmer2/src/cmd/phylo/args.rs @@ -37,12 +37,13 @@ impl From for DistanceMetric { /// Partial transfer of `obikmer`'s `phylo` command: the plain distance-metric /// path (`--metric`/NJ/UPGMA), annex construction (`--sibling-annex`), -/// entropy reporting (`--shannon`) and SNP pseudo-alignment sampling +/// annex diagnostics (`--sibling-stats`, `--sibling-hist`), entropy +/// reporting (`--shannon`) and SNP pseudo-alignment sampling /// (`--pseudo-alignment`, `--subsample`, `--free-loss`, `--no-ambiguity`, /// `--entropy`/`--entropy-sd`) — everything else sibling-annex-based -/// (`--sankoff`, `--tnt`/`--phyg`/`--iqtree`, stats, ...) stays in `obikmer` -/// until the rest of `obikphylo::siblings` is reconnected (see the project -/// memory on this). +/// (`--sankoff`, `--tnt`/`--phyg`/`--iqtree`, raw SNP distance, family +/// overlap, ...) stays in `obikmer` until the rest of +/// `obikphylo::siblings` is reconnected (see the project memory on this). #[derive(Args)] pub struct PhyloArgs { /// Index directory @@ -64,6 +65,21 @@ pub struct PhyloArgs { #[arg(long)] pub sibling_annex: bool, + /// Tally the sibling-count distribution (CSV) of an already-built annex + /// (run with `--sibling-annex` first, in this invocation or an earlier + /// one). A separate, occasional diagnostic pass — not run every time the + /// annex itself is (re)built. + #[arg(long)] + pub sibling_stats: bool, + + /// Print just the global family-size histogram (1-4 members) of an + /// already-built annex — the `global` row `--sibling-stats` also + /// writes, but without the per-genome breakdown, so it skips + /// `--sibling-stats`'s cross-partition resolution entirely (annex bits + /// only). + #[arg(long)] + pub sibling_hist: 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 9e0958fb..7e275d47 100644 --- a/src/obikmer2/src/cmd/phylo/mod.rs +++ b/src/obikmer2/src/cmd/phylo/mod.rs @@ -77,6 +77,66 @@ pub fn run(args: PhyloArgs) { rep.push(t.stop()); } + // ── Sibling-count distribution (`--sibling-stats`) ────────────────────────── + if args.sibling_stats { + let t = Stage::start("sibling_stats"); + let stats = cache.sibling_annex_stats().unwrap_or_else(|e| { + eprintln!("error computing sibling-annex stats: {e}"); + std::process::exit(1); + }); + rep.push(t.stop()); + + let path = args.output.as_ref() + .map(|p| format!("{}_siblings.csv", p.display())) + .unwrap_or_else(|| "siblings.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); + })); + // One row per genome (4 columns, family size 1-4: number of + // families of that size for which the genome carries at least one + // member), plus a `global` row — the actual deduplicated + // family-size histogram (`stats.counts`), NOT a sum of the + // per-genome columns (a family shared by several genomes would + // otherwise be counted once per genome it appears in). + writeln!(f, "genome,1,2,3,4").unwrap(); + for (label, counts) in labels.iter().zip(stats.per_genome.iter()) { + writeln!(f, "{label},{},{},{},{}", counts[0], counts[1], counts[2], counts[3]).unwrap(); + } + writeln!( + f, "global,{},{},{},{}", + stats.counts[0], stats.counts[1], stats.counts[2], stats.counts[3], + ).unwrap(); + info!("sibling-count distribution → {path}"); + } + + // ── Family-size histogram (`--sibling-hist`) ──────────────────────────────── + if args.sibling_hist { + let t = Stage::start("sibling_hist"); + let counts = cache.sibling_family_size_histogram().unwrap_or_else(|e| { + eprintln!("error computing sibling family-size histogram: {e}"); + std::process::exit(1); + }); + rep.push(t.stop()); + + let path = args.output.as_ref() + .map(|p| format!("{}_sibling_hist.csv", p.display())) + .unwrap_or_else(|| "sibling_hist.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); + })); + writeln!(f, "size,count").unwrap(); + for (size, count) in counts.iter().enumerate() { + writeln!(f, "{},{count}", size + 1).unwrap(); + } + let total: u64 = counts.iter().sum(); + info!( + "family-size histogram → {path} (total {total} famil{})", + if total == 1 { "y" } else { "ies" } + ); + } + // ── Shannon entropy report (`--shannon`) ──────────────────────────────────── if args.shannon { let path = args.output.as_ref() diff --git a/src/obikphylo/src/siblings/algorithms/mod.rs b/src/obikphylo/src/siblings/algorithms/mod.rs index 416a00c2..9c0d95d5 100644 --- a/src/obikphylo/src/siblings/algorithms/mod.rs +++ b/src/obikphylo/src/siblings/algorithms/mod.rs @@ -19,17 +19,20 @@ mod entropy; mod family_scan; mod masking; mod minorant_selection; +mod stats; mod subsample; use obikidxcache::index_cache::IndexCache; pub use alignment::SnpAlignment; +pub use stats::SiblingAnnexStats; pub use subsample::EntropyBias; pub(crate) use alignment::snp_pseudo_alignment; pub(crate) use annex::build_layer_sibling_annex; pub(crate) use entropy::{ensure_layer_entropy_annex, family_entropy, family_entropy_4, iter_full_entropy}; pub(crate) use family_scan::{Selection, scan_layer_families}; +pub(crate) use stats::{sibling_annex_stats, sibling_family_size_histogram}; pub(crate) use subsample::sample_index; /// Whether every layer number in `cache` fits in a `FamilyMask` field diff --git a/src/obikphylo/src/siblings/algorithms/stats.rs b/src/obikphylo/src/siblings/algorithms/stats.rs new file mode 100644 index 00000000..e8417ef8 --- /dev/null +++ b/src/obikphylo/src/siblings/algorithms/stats.rs @@ -0,0 +1,133 @@ +//! Sibling-annex diagnostics: family-size (1-4 members) distribution, read +//! back from an already-built annex — a separate, occasional pass, never +//! fused into [`super::build_layer_sibling_annex`] itself. +//! +//! Every count here is **per family, not per slot**: a family with `F` +//! members occupies `F` annex slots (one per observed member), all sharing +//! the same mask. Counting every slot would count each family up to 4 times +//! over; only the minorant's slot is tallied (`FamilyMask::is_minorant()`, +//! already stored — no re-derivation, no lookup). + +use obikidxcache::index_cache::IndexCache; +use obikindex::{OKIError, OKIResult}; + +use super::family_scan::{Selection, scan_layer_families}; +use super::is_fast_mode; +use crate::siblings::{ANNEX_FILE_NAME, SiblingAnnex}; + +/// Layers scanned at once per [`sibling_family_size_histogram`]'s +/// `par_map_reduce` call — same bounding reasoning as +/// `distance::MAX_CONCURRENT_LAYERS`, a separate constant because the two +/// live in different modules with no shared dependency to hang a common one +/// off. +const MAX_CONCURRENT_LAYERS: usize = 8; + +/// `counts[s]` = number of *families* of size `s + 1` (index 0 = size 1, +/// i.e. 0 siblings, ... index 3 = size 4, i.e. 3 siblings). `per_genome[g][s]` +/// = number of families of size `s + 1` for which genome `g` (index into +/// `IndexCache::meta().genomes()`) carries at least one member. `pub`: part +/// of the public signature of +/// [`crate::siblings::extensions::SiblingExt::sibling_annex_stats`]. +#[derive(Debug, Clone, Default)] +pub struct SiblingAnnexStats { + pub counts: [u64; 4], + pub per_genome: Vec<[u64; 4]>, +} + +/// Every cached layer's own annex file is independent — no cross-partition +/// resolution, no shared state between layers — so this is a genuine +/// `IndexCache::par_map_reduce`: each layer's own `[u64; 4]` partial (annex +/// bits only, `FamilyMask::is_minorant()` + `siblings()`) summed across +/// every cached layer. Unlike [`sibling_annex_stats`], which drives +/// `scan_layer_families`'s own internal partition-grouped parallelism and +/// must stay a plain sequential loop across layers (see +/// `IndexCache::par_map_reduce`'s own docs on why stacking parallel layers +/// on top of that would fight it for cache/mmap locality) — this function +/// has no such inner parallelism to protect. +/// +/// Existence of every cached layer's annex file is checked up front, before +/// the parallel map runs, so a missing one is reported as an error rather +/// than silently read back as zero. +pub(crate) fn sibling_family_size_histogram(cache: &IndexCache) -> OKIResult<[u64; 4]> { + for (part, l, layer) in cache.iter_indexed() { + let path = layer.dir().join(ANNEX_FILE_NAME); + if !path.exists() { + return Err(OKIError::InvalidInput(format!( + "no sibling annex at {} (partition {part}, layer {l}) — run --sibling-annex first", + path.display() + ))); + } + } + + Ok(cache.par_map_reduce( + MAX_CONCURRENT_LAYERS, + || [0u64; 4], + |layer| { + let annex = SiblingAnnex::open(&layer.dir().join(ANNEX_FILE_NAME)) + .expect("existence already checked before this parallel map ran"); + let mut counts = [0u64; 4]; + for slot in 0..annex.len() { + let Some(mask) = annex.get(slot) else { continue }; + if !mask.is_minorant() { + continue; // family tallied once, at its minorant + } + counts[mask.siblings() as usize] += 1; + } + counts + }, + |a, b| { + let mut sum = a; + for i in 0..4 { + sum[i] += b[i]; + } + sum + }, + )) +} + +/// Tally the family-size distribution of an already-built annex (globally, +/// and per genome) — needs each family's per-genome carry, so unlike +/// [`sibling_family_size_histogram`] it pays `scan_layer_families`'s full +/// cross-partition resolution cost, one layer at a time (not +/// `par_map_reduce`: `scan_layer_families` already parallelises internally +/// at a finer grain, across partitions — stacking another parallel layer +/// loop on top would fight it for locality instead of helping, exactly the +/// same reasoning `SiblingExt::build_sibling_annex`/`shannon_entropy_csv`/ +/// `snp_pseudo_alignment` already follow). +pub(crate) fn sibling_annex_stats(cache: &IndexCache) -> OKIResult { + let n_genomes = cache.meta().genomes().len(); + let fast_mode = is_fast_mode(cache); + + let mut stats = SiblingAnnexStats { + per_genome: vec![[0u64; 4]; n_genomes], + ..Default::default() + }; + + for part in cache.partitions() { + let n_layer = cache.n_layer(part).unwrap_or(0); + for l in 0..n_layer { + scan_layer_families( + cache, + part, + l, + n_genomes, + fast_mode, + &Selection::All, + |_family_idx, mask, genome_mask| { + // "Genome g represents this family" means g carries + // *any* of its members, not just the minorant's own — + // `genome_mask[g] != 0` is exactly that. + let s = mask.siblings() as usize; + stats.counts[s] += 1; + for (g, &m) in genome_mask.iter().enumerate() { + if m != 0 { + stats.per_genome[g][s] += 1; + } + } + }, + )?; + } + } + + Ok(stats) +} diff --git a/src/obikphylo/src/siblings/extensions/sibling_ext.rs b/src/obikphylo/src/siblings/extensions/sibling_ext.rs index bf335a99..ace56266 100644 --- a/src/obikphylo/src/siblings/extensions/sibling_ext.rs +++ b/src/obikphylo/src/siblings/extensions/sibling_ext.rs @@ -14,8 +14,9 @@ use obikindex::{OKIError, OKIResult}; use obisys::progress_bar; use crate::siblings::algorithms::{ - EntropyBias, Selection, SnpAlignment, build_layer_sibling_annex, family_entropy, - family_entropy_4, is_fast_mode, scan_layer_families, snp_pseudo_alignment, + EntropyBias, Selection, SiblingAnnexStats, SnpAlignment, build_layer_sibling_annex, + family_entropy, family_entropy_4, is_fast_mode, scan_layer_families, + sibling_annex_stats, sibling_family_size_histogram, snp_pseudo_alignment, }; use crate::siblings::extensions::SiblingBuilder; use crate::siblings::ENTROPY_ANNEX_FILE_NAME; @@ -76,6 +77,23 @@ pub trait SiblingExt { excluded: &[bool], entropy_bias: Option, ) -> OKIResult; + + /// The global family-size histogram alone (`SiblingAnnexStats::counts`, + /// no `per_genome`) — reads only the already-built annex + /// (`FamilyMask::siblings()` + `is_minorant()`, already-stored bits, + /// mmap'd), nothing else: no cross-partition resolution, no + /// `scan_layer_families`. A genuine `IndexCache::par_map_reduce` (see + /// its own docs) — every cached layer's annex file is independent, so + /// layers scan concurrently, unlike [`sibling_annex_stats`](Self::sibling_annex_stats) + /// below. + fn sibling_family_size_histogram(&self) -> OKIResult<[u64; 4]>; + + /// Tally the family-size distribution of an already-built annex + /// (globally, and per genome) — pays the full per-genome + /// cross-partition resolution cost to get the per-genome breakdown; + /// use [`sibling_family_size_histogram`](Self::sibling_family_size_histogram) + /// instead when only the global histogram is needed. + fn sibling_annex_stats(&self) -> OKIResult; } impl SiblingExt for IndexCache { @@ -200,4 +218,12 @@ impl SiblingExt for IndexCache { ) -> OKIResult { snp_pseudo_alignment(self, n, free_loss, no_ambiguity, excluded, entropy_bias) } + + fn sibling_family_size_histogram(&self) -> OKIResult<[u64; 4]> { + sibling_family_size_histogram(self) + } + + fn sibling_annex_stats(&self) -> OKIResult { + sibling_annex_stats(self) + } } diff --git a/src/obikphylo/src/siblings/mod.rs b/src/obikphylo/src/siblings/mod.rs index 3dcdb2c3..00dbf759 100644 --- a/src/obikphylo/src/siblings/mod.rs +++ b/src/obikphylo/src/siblings/mod.rs @@ -29,7 +29,7 @@ mod siblingannex; pub(crate) use entropy_annex::{ENTROPY_ANNEX_FILE_NAME, EntropyAnnex, EntropyAnnexBuilder}; pub(crate) use siblingannex::{FamilyMask, SiblingAnnex, SiblingAnnexBuilder}; -pub use algorithms::{EntropyBias, SnpAlignment}; +pub use algorithms::{EntropyBias, SiblingAnnexStats, SnpAlignment}; pub use extensions::SiblingExt; pub(crate) const ANNEX_FILE_NAME: &str = "siblings.psib";