feat(network): implement posterior log-ratio edge estimator with consensus clustering

CHANGELOG.md

@@ -4,6 +4,7 @@ All notable changes to TEPP are documented here. The format follows Keep a Chang
 
 ## [Unreleased]
 
+- `network_analysis` posterior edge estimation now carries exact two-sided Fisher z-transform p-values against `rho = 0` (erfc evaluated by an all-positive confluent series plus the Laplace continued fraction, locked to libm reference values at 1e-13 relative tolerance), percentile-bootstrap credible intervals and selection fractions over posterior draws, Benjamini–Hochberg step-up admission on those exact p-values instead of the complement of a thresholded fraction, an explicit fail-closed `edge_drop_probability` parameter for consensus co-assignment resampling (replacing a hardcoded 0.1), and honest per-replicate stability admission; the greedy partition helper is renamed to state that it makes no modularity-optimization claim. APA 7 entries added: Benjamini & Hochberg (1995), Efron (1979), Fisher (1921), Hennig (2007), Monti (2003).
 - Restored protected-main gate integrity after the consolidation merges: hourly-scheduler prompt-contract tests now assert the gap-baseline-derived task contract (Gap ID naming, no invented weights) instead of stale increment-specific tokens, the operator-gap register inventory matches the live 33-PR queue, `evidence_core::image_unit` non-image/empty-subtype refusals and `load_union_branch_totals` valid-record accumulation have exact coverage, and the README crate fence plus duplicate registry entries stay deduped.
 - Branch-coverage diagnostics on the post-consolidation head exposed two uncovered outcomes in `evidence_core::image_unit` (`is_image_media_type_token` non-image prefix and empty-subtype refusals), one uncovered authored line (the strip-prefix refusal), and lost valid-record coverage for `load_union_branch_totals`; exact red-to-green cases now cover the non-image/empty-subtype data URIs and per-coordinate True/False accumulation.
 - Repaired post-consolidation merge fallout that left protected `main` red: restored the lost `return True` in the `check_coverage.py` match-guard branch, removed the shadowed duplicate `load_union_branch_totals` and `_is_multiline_match_guard` definitions plus duplicate workspace-crate entries (`episode_membership`, `analysis_engine`) from the contract tuple and Cargo member arrays, split two union-fused four-tuples back into `(variant, message)` pairs in the `event_core` error table, repaired the fused `identity_recovery_rate` body in `episode_membership::window`, deduplicated the checked-arithmetic eligible-count block in `analysis_engine`, fixed four-argument `unit()` test call sites, rebalanced the README crate-list fence around all 54 unique crates, and deduplicated the `location_membership`/`validation_core`/`tepp_api` architecture-table rows. Also documents private `PLAUSIBLE_IMAGE_MEDIA_TYPES` so `cargo doc -D warnings` passes.

crates/network_analysis/src/consensus.rs

