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

@@ -37,12 +37,13 @@ impl From<MetricArg> 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`
+/// (`--sankoff`, `--tnt`/`--phyg`/`--iqtree`, raw SNP distance, family
-/// until the rest of `obikphylo::siblings` is reconnected (see the project
+/// overlap, ...) stays in `obikmer` until the rest of
-/// memory on this).
+/// `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

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()

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

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<SiblingAnnexStats> {
+    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)
+}

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,
+    EntropyBias, Selection, SiblingAnnexStats, SnpAlignment, build_layer_sibling_annex,
-    family_entropy_4, is_fast_mode, scan_layer_families, snp_pseudo_alignment,
+    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<EntropyBias>,
     ) -> OKIResult<SnpAlignment>;
+
+    /// 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<SiblingAnnexStats>;
 }
 
 impl SiblingExt for IndexCache {
@@ -200,4 +218,12 @@ impl SiblingExt for IndexCache {
     ) -> OKIResult<SnpAlignment> {
         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<SiblingAnnexStats> {
+        sibling_annex_stats(self)
+    }
 }

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";