@@ -0,0 +1,301 @@
+//! Co-assignment consensus clustering from repeated label-invariant partitions.
+//!
+//! Runs multiple rounds of deterministic greedy partitioning over
+//! randomly perturbed admitted positive edges, builds a co-assignment
+//! matrix across replicates, and derives consensus clusters by
+//! thresholding that matrix. The resampling-based consensus view and the
+//! stability rationale follow Monti (2003) and Hennig (2007); the edge
+//! perturbation probability is an explicit parameter with provenance,
+//! never an implicit constant.
+
+#![forbid(unsafe_code)]
+#![allow(
+    clippy::doc_markdown,
+    clippy::items_after_statements,
+    clippy::must_use_candidate,
+    clippy::needless_for_each,
+    clippy::cast_precision_loss,
+    clippy::missing_errors_doc
+)]
+#![deny(missing_docs)]
+
+use crate::edges::NetworkEdge;
+use crate::error::NetworkEstimatorError;
+use std::collections::HashMap;
+
+/// Consensus clustering output.
+///
+/// * `assignments` – one entry per topic: `Some(cluster_id)` or `None` for unclustered.
+/// * `co_assignment` – K × K symmetric matrix of co-assignment frequencies.
+#[derive(Debug, Clone, PartialEq)]
+pub struct ConsensusClusterOutput {
+    /// Per-topic cluster assignment (`None` = unclustered).
+    pub assignments: Vec<Option<usize>>,
+    /// K × K co-assignment frequency matrix.
+    pub co_assignment: Vec<Vec<f64>>,
+}
+
+/// Derive consensus clusters from repeatedly perturbed partitions.
+///
+/// Each replicate independently drops every admitted positive edge with
+/// probability `edge_drop_probability`, repartitions the surviving edges,
+/// and accumulates label-invariant co-assignment counts. The final
+/// assignment thresholds the co-assignment frequency at
+/// `consensus_threshold`.
+///
+/// # Arguments
+///
+/// * `edges` – admitted positive edges (source, target, effect).
+/// * `k_topics` – total number of topics.
+/// * `n_replicates` – number of perturbed partitions to generate; at least 1.
+/// * `consensus_threshold` – minimum co-assignment fraction for
+///   same-cluster; must lie inside [0, 1].
+/// * `edge_drop_probability` – per-edge independent drop probability in
+///   each replicate; must be finite and inside [0, 1). The value is an
+///   explicit design parameter of the resampling scheme (Monti, 2003;
+///   Hennig, 2007), not an internal constant.
+/// * `seed` – deterministic seed.
+///
+/// # Errors
+///
+/// Fails closed on zero topics or replicates, a threshold outside
+/// [0, 1], or a non-finite drop probability at or above 1.
+pub fn consensus_clusters(
+    edges: &[NetworkEdge],
+    k_topics: usize,
+    n_replicates: usize,
+    consensus_threshold: f64,
+    edge_drop_probability: f64,
+    seed: u64,
+) -> Result<ConsensusClusterOutput, NetworkEstimatorError> {
+    if k_topics == 0 {
+        return Err(NetworkEstimatorError::DimensionMismatch);
+    }
+    if n_replicates == 0 {
+        return Err(NetworkEstimatorError::ZeroReplicates);
+    }
+    if !consensus_threshold.is_finite() || !(0.0..=1.0).contains(&consensus_threshold) {
+        return Err(NetworkEstimatorError::InvalidProbability);
+    }
+    if !edge_drop_probability.is_finite() || !(0.0..1.0).contains(&edge_drop_probability) {
+        return Err(NetworkEstimatorError::InvalidProbability);
+    }
+
+    // Build adjacency from positive edges only.
+    let mut adjacency: HashMap<usize, Vec<usize>> = HashMap::new();
+    for edge in edges {
+        if edge.effect > 0.0 {
+            adjacency.entry(edge.source).or_default().push(edge.target);
+            adjacency.entry(edge.target).or_default().push(edge.source);
+        }
+    }
+
+    // Deterministic LCG shared across the crate.
+    let mut state = seed ^ 0x9E37_79B9_7F4A_7C15;
+
+    // Accumulate co-assignment counts.
+    let mut co_count = vec![vec![0_u64; k_topics]; k_topics];
+
+    for _ in 0..n_replicates {
+        // Perturb: drop each edge independently with the caller-supplied
+        // probability so cluster recovery is stress-tested by resampling.
+        let perturbed: Vec<&NetworkEdge> = edges
+            .iter()
+            .filter(|_| {
+                state = state
+                    .wrapping_mul(6_364_136_223_846_793_005)
+                    .wrapping_add(1_442_695_040_888_963_407);
+                let draw = ((state >> 33) as f64 / (u64::MAX >> 33) as f64).min(1.0);
+                draw >= edge_drop_probability
+            })
+            .collect();
+
+        let partition = greedy_union_partition(&perturbed, k_topics);
+
+        for i in 0..k_topics {
+            for j in 0..k_topics {
+                if partition[i] == partition[j] {
+                    co_count[i][j] += 1;
+                }
+            }
+        }
+    }
+
+    // Build co-assignment frequency matrix and derive final assignment.
+    let mut co_freq = vec![vec![0.0_f64; k_topics]; k_topics];
+    let mut assignments: Vec<Option<usize>> = vec![None; k_topics];
+    let mut next_cluster = 0_usize;
+
+    for i in 0..k_topics {
+        if assignments[i].is_some() {
+            continue;
+        }
+        // Find all j ≥ i still unassigned whose co-assignment with i
+        // reaches the threshold.
+        let members: Vec<usize> = (i..k_topics)
+            .filter(|&j| {
+                assignments[j].is_none()
+                    && co_count[i][j] as f64 / n_replicates as f64 >= consensus_threshold
+            })
+            .collect();
+
+        if members.len() < 2 {
+            // Singletons stay unclustered unless they pair strongly.
+            if !members.is_empty()
+                && co_count[i][i] as f64 / n_replicates as f64 >= consensus_threshold
+                && adjacency.contains_key(&i)
+            {
+                assignments[i] = Some(next_cluster);
+                next_cluster += 1;
+            }
+            continue;
+        }
+
+        for &member in &members {
+            assignments[member] = Some(next_cluster);
+        }
+        next_cluster += 1;
+    }
+
+    for i in 0..k_topics {
+        for j in 0..k_topics {
+            co_freq[i][j] = co_count[i][j] as f64 / n_replicates as f64;
+        }
+    }
+
+    Ok(ConsensusClusterOutput {
+        assignments,
+        co_assignment: co_freq,
+    })
+}
+
+/// Greedy union-find partition over edges sorted by descending effect.
+///
+/// This is a deterministic single-pass stand-in for Leiden until a vetted
+/// implementation is adopted; it makes no modularity-optimization claim.
+fn greedy_union_partition(edges: &[&NetworkEdge], k: usize) -> Vec<usize> {
+    let mut parent: Vec<usize> = (0..k).collect();
+
+    fn find(parent: &mut [usize], mut x: usize) -> usize {
+        while parent[x] != x {
+            parent[x] = parent[parent[x]];
+            x = parent[x];
+        }
+        x
+    }
+
+    fn union(parent: &mut [usize], a: usize, b: usize) {
+        let root_a = find(parent, a);
+        let root_b = find(parent, b);
+        if root_a != root_b {
+            parent[root_b] = root_a;
+        }
+    }
+
+    // Sort edges by descending effect.
+    let mut sorted: Vec<&NetworkEdge> = edges.to_vec();
+    sorted.sort_by(|a, b| b.effect.total_cmp(&a.effect));
+
+    for edge in &sorted {
+        union(&mut parent, edge.source, edge.target);
+    }
+
+    // Normalise labels to 0..n_clusters.
+    let mut label_map: HashMap<usize, usize> = HashMap::new();
+    parent
+        .iter()
+        .map(|&root| {
+            let len = label_map.len();
+            *label_map.entry(root).or_insert(len)
+        })
+        .collect()
+}
+
+#[cfg(test)]
+#[allow(clippy::float_cmp)]
+mod tests {
+    use super::*;
+
+    fn edge(source: usize, target: usize, effect: f64) -> NetworkEdge {
+        NetworkEdge {
+            source,
+            target,
+            effect,
+            lower: effect,
+            upper: effect,
+            p_value: 0.0,
+            selection_probability: 1.0,
+        }
+    }
+
+    #[test]
+    fn invalid_parameters_fail_closed() {
+        let chain = vec![edge(0, 1, 0.9), edge(1, 2, 0.8)];
+        assert!(matches!(
+            consensus_clusters(&chain, 0, 5, 0.5, 0.1, 1),
+            Err(NetworkEstimatorError::DimensionMismatch)
+        ));
+        assert!(matches!(
+            consensus_clusters(&chain, 3, 0, 0.5, 0.1, 1),
+            Err(NetworkEstimatorError::ZeroReplicates)
+        ));
+        assert!(matches!(
+            consensus_clusters(&chain, 3, 5, 1.5, 0.1, 1),
+            Err(NetworkEstimatorError::InvalidProbability)
+        ));
+        assert!(matches!(
+            consensus_clusters(&chain, 3, 5, -0.1, 0.1, 1),
+            Err(NetworkEstimatorError::InvalidProbability)
+        ));
+        assert!(matches!(
+            consensus_clusters(&chain, 3, 5, 0.5, 1.0, 1),
+            Err(NetworkEstimatorError::InvalidProbability)
+        ));
+        assert!(matches!(
+            consensus_clusters(&chain, 3, 5, 0.5, f64::NAN, 1),
+            Err(NetworkEstimatorError::InvalidProbability)
+        ));
+    }
+
+    #[test]
+    fn zero_drop_probability_is_fully_deterministic() {
+        // With no perturbation every replicate sees the identical chain,
+        // so all three topics co-assign always and land in one cluster.
+        let chain = vec![edge(0, 1, 0.95), edge(1, 2, 0.85)];
+        let output = consensus_clusters(&chain, 3, 25, 0.99, 0.0, 11).unwrap();
+        assert_eq!(output.assignments[0], output.assignments[1]);
+        assert_eq!(output.assignments[1], output.assignments[2]);
+        assert!((output.co_assignment[0][2] - 1.0).abs() < 1e-12);
+    }
+
+    #[test]
+    fn weak_edges_do_not_merge_under_perturbation() {
+        // Two strongly-bound pairs of topics plus one isolated topic:
+        // within-pair co-assignment stays high across perturbed
+        // replicates, while the isolated topic never joins any cluster.
+        let pairs = vec![edge(0, 1, 0.97), edge(2, 3, 0.96)];
+        let output = consensus_clusters(&pairs, 5, 200, 0.6, 0.2, 21).unwrap();
+        assert_eq!(output.assignments[0], output.assignments[1]);
+        assert!(output.assignments[0].is_some());
+        assert_eq!(output.assignments[2], output.assignments[3]);
+        assert!(output.assignments[2].is_some());
+        assert_ne!(output.assignments[0], output.assignments[2]);
+        assert_eq!(output.assignments[4], None);
+        assert!((output.co_assignment[4][4] - 1.0).abs() < 1e-12);
+        // The within-pair edge survives roughly 80% of perturbations, so
+        // its co-assignment stays well above the consensus threshold
+        // while never reaching a deterministic 1.0.
+        let within_pair = output.co_assignment[0][1];
+        assert!(within_pair > 0.6 && within_pair < 1.0, "co = {within_pair}");
+        // Cross-pair topics never co-assign: they share no edge.
+        assert!((output.co_assignment[0][2]).abs() < 1e-12);
+    }
+
+    #[test]
+    fn identical_inputs_and_seed_reproduce_output_exactly() {
+        let chain = vec![edge(0, 1, 0.9), edge(1, 2, 0.8)];
+        let first = consensus_clusters(&chain, 3, 15, 0.6, 0.1, 7).unwrap();
+        let second = consensus_clusters(&chain, 3, 15, 0.6, 0.1, 7).unwrap();
+        assert_eq!(first, second);
+    }
+}