diff --git a/CHANGELOG.md b/CHANGELOG.md index 062a69412..c2681bf10 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -38,6 +38,8 @@ All notable changes to TEPP are documented here. The format follows Keep a Chang ## [Unreleased] +- `network_analysis` replaces the greedy union-find consensus stand-in with Traag, Waltman, and van Eck (2019) Leiden modularity (γ = 1): fast local moving, refinement that keeps communities internally connected, and aggregation. Two cliques joined by a weak bridge stay two communities; union-find glued them. Isolated topics stay unclustered. Co-assignment follows Monti, Tamayo, Mesirov, and Golub (2003) and Hennig (2007) with an explicit drop probability; the drop stream is independent of the per-replicate Leiden RNG, admitted edges are sorted, and out-of-range endpoints stay unclustered. This is not a graphical lasso, not a causal cluster, and not an export workflow. + - `event_core` adds bounded Allen interval-consistency classification, atomic path-consistency closure, contradiction/resource refusals, and an explicit dependency-error fallback without claiming unrestricted global satisfiability. - `psychometric_core` recovers the Driver, Oud, and Voelkle (2017, Table 2, p. 12 `MANIFESTTRAITVAR`; §7.1, p. 19; p. 16 `MANIFESTTRAITVARstd`; footnote 4; 2017-era ctsem `summary.ctsemFit.R`; JSS PDF re-opened 2026-08-27T14:20Z from https://www.jstatsoft.org/index.php/jss/article/download/v077i05/1104) scalar standardised manifest-trait variance on current main after `0ce16e8` dropped the pre-consolidation code while research notes already named the map (register items 83–84). Table 2 names `MANIFESTTRAITVAR` `Ψ_τ` the additional time-invariant variance-covariance on the measurement level and sets it `NULL` when there is no manifest trait. Equation 5 writes `Γ ~ N(τ, Ψ)` and names that covariance the manifest traits. Section 7.1 names manifest traits stable individual differences in indicator levels, distinct from process-level `TRAITVAR` `φ_ξ`. Page 16 prints standardised matrices with the suffix `std` when appropriate. The printed example on p. 16 is `discreteDRIFTstd`, not `MANIFESTTRAITVARstd`. Footnote 4 standardises using only the relevant variance, not the total. The relevant variance for that named indicator-level correlation is `MANIFESTTRAITVAR`, not process-level `TRAITVAR` and not residual `MANIFESTVAR` `θ`. The 2017-era source forms `MANIFESTTRAITVARstd` only when `MANIFESTTRAITVAR != 0`, as `solve(sqrt(diag(MANIFESTTRAITVAR) + ridging)) %&% MANIFESTTRAITVAR` when `verbose = TRUE`. OpenMx `%&%` is `t(A) %*% B %*% A`. Unlike `TRAITVARstd`, that formation adds `diag(c(ridging), n.manifest)`. The default `ridging = FALSE` adds 0, not `0.0001`; that ridge is a numerical hack and is not this exact map. The scalar correlation is `ψ / ψ = 1` after strictly positive `MANIFESTTRAITVAR`. Form strictly positive `ψ` first, then `1 / √ψ`, then `(1 / √ψ) ψ (1 / √ψ)`. Unstandardised `MANIFESTTRAITVAR` is defined for a zero trait; standardised `MANIFESTTRAITVAR` is not. Zero `MANIFESTTRAITVAR` skips forming `MANIFESTTRAITVARstd` in the 2017-era source and fails closed here. Indicator-level trait variance is an event-time structural quantity, so a non-event clock fails closed. `MANIFESTTRAITVAR` does not require stable `a < 0`. Distinct positive `ψ` recover the same 1. `trait / trait = 1` is `TRAITVARstd` and recovers the same number and remains a distinct named quantity. `θ` is `MANIFESTVAR` and is measurement error, not this correlation. Meredith (1993) remains unread (web search 2026-08-27T14:20Z: Springer/Cambridge Core paywalled; Unpaywall historically `is_oa: false`; Springer `content/pdf` is an HTML stub). Mislevy (1991, *Psychometrika, 56*, 177–196) remains unread on the same terms (DOI `10.1007/bf02294457`). Still not a Kalman filter, not a matrix `expm`, not ESEM estimation, not DSEM, and not ctsem estimation. diff --git a/crates/network_analysis/src/consensus.rs b/crates/network_analysis/src/consensus.rs index f18efa9d7..00b903128 100644 --- a/crates/network_analysis/src/consensus.rs +++ b/crates/network_analysis/src/consensus.rs @@ -1,12 +1,13 @@ -//! Co-assignment consensus clustering from repeated label-invariant partitions. +//! Co-assignment consensus clustering from repeated Leiden 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 +//! Runs multiple rounds of Traag, Waltman, and van Eck (2019) Leiden +//! modularity on 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, Tamayo, Mesirov, and Golub (2003) +//! and Hennig (2007); the edge //! perturbation probability is an explicit parameter with provenance, -//! never an implicit constant. +//! never an implicit constant. Union-find is not used. #![forbid(unsafe_code)] #![allow( @@ -21,6 +22,7 @@ use crate::edges::NetworkEdge; use crate::error::NetworkEstimatorError; +use crate::leiden::leiden_partition; use std::collections::HashMap; /// Consensus clustering output. @@ -35,13 +37,13 @@ pub struct ConsensusClusterOutput { pub co_assignment: Vec>, } -/// Derive consensus clusters from repeatedly perturbed partitions. +/// Derive consensus clusters from repeatedly perturbed Leiden 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`. +/// probability `edge_drop_probability`, repartitions the surviving edges +/// with Leiden modularity (Traag et al., 2019), and accumulates +/// label-invariant co-assignment counts. The final assignment thresholds +/// the co-assignment frequency at `consensus_threshold`. /// /// # Arguments /// @@ -52,8 +54,8 @@ pub struct ConsensusClusterOutput { /// 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. +/// explicit design parameter of the resampling scheme (Monti, Tamayo, +/// Mesirov, & Golub, 2003; Hennig, 2007), not an internal constant. /// * `seed` – deterministic seed. /// /// # Errors @@ -81,39 +83,51 @@ pub fn consensus_clusters( return Err(NetworkEstimatorError::InvalidProbability); } - // Build adjacency from positive edges only. + // Admitted positive in-range edges only, sorted so input permutation + // cannot change the drop stream. Out-of-range and self-loop endpoints + // never enter adjacency or Leiden (they are not topics in this k). + let mut admitted: Vec<&NetworkEdge> = edges + .iter() + .filter(|edge| { + edge.effect > 0.0 + && edge.source < k_topics + && edge.target < k_topics + && edge.source != edge.target + }) + .collect(); + admitted.sort_by_key(|edge| (edge.source, edge.target)); + let mut adjacency: HashMap> = 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); - } + for edge in &admitted { + 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; + // Drop-stream LCG is independent of Leiden's per-replicate stream so + // graph-size RNG consumption cannot shift later replicates' samples. + let mut drop_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 { + for replicate in 0..n_replicates { // Perturb: drop each surviving positive edge independently with // the caller-supplied probability so cluster recovery is // stress-tested by resampling. Negative-effect edges are not // part of the partition graph at all, matching the adjacency. - let perturbed: Vec<&NetworkEdge> = edges + let perturbed: Vec<&NetworkEdge> = admitted .iter() - .filter(|edge| edge.effect > 0.0) + .copied() .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); + let draw_bits = crate::leiden::step_rng(&mut drop_state); + let draw = ((draw_bits >> 33) as f64 / (u64::MAX >> 33) as f64).min(1.0); draw >= edge_drop_probability }) .collect(); - let partition = greedy_union_partition(&perturbed, k_topics); + let mut leiden_state = + seed ^ 0xD1B5_4A32_D192_ED03 ^ (replicate as u64).wrapping_mul(0x9E37_79B9_7F4A_7C15); + let partition = leiden_partition(&perturbed, k_topics, &mut leiden_state); for i in 0..k_topics { for j in 0..k_topics { @@ -172,48 +186,6 @@ pub fn consensus_clusters( }) } -/// 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 { - let mut parent: Vec = (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 = 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 { @@ -278,8 +250,8 @@ mod tests { #[test] fn redundant_edges_between_clustered_topics_are_harmless() { - // A duplicate edge forces the union step onto an already-merged - // pair; output must stay identical to the simple-chain case. + // A duplicate edge must not change the Leiden partition of a + // strongly connected triple. let chain = vec![edge(0, 1, 0.95), edge(1, 2, 0.9)]; let duplicated = vec![edge(0, 1, 0.95), edge(1, 2, 0.9), edge(0, 2, 0.8)]; let plain = consensus_clusters(&chain, 3, 20, 0.99, 0.0, 31).unwrap(); @@ -367,4 +339,50 @@ mod tests { let second = consensus_clusters(&chain, 3, 15, 0.6, 0.1, 7).unwrap(); assert_eq!(first, second); } + + #[test] + fn two_cliques_with_a_weak_bridge_stay_two_consensus_clusters() { + // Operator-visible GAP-009 remainder: union-find glued both + // triangles; Leiden modularity plus co-assignment keeps them + // apart under a high consensus threshold and no drop noise. + let planted = vec![ + edge(0, 1, 1.0), + edge(1, 2, 1.0), + edge(0, 2, 1.0), + edge(3, 4, 1.0), + edge(4, 5, 1.0), + edge(3, 5, 1.0), + edge(2, 3, 0.01), + ]; + let output = consensus_clusters(&planted, 6, 20, 0.8, 0.0, 17).unwrap(); + assert_eq!(output.assignments[0], output.assignments[1]); + assert_eq!(output.assignments[1], output.assignments[2]); + assert_eq!(output.assignments[3], output.assignments[4]); + assert_eq!(output.assignments[4], output.assignments[5]); + assert_ne!(output.assignments[0], output.assignments[3]); + assert!(output.assignments[0].is_some()); + assert!(output.assignments[3].is_some()); + assert!(output.co_assignment[0][2] > 0.8); + assert!(output.co_assignment[0][3] < 0.8); + } + + #[test] + fn out_of_range_endpoints_do_not_create_singleton_clusters() { + // An edge whose target is outside k_topics is not a topic–topic + // association. Topic 0 must stay unclustered rather than gaining + // a singleton from that stray adjacency. + let stray = vec![edge(0, 2, 0.99), edge(1, 1, 0.8)]; + let output = consensus_clusters(&stray, 2, 12, 0.5, 0.0, 3).unwrap(); + assert_eq!(output.assignments[0], None); + assert_eq!(output.assignments[1], None); + } + + #[test] + fn permuting_equal_edges_does_not_change_consensus() { + let ordered = vec![edge(0, 1, 0.9), edge(1, 2, 0.8), edge(2, 3, 0.85)]; + let shuffled = vec![edge(2, 3, 0.85), edge(0, 1, 0.9), edge(1, 2, 0.8)]; + let first = consensus_clusters(&ordered, 4, 30, 0.6, 0.2, 99).unwrap(); + let second = consensus_clusters(&shuffled, 4, 30, 0.6, 0.2, 99).unwrap(); + assert_eq!(first, second); + } } diff --git a/crates/network_analysis/src/leiden.rs b/crates/network_analysis/src/leiden.rs new file mode 100644 index 000000000..975871905 --- /dev/null +++ b/crates/network_analysis/src/leiden.rs @@ -0,0 +1,711 @@ +//! Leiden community detection for admitted positive topic–topic edges. +//! +//! Replaces the greedy union-find stand-in. Each replicate partition is +//! Traag, Waltman, and van Eck (2019): fast local moving of nodes under +//! Newman–Girvan modularity (γ = 1), refinement that keeps communities +//! internally connected, and aggregation. Louvain is rejected because it +//! can emit internally disconnected communities. Union-find is rejected +//! because it merges every surviving edge into one component and makes +//! no modularity claim. +//! +//! Isolated topics remain singleton labels. The consensus layer still +//! maps topics that never appear in the perturbation graph to +//! unclustered. This module does not claim a causal cluster, a graphical +//! lasso, or an end-to-end export workflow. + +#![allow( + clippy::doc_markdown, + clippy::items_after_statements, + clippy::must_use_candidate, + clippy::cast_precision_loss, + clippy::cast_possible_truncation, + clippy::missing_errors_doc, + clippy::too_many_lines +)] + +use crate::edges::NetworkEdge; + +/// Deterministic LCG shared with the consensus perturbation stream. +pub(crate) fn step_rng(state: &mut u64) -> u64 { + *state = state + .wrapping_mul(6_364_136_223_846_793_005) + .wrapping_add(1_442_695_040_888_963_407); + *state +} + +fn shuffle(state: &mut u64, items: &mut [T]) { + for index in (1..items.len()).rev() { + let bound = (index + 1) as u64; + let other = (step_rng(state) % bound) as usize; + items.swap(index, other); + } +} + +struct Graph { + node_count: usize, + adjacency: Vec>, + strength: Vec, + total_weight: f64, +} + +impl Graph { + fn from_edges(edges: &[&NetworkEdge], node_count: usize) -> Self { + let mut adjacency = vec![Vec::new(); node_count]; + let mut strength = vec![0.0_f64; node_count]; + for edge in edges { + if edge.source >= node_count || edge.target >= node_count { + continue; + } + if edge.source == edge.target || edge.effect <= 0.0 { + continue; + } + adjacency[edge.source].push((edge.target, edge.effect)); + adjacency[edge.target].push((edge.source, edge.effect)); + strength[edge.source] += edge.effect; + strength[edge.target] += edge.effect; + } + for neighbors in &mut adjacency { + neighbors.sort_by(|left, right| left.0.cmp(&right.0).then(right.1.total_cmp(&left.1))); + neighbors.dedup_by(|later, kept| { + if later.0 == kept.0 { + kept.1 += later.1; + true + } else { + false + } + }); + } + let total_weight = strength.iter().sum::() / 2.0; + Self { + node_count, + adjacency, + strength, + total_weight, + } + } + + fn weight_to_community(&self, node: usize, membership: &[usize], community: usize) -> f64 { + self.adjacency[node] + .iter() + .filter(|(neighbor, _)| membership[*neighbor] == community) + .map(|(_, weight)| *weight) + .sum() + } +} + +/// Partition `k` topics with Leiden modularity on the surviving edges. +pub(crate) fn leiden_partition( + edges: &[&NetworkEdge], + topic_count: usize, + rng: &mut u64, +) -> Vec { + if topic_count == 0 { + return Vec::new(); + } + let mut graph = Graph::from_edges(edges, topic_count); + let mut membership: Vec = (0..topic_count).collect(); + if graph.total_weight <= 0.0 { + return membership; + } + + let mut leaf_members: Vec> = (0..topic_count).map(|node| vec![node]).collect(); + for _ in 0..topic_count { + let moved = move_nodes_fast(&graph, &mut membership, rng); + split_disconnected(&graph, &mut membership); + let refined = refine_partition(&graph, &membership, rng); + let (aggregate, next_membership, next_leaves) = + aggregate_graph(&graph, &membership, &refined, &leaf_members); + if aggregation_level_is_stable(aggregate.node_count, graph.node_count, moved) { + break; + } + if aggregate.node_count <= 1 { + membership = next_membership; + leaf_members = next_leaves; + break; + } + graph = aggregate; + membership = next_membership; + leaf_members = next_leaves; + } + + project_original_labels(&leaf_members, &membership, topic_count) +} + +/// Stop only when local moving made no change and refinement did not reduce +/// the graph. A move that leaves the aggregate the same size continues so +/// the next Leiden pass starts from the updated partition (Traag, Waltman, +/// & van Eck, 2019). +fn aggregation_level_is_stable(aggregate_nodes: usize, graph_nodes: usize, moved: bool) -> bool { + aggregate_nodes == graph_nodes && !moved +} + +fn project_original_labels( + leaf_members: &[Vec], + membership: &[usize], + topic_count: usize, +) -> Vec { + let mut labels = vec![0_usize; topic_count]; + for (supernode, members) in leaf_members.iter().enumerate() { + let community = membership[supernode]; + for &node in members { + labels[node] = community; + } + } + dense_labels(&labels) +} + +fn dense_labels(membership: &[usize]) -> Vec { + let mut remap = vec![ + None; + membership + .iter() + .copied() + .max() + .map_or(0, |value| value + 1) + ]; + let mut next = 0_usize; + membership + .iter() + .map(|&community| { + *remap[community].get_or_insert_with(|| { + let label = next; + next += 1; + label + }) + }) + .collect() +} + +fn community_totals(graph: &Graph, membership: &[usize]) -> Vec { + let community_count = membership + .iter() + .copied() + .max() + .map_or(0, |value| value + 1); + let mut totals = vec![0.0_f64; community_count]; + for (node, &community) in membership.iter().enumerate() { + totals[community] += graph.strength[node]; + } + totals +} + +fn modularity_delta( + strength: f64, + weight_from: f64, + weight_to: f64, + total_from: f64, + total_to: f64, + total_weight: f64, +) -> f64 { + // Moving node i from community c to d, γ = 1: + // ΔQ = (k_{i→d} − k_{i→c}) / m − k_i (Σ_d − Σ_c + k_i) / (2 m²) + // `total_from` still includes `strength`. + let two_m = 2.0 * total_weight; + (weight_to - weight_from) / total_weight + - strength * (total_to - total_from + strength) / (two_m * total_weight) +} + +fn move_nodes_fast(graph: &Graph, membership: &mut [usize], rng: &mut u64) -> bool { + if graph.node_count <= 1 { + return false; + } + if graph.total_weight <= 0.0 { + return false; + } + let mut totals = community_totals(graph, membership); + let mut queue: Vec = (0..graph.node_count).collect(); + shuffle(rng, &mut queue); + let mut queued = vec![true; graph.node_count]; + let mut any_move = false; + let mut cursor = 0_usize; + while cursor < queue.len() { + let node = queue[cursor]; + cursor += 1; + queued[node] = false; + let from = membership[node]; + let weight_from = graph.weight_to_community(node, membership, from); + let mut best_community = from; + let mut best_delta = 0.0_f64; + let mut seen = Vec::new(); + for &(neighbor, _) in &graph.adjacency[node] { + let candidate = membership[neighbor]; + if candidate == from || seen.contains(&candidate) { + continue; + } + seen.push(candidate); + let weight_to = graph.weight_to_community(node, membership, candidate); + let delta = modularity_delta( + graph.strength[node], + weight_from, + weight_to, + totals[from], + totals[candidate], + graph.total_weight, + ); + if delta > best_delta { + best_delta = delta; + best_community = candidate; + } + } + if best_community == from { + continue; + } + membership[node] = best_community; + totals[from] -= graph.strength[node]; + totals[best_community] += graph.strength[node]; + any_move = true; + for &(neighbor, _) in &graph.adjacency[node] { + if membership[neighbor] != best_community && !queued[neighbor] { + queued[neighbor] = true; + queue.push(neighbor); + } + } + } + any_move +} + +fn split_disconnected(graph: &Graph, membership: &mut [usize]) { + let community_count = membership + .iter() + .copied() + .max() + .map_or(0, |value| value + 1); + let mut next_label = community_count; + for community in 0..community_count { + let members: Vec = membership + .iter() + .enumerate() + .filter_map(|(node, &label)| (label == community).then_some(node)) + .collect(); + if members.len() <= 1 { + continue; + } + let mut index_of = vec![0_usize; graph.node_count]; + for (index, &node) in members.iter().enumerate() { + index_of[node] = index; + } + let mut visited = vec![false; members.len()]; + let mut first_component = true; + for start_index in 0..members.len() { + if visited[start_index] { + continue; + } + let mut stack = vec![start_index]; + visited[start_index] = true; + let mut component = Vec::new(); + while let Some(index) = stack.pop() { + component.push(members[index]); + for &(neighbor, _) in &graph.adjacency[members[index]] { + if membership[neighbor] != community { + continue; + } + let neighbor_index = index_of[neighbor]; + if !visited[neighbor_index] { + visited[neighbor_index] = true; + stack.push(neighbor_index); + } + } + } + if first_component { + first_component = false; + continue; + } + for node in component { + membership[node] = next_label; + } + next_label += 1; + } + } +} + +fn refine_partition(graph: &Graph, membership: &[usize], rng: &mut u64) -> Vec { + let mut refined: Vec = (0..graph.node_count).collect(); + let community_count = membership + .iter() + .copied() + .max() + .map_or(0, |value| value + 1); + for community in 0..community_count { + let members: Vec = membership + .iter() + .enumerate() + .filter_map(|(node, &label)| (label == community).then_some(node)) + .collect(); + if members.len() <= 1 { + continue; + } + let mut sub_membership: Vec = (0..members.len()).collect(); + let mut sub_adj = vec![Vec::new(); members.len()]; + let mut sub_strength = vec![0.0_f64; members.len()]; + let mut index_of = vec![usize::MAX; graph.node_count]; + for (index, &node) in members.iter().enumerate() { + index_of[node] = index; + } + for (index, &node) in members.iter().enumerate() { + for &(neighbor, weight) in &graph.adjacency[node] { + if membership[neighbor] != community { + continue; + } + let neighbor_index = index_of[neighbor]; + sub_adj[index].push((neighbor_index, weight)); + sub_strength[index] += weight; + } + } + let sub_total = sub_strength.iter().sum::() / 2.0; + if sub_total <= 0.0 { + continue; + } + let subgraph = Graph { + node_count: members.len(), + adjacency: sub_adj, + strength: sub_strength, + total_weight: sub_total, + }; + move_nodes_fast(&subgraph, &mut sub_membership, rng); + split_disconnected(&subgraph, &mut sub_membership); + let offset = refined.iter().copied().max().map_or(0, |value| value + 1); + for (index, &node) in members.iter().enumerate() { + refined[node] = offset + sub_membership[index]; + } + } + dense_labels(&refined) +} + +fn aggregate_graph( + graph: &Graph, + parent: &[usize], + refined: &[usize], + leaf_members: &[Vec], +) -> (Graph, Vec, Vec>) { + let refined_count = refined.iter().copied().max().map_or(0, |value| value + 1); + let mut next_leaves = vec![Vec::new(); refined_count]; + for node in 0..graph.node_count { + next_leaves[refined[node]].extend_from_slice(&leaf_members[node]); + } + let mut adjacency = vec![Vec::new(); refined_count]; + let mut strength = vec![0.0_f64; refined_count]; + for node in 0..graph.node_count { + strength[refined[node]] += graph.strength[node]; + let source = refined[node]; + for &(neighbor, weight) in &graph.adjacency[node] { + let target = refined[neighbor]; + if source < target { + adjacency[source].push((target, weight)); + adjacency[target].push((source, weight)); + } + } + } + for neighbors in &mut adjacency { + neighbors.sort_by_key(|left| left.0); + neighbors.dedup_by(|later, kept| { + if later.0 == kept.0 { + kept.1 += later.1; + true + } else { + false + } + }); + } + let aggregate = Graph { + node_count: refined_count, + adjacency, + strength, + total_weight: graph.total_weight, + }; + let mut parent_of_refined = vec![0_usize; refined_count]; + for node in 0..graph.node_count { + parent_of_refined[refined[node]] = parent[node]; + } + let next_membership = dense_labels(&parent_of_refined); + (aggregate, next_membership, next_leaves) +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::edges::NetworkEdge; + + 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, + } + } + + fn refs(owned: &[NetworkEdge]) -> Vec<&NetworkEdge> { + owned.iter().collect() + } + + fn communities_connected(edges: &[&NetworkEdge], labels: &[usize]) -> bool { + let graph = Graph::from_edges(edges, labels.len()); + let community_count = labels.iter().copied().max().map_or(0, |value| value + 1); + for community in 0..community_count { + let members: Vec = labels + .iter() + .enumerate() + .filter_map(|(node, &label)| (label == community).then_some(node)) + .collect(); + if members.len() <= 1 { + continue; + } + let mut index_of = vec![0_usize; labels.len()]; + for (index, &node) in members.iter().enumerate() { + index_of[node] = index; + } + let mut seen = vec![false; members.len()]; + let mut stack = vec![0_usize]; + seen[0] = true; + let mut visited = 1_usize; + while let Some(index) = stack.pop() { + for &(neighbor, _) in &graph.adjacency[members[index]] { + if labels[neighbor] != community { + continue; + } + let neighbor_index = index_of[neighbor]; + if !seen[neighbor_index] { + seen[neighbor_index] = true; + visited += 1; + stack.push(neighbor_index); + } + } + } + if visited != members.len() { + return false; + } + } + true + } + + #[test] + fn empty_and_edgeless_graphs_stay_singletons() { + let mut rng = 1_u64; + assert!(leiden_partition(&[], 0, &mut rng).is_empty()); + let labels = leiden_partition(&[], 4, &mut rng); + assert_eq!(labels, vec![0, 1, 2, 3]); + assert!(dense_labels(&[]).is_empty()); + let totals = community_totals(&Graph::from_edges(&[], 0), &[]); + assert!(totals.is_empty()); + } + + #[test] + fn out_of_range_self_and_nonpositive_edges_are_ignored() { + let owned = vec![ + edge(0, 9, 1.0), + edge(9, 0, 1.0), + edge(1, 1, 1.0), + edge(0, 1, 0.0), + edge(0, 1, -0.4), + ]; + let edges = refs(&owned); + let mut rng = 2_u64; + let labels = leiden_partition(&edges, 2, &mut rng); + assert_eq!(labels, vec![0, 1]); + } + + #[test] + fn two_nodes_with_a_positive_edge_share_a_community() { + let owned = vec![edge(0, 1, 0.9)]; + let edges = refs(&owned); + let mut rng = 3_u64; + let labels = leiden_partition(&edges, 2, &mut rng); + assert_eq!(labels[0], labels[1]); + assert!(communities_connected(&edges, &labels)); + } + + #[test] + fn two_cliques_joined_by_a_weak_bridge_stay_apart() { + // Union-find would glue both triangles into one component. Leiden + // modularity recovers the planted pair of communities. + let owned = vec![ + edge(0, 1, 1.0), + edge(1, 2, 1.0), + edge(0, 2, 1.0), + edge(3, 4, 1.0), + edge(4, 5, 1.0), + edge(3, 5, 1.0), + edge(2, 3, 0.01), + ]; + let edges = refs(&owned); + let mut rng = 11_u64; + let labels = leiden_partition(&edges, 6, &mut rng); + assert_eq!(labels[0], labels[1]); + assert_eq!(labels[1], labels[2]); + assert_eq!(labels[3], labels[4]); + assert_eq!(labels[4], labels[5]); + assert_ne!(labels[0], labels[3]); + assert!(communities_connected(&edges, &labels)); + } + + #[test] + fn duplicate_edges_are_summed_and_remain_connected() { + let owned = vec![edge(0, 1, 0.4), edge(0, 1, 0.4), edge(1, 2, 0.8)]; + let edges = refs(&owned); + let mut rng = 7_u64; + let labels = leiden_partition(&edges, 3, &mut rng); + assert_eq!(labels[0], labels[1]); + assert_eq!(labels[1], labels[2]); + assert!(communities_connected(&edges, &labels)); + } + + #[test] + fn identical_seed_reproduces_the_partition() { + let owned = vec![ + edge(0, 1, 0.9), + edge(1, 2, 0.85), + edge(3, 4, 0.88), + edge(2, 3, 0.2), + ]; + let edges = refs(&owned); + let mut first_rng = 21_u64; + let mut second_rng = 21_u64; + let first = leiden_partition(&edges, 5, &mut first_rng); + let second = leiden_partition(&edges, 5, &mut second_rng); + assert_eq!(first, second); + assert_eq!(first.len(), 5); + assert_ne!(first[4], first[0]); + assert!(communities_connected(&edges, &first)); + } + + #[test] + fn isolated_topic_keeps_a_distinct_label() { + let owned = vec![edge(0, 1, 0.95)]; + let edges = refs(&owned); + let mut rng = 5_u64; + let labels = leiden_partition(&edges, 3, &mut rng); + assert_eq!(labels[0], labels[1]); + assert_ne!(labels[2], labels[0]); + assert!(communities_connected(&edges, &labels)); + } + + #[test] + fn disconnected_membership_is_split_into_connected_components() { + let owned = vec![edge(0, 1, 1.0), edge(2, 3, 1.0)]; + let edges = refs(&owned); + let graph = Graph::from_edges(&edges, 4); + let mut membership = vec![0, 0, 0, 0]; + split_disconnected(&graph, &mut membership); + assert_eq!(membership[0], membership[1]); + assert_eq!(membership[2], membership[3]); + assert_ne!(membership[0], membership[2]); + } + + #[test] + fn single_node_graph_does_not_move() { + let graph = Graph { + node_count: 1, + adjacency: vec![Vec::new()], + strength: vec![1.0], + total_weight: 0.5, + }; + let mut membership = vec![0]; + let mut rng = 9_u64; + assert!(!move_nodes_fast(&graph, &mut membership, &mut rng)); + assert_eq!(membership, vec![0]); + } + + #[test] + fn edgeless_multi_node_graph_does_not_move() { + let graph = Graph::from_edges(&[], 3); + let mut membership = vec![0, 1, 2]; + let mut rng = 4_u64; + assert!(!move_nodes_fast(&graph, &mut membership, &mut rng)); + assert_eq!(membership, vec![0, 1, 2]); + } + + #[test] + fn refine_skips_edgeless_communities_and_merges_a_triangle() { + let owned = vec![edge(0, 1, 1.0), edge(1, 2, 1.0), edge(0, 2, 1.0)]; + let edges = refs(&owned); + let graph = Graph::from_edges(&edges, 4); + let membership = vec![0, 0, 0, 1]; + let mut rng = 8_u64; + let refined = refine_partition(&graph, &membership, &mut rng); + assert_eq!(refined[0], refined[1]); + assert_eq!(refined[1], refined[2]); + assert_ne!(refined[3], refined[0]); + + let empty_pair = Graph::from_edges(&[], 2); + let pair_membership = vec![0, 0]; + let skipped = refine_partition(&empty_pair, &pair_membership, &mut rng); + assert_eq!(skipped.len(), 2); + } + + #[test] + fn modularity_delta_is_negative_when_leaving_a_stronger_community() { + let delta = modularity_delta(1.0, 2.0, 0.01, 3.0, 2.0, 4.0); + assert!(delta < 0.0); + } + + #[test] + fn shuffle_of_a_single_item_is_a_no_op() { + let mut rng = 13_u64; + let mut items = vec![7_usize]; + shuffle(&mut rng, &mut items); + assert_eq!(items, vec![7]); + let mut empty: Vec = Vec::new(); + shuffle(&mut rng, &mut empty); + assert!(empty.is_empty()); + } + + #[test] + fn aggregate_sums_parallel_supernode_edges_and_keeps_distinct_neighbors() { + // Two crossing edges between supernodes 0 and 1 plus one edge + // from 1 to 2. Aggregation must merge the parallel pair and keep + // the distinct neighbour. + let owned = vec![ + edge(0, 1, 1.0), + edge(2, 3, 1.0), + edge(4, 5, 1.0), + edge(0, 2, 0.3), + edge(1, 3, 0.4), + edge(2, 4, 0.2), + ]; + let edges = refs(&owned); + let graph = Graph::from_edges(&edges, 6); + let refined = vec![0, 0, 1, 1, 2, 2]; + let parent = refined.clone(); + let leaf_members: Vec> = (0..6).map(|node| vec![node]).collect(); + let (aggregate, next_membership, next_leaves) = + aggregate_graph(&graph, &parent, &refined, &leaf_members); + assert_eq!(aggregate.node_count, 3); + assert_eq!(aggregate.adjacency[0].len(), 1); + assert_eq!(aggregate.adjacency[0][0].0, 1); + assert!((aggregate.adjacency[0][0].1 - 0.7).abs() < 1e-15); + assert_eq!(aggregate.adjacency[1].len(), 2); + assert_eq!(aggregate.adjacency[1][0].0, 0); + assert_eq!(aggregate.adjacency[1][1].0, 2); + assert!((aggregate.adjacency[1][0].1 - 0.7).abs() < 1e-15); + assert!((aggregate.adjacency[1][1].1 - 0.2).abs() < 1e-15); + assert_eq!(next_membership, vec![0, 1, 2]); + assert_eq!(next_leaves[0], vec![0, 1]); + assert_eq!(next_leaves[1], vec![2, 3]); + assert_eq!(next_leaves[2], vec![4, 5]); + } + + #[test] + fn connectedness_oracle_rejects_a_disconnected_community() { + let owned = vec![edge(0, 1, 1.0), edge(2, 3, 1.0)]; + let edges = refs(&owned); + assert!(!communities_connected(&edges, &[0, 0, 0, 0])); + assert!(communities_connected(&edges, &[0, 0, 1, 1])); + assert!(communities_connected(&edges, &[0, 1, 2, 3])); + } + + #[test] + fn aggregation_continues_when_nodes_moved_on_a_same_size_level() { + // Unique-branch coverage of the `!moved` false arm: refinement + // did not shrink the graph, but local moving still changed the + // partition, so Leiden must not stop (Traag et al., 2019). + assert!(aggregation_level_is_stable(4, 4, false)); + assert!(!aggregation_level_is_stable(4, 4, true)); + assert!(!aggregation_level_is_stable(2, 4, false)); + assert!(!aggregation_level_is_stable(1, 4, true)); + assert!(aggregation_level_is_stable(1, 1, false)); + assert!(!aggregation_level_is_stable(1, 1, true)); + } +} diff --git a/crates/network_analysis/src/lib.rs b/crates/network_analysis/src/lib.rs index dfa99347b..886f44c91 100644 --- a/crates/network_analysis/src/lib.rs +++ b/crates/network_analysis/src/lib.rs @@ -14,6 +14,7 @@ mod cluster; mod consensus; mod edges; mod error; +mod leiden; mod simplex; mod stability; diff --git a/docs/TRACEABILITY.md b/docs/TRACEABILITY.md index 2b783c2ab..489bf54d1 100644 --- a/docs/TRACEABILITY.md +++ b/docs/TRACEABILITY.md @@ -75,7 +75,7 @@ The full APA 7th standards/literature register remains `docs/research/standards- | report template/section/copied/style/modality method effects | ADR 0004/0012; PRD/TRD | simulation truth factors implemented; `corpus_background` background-versus-unique-content identity on the active PR; estimator-side method model remains future | partial | | report template/section/copied/style/modality method effects | ADR 0004/0012; PRD/TRD | simulation truth factors implemented; `prompt_source` prompt-versus-unique-content identity on the active PR; estimator-side method model remains future | partial | | candidate K statistical/Pareto gates | ADR 0012; research | `model_selection` fits each candidate `K` with the CPU `f64` reference and scores the actual mixture likelihood plus Schwarz's (1978) `ℓ − (p ln N)/2` penalty before the Pareto gate; candidate blinding, blinded LLM review, GPU, and backend comparison remain accepted-target | active-PR | -| compositional topic correlation / stable clustering | ADR 0005/0012; research | future `network_analysis` | accepted-target | +| compositional topic correlation / stable clustering | ADR 0005/0012; research | `network_analysis` ILR posterior edges, BH admission, bootstrap stability, and Traag et al. (2019) Leiden consensus on this PR (`docs/research/leiden-consensus.md`); graphical lasso and export workflow remain accepted-target | partial | | posterior ESEM / longitudinal invariance / DSEM | ADR 0005 | `psychometric_core` construct/input gates, true-loading OLS recovery, posterior-draw point-estimate averaging, Rubin `T` on draw-level OLS loadings, CWC within/between OLS plus the contextual effect, event-time log-rate, constant- and time-varying-predictor discrete effects (Voelkle Eqs. 12 and 14), exact scalar discrete process noise (Driver et al., 2017, Eq. 3), lagged latent covariance and unconditional latent variance (Driver et al., 2017, Eq. 3–4), stationary within-subject variance (Driver et al., 2017, Eq. 4 as `Δt → ∞`; `asymDIFFUSION`), trait-plus-state variance (Driver et al., 2017, §4.3 `TRAITVAR`; not process noise), observed-indicator variance and lagged observed covariance (Driver et al., 2017, Eq. 5; Table 2 `MANIFESTVAR` is `Θ`, not `Var(y)`; `MANIFESTTRAITVAR` is not `MANIFESTVAR`; `Θ` does not enter lagged observed covariance; observed-indicator mean is `τ + λ μ`; `MANIFESTMEANS` is not `E(y)`; `CINT` is not `MANIFESTMEANS`; discrete latent mean is `exp(a Δt) μ_0 + (exp(a Δt) − 1)/a κ`; `T0MEANS` is not `μ_t`; evolved observed mean is `τ + λ μ_t`; `τ + λ μ_0` is not `E(y_t)`; contemporaneous `TDPREDEFFECT` impulse is `m x`, not `CINT`, not `TIPREDEFFECT`, and not Voelkle Eq. 14; Eq. 5 of that contemporaneous impulse is `τ + λ(μ_t + m x)`, and `τ + λ μ_t` is not that observed mean; time-independent `TIPREDEFFECT` increment is `A^{-1}[e^{A Δt} − I] B z`, not `CINT`, not `M x`, not Voelkle Eq. 14, and not the coefficient `B`; Eq. 5 of that increment is `τ + λ(μ_t + A^{-1}[e^{A Δt} − I] B z)`, and `τ + λ μ_t` is not that observed mean; `τ + λ(μ_t + m x)` is not that observed mean; `τ + λ(μ_t + e^{a(t−u)} m x)` is not that observed mean when `u ≠ t`; within-interval `TDPREDEFFECT` carry is `e^{A(t−u)} M x` for `t0 < u < t`, not the contemporaneous Dirac, not `CINT`, not `TIPREDEFFECT`, and not Voelkle Eq. 14; Eq. 5 of that carry is `τ + λ(μ_t + e^{a(t−u)} m x)`, and `τ + λ μ_t` is not that observed mean; `τ + λ(μ_t + m x)` is not that carried observed mean when `u ≠ t`; §7.2 level-change `CINT` is `κ = −a m x` (`a < 0`; not the dissipating Dirac, not a free `CINT`, not `TIPREDEFFECT`; Eq. 3 of that setting is `(1 − e^{a Δt}) m x`); §7.2 extra-process contribution is `a_{ηξ} x (e^{ε Δt} − e^{a Δt}) / (ε − a)` (not `κ = −a m x`, not `(1 − e^{a Δt}) m x`, not the dissipating Dirac; `ε ≥ 0` fails closed; Eq. 5 of that contribution is `τ + λ(μ_t + a_{ηξ} x (e^{ε Δt} − e^{a Δt}) / (ε − a)`; extra `LAMBDA` is 0; `τ + λ μ_t` is not that observed mean; after-t0 extra-process `TDPREDEFFECT` uses `t − u` with `t0 < u < t` while `μ_t` uses `Δt`; that after-t0 observed mean is not the first-occasion extra-process observed mean; §7.2 `asymTIPREDEFFECT` is `-B z / a` for `a < 0` and is not `B`, not `A^{-1}[e^{A Δt} − I] B z`, not `CINT`, and not `M x`; §7.2 `addedTIPREDVAR` is `(B / a)² v` and is not `TRAITVAR`, not `asymDIFFUSION`, and not `-B z / a`; Table 2 `asymCINT` is `-κ / a` for `a < 0` and is not `κ`, not `A^{-1}[e^{A Δt} − I] κ`, not `T0MEANS`, and not `-B z / a`; p. 16 stationary `T0MEANS` is `-κ / a + −B z / a` and is not free `T0MEANS`, not `asymCINT` alone, not `asymTIPREDEFFECT` alone, and not the finite-interval discrete latent mean; Eq. 5 of that constrained mean is `τ + λ(−κ / a + −B z / a)`; `τ + λ μ_0` is not that observed mean; `MANIFESTMEANS` is not `E(y_0)`; the constrained latent mean is not `E(y_0)`; stationary `T0VAR` is `trait + −q / (2 a) + (B / a)² v` (not free `T0VAR`, not `asymDIFFUSION` alone, not `TRAITVAR` alone, not `addedTIPREDVAR` alone, and not the finite-interval discrete latent variance. Eq. 5 of that constrained variance is `λ²(trait + −q / (2 a) + (B / a)² v) + θ + ψ` (JSS PDF re-opened 2026-08-22T03:20Z; form the stationary latent variance first, then `λ² p + θ + ψ`; `λ² p_0` is not that observed variance; `λ²(−q / (2 a)) + θ` is not that observed variance when `TRAITVAR` or `addedTIPREDVAR` is nonzero; `MANIFESTVAR` is not `Var(y_0)`; the constrained latent variance is not `Var(y_0)`)); lagged stationary `T0VAR` is `trait + e^{a Δt}(−q / (2 a)) + (B / a)² v` (trait and `addedTIPREDVAR` do not decay; contemporaneous `T0VAR` is not that lagged map; decaying the constrained total as if it were all state is not that lagged map; Eq. 5 of that lagged covariance is `λ²(trait + e^{a Δt}(−q / (2 a)) + (B / a)² v) + ψ`; `Θ` does not enter; contemporaneous `Var(y_0)` is not that lagged observed covariance; the lagged latent covariance is not that observed covariance); later-occasion stationary `T0VAR` is `trait + e^{2 a Δt}(−q / (2 a)) + Q_Δt + (B / a)² v` (trait and `addedTIPREDVAR` do not enter `Q_Δt`; under stationarity that composition equals contemporaneous `T0VAR`; evolving the constrained total as if it were all state is not that later map; the lagged covariance omits `Q_Δt`; `Q_Δt` is not that later map; Eq. 5 of that later-occasion variance is `λ²(trait + e^{2 a Δt}(−q / (2 a)) + Q_Δt + (B / a)² v) + θ + ψ`; lagged observed covariance omits `Q_Δt` and `θ`; `MANIFESTVAR` is not `Var(y_t)`; the later-occasion latent variance is not `Var(y_t)`); predetermined later-occasion `T0VAR` is `trait + e^{2 a Δt} p_0 + Q_Δt + (B / a)² v` (free `T0VAR` `p_0` is not that later map; setting `p_0 = −q / (2 a)` recovers the stationary later-occasion map; stationary later variance uses `−q / (2 a)` in place of `p_0` and is not that later map when `p_0` is free; evolving `trait + p_0 + (B / a)² v` as if it were all state is not that later map; Eq. 5 of that predetermined later-occasion variance is `λ²(trait + e^{2 a Δt} p_0 + Q_Δt + (B / a)² v) + θ + ψ`; `MANIFESTVAR` is not `Var(y_t)`; the predetermined later-occasion latent variance is not `Var(y_t)`; stationary later observed variance is not that observed variance when `p_0` is free); predetermined lagged `T0VAR` is `trait + e^{a Δt} p_0 + (B / a)² v` (free `T0VAR` `p_0` is not that lagged map; setting `p_0 = −q / (2 a)` recovers the stationary lagged map; stationary lagged covariance uses `−q / (2 a)` in place of `p_0` and is not that lagged map when `p_0` is free; evolving `trait + p_0 + (B / a)² v` as if it were all state is not that lagged map; later-occasion variance includes `Q_Δt` and is not that lagged map; Eq. 5 of that predetermined lagged covariance is `λ²(trait + e^{a Δt} p_0 + (B / a)² v) + ψ`; `MANIFESTVAR` does not enter; the predetermined lagged latent covariance is not that observed covariance; predetermined later observed variance includes `Q_Δt` and `θ` and is not that lagged observed covariance; stationary lagged observed covariance is not that observed covariance when `p_0` is free; the predetermined first-occasion variance of §4.3 predetermined `T0VAR` is `trait + p_0 + (B / a)² v`; free `p_0` is not that map; stationary first-occasion variance uses `−q / (2 a)` in place of `p_0` and is not that map when `p_0` is free; lagged covariance decays the state and is not that map; later-occasion variance includes `Q_Δt` and is not that map; Eq. 5 of that predetermined first-occasion variance is `λ²(trait + p_0 + (B / a)² v) + θ + ψ`; `MANIFESTVAR` is not that first-occasion observed variance; the predetermined first-occasion latent variance is not that observed variance; stationary first-occasion observed variance is not that observed variance when `p_0` is free; predetermined later observed variance includes `Q_Δt` and is not that first-occasion observed variance; later-start lagged covariance of predetermined `T0VAR` is `trait + e^{a s}(e^{2 a u} p_0 + Q_u) + (B / a)² v` (Driver et al., 2017, §4.3 `startoffset`; Eq. 4; JSS PDF re-opened 2026-08-23T10:27Z; first-occasion lagged omits `e^{a s} Q_u`; later-occasion variance does not lag; stationary lagged uses `−q / (2 a)`; decaying the later total is not that map; Eq. 5 of that later-start lagged covariance is `λ²` of it plus `ψ`; `Θ` does not enter; first-occasion lagged observed omits `e^{a s} Q_u`; later observed variance includes `Q_u` and `θ`; later-start later-occasion variance of predetermined `T0VAR` is `trait + e^{2 a s}(e^{2 a u} p_0 + Q_u) + Q_s + (B / a)² v` (Driver et al., 2017, §4.3 `startoffset`; Eq. 3–4 Chapman–Kolmogorov `Q_{u+s} = e^{2 a s} Q_u + Q_s`; JSS PDF re-opened 2026-08-23T11:05Z; later-occasion variance at `u` omits `Q_s`; later-start lagged covariance omits `Q_s`; stationary later uses `−q / (2 a)`; evolving the later total as if it were all state is not that map; ignoring `startoffset` omits `e^{2 a s} Q_u`; Eq. 5 of that later-start later-occasion variance is `λ²` of it plus `θ + ψ`; `MANIFESTVAR` is not that observed variance; p. 16 `discreteDRIFTstd` is `e^{a Δt}` after strictly positive `asymDIFFUSION` `-q / (2 a)` (footnote 4; unstandardised `e^{a Δt}` is defined for growing `a ≥ 0` and for zero diffusion and is not `discreteDRIFTstd`; the §7.1 trait-plus-state autocorrelation uses `TRAITVAR` and is not `discreteDRIFTstd`; p. 16 `discreteDIFFUSIONstd` is `Q_Δt / (−q / (2 a))` after strictly positive `asymDIFFUSION` `-q / (2 a)` (footnote 4; unstandardised `Q_Δt` is defined for growing `a ≥ 0` and for zero diffusion and is not `discreteDIFFUSIONstd`; the continuous standardisation `−2 a` is not `discreteDIFFUSIONstd`; `Q_Δt / (trait + p + added)` uses `TRAITVAR` and is not `discreteDIFFUSIONstd`; `TRAITVAR` is not the standardisation variance; p. 16 `DIFFUSIONstd` is `q / (−q / (2 a)) = −2 a` after strictly positive `asymDIFFUSION` `-q / (2 a)` (Driver et al., 2017, p. 16; Eq. 4; footnote 4; JSS PDF re-opened 2026-08-23T13:20Z; unstandardised `q` is defined for growing `a ≥ 0` and for zero diffusion and is not `DIFFUSIONstd`; the discrete standardisation `Q_Δt / (−q / (2 a))` depends on `Δt` and is not `DIFFUSIONstd`; `q / (trait + p + added)` uses `TRAITVAR` and is not `DIFFUSIONstd`; `TRAITVAR` is not the standardisation variance; p. 16 `DRIFTstd` is the continuous auto-effect after strictly positive `asymDIFFUSION` `-q / (2 a)` (Driver et al., 2017, p. 16; Eq. 1; footnote 4; JSS PDF re-opened 2026-08-23T13:28Z); unstandardised `a` is defined for growing `a ≥ 0` and for zero diffusion and is not `DRIFTstd`; the discrete standardisation `e^{a Δt}` depends on the event interval and is not `DRIFTstd`; `a p / (trait + p + added)` uses `TRAITVAR` and is not `DRIFTstd`; `TRAITVAR` is not the standardisation variance); p. 16 `asymTIPREDEFFECTstd` is `(-B / a) · √v / √(-q / (2 a))` after strictly positive `asymDIFFUSION` `-q / (2 a)` and strictly positive predictor variance `v` (Driver et al., 2017, p. 16; §7.2; footnote 4; JSS PDF re-opened 2026-08-23T14:25Z; unstandardised `-B / a` is defined for a zero coefficient and for zero predictor variance and is not `asymTIPREDEFFECTstd`; the finite-interval standardisation `A^{-1}[e^{A Δt} − I] B · √v / √p` depends on the event interval and is not `asymTIPREDEFFECTstd`; `(-B / a) · √v / √(trait + p + added)` uses `TRAITVAR` and is not `asymTIPREDEFFECTstd`; `TRAITVAR` is not the standardisation variance); p. 16 `TIPREDEFFECTstd` is `B · √v / √(-q / (2 a))` after strictly positive `asymDIFFUSION` `-q / (2 a)` and strictly positive predictor variance `v` (Driver et al., 2017, p. 16; §7.2; footnote 4; JSS PDF re-opened 2026-08-23T16:21Z; unstandardised `B` is defined for a zero coefficient and for zero predictor variance and is not `TIPREDEFFECTstd`; the asymptotic standardisation `(-B / a) · √v / √p` is the total change and is not `TIPREDEFFECTstd`; the finite-interval standardisation `A^{-1}[e^{A Δt} − I] B · √v / √p` depends on the event interval and is not `TIPREDEFFECTstd`; `B · √v / √(trait + p + added)` uses `TRAITVAR` and is not `TIPREDEFFECTstd`; `TRAITVAR` is not the standardisation variance); Table 3 `T0TIPREDEFFECTstd` is `t0_b · √v / √p_0` after strictly positive free `T0VAR` `p_0` and strictly positive predictor variance `v` (Driver et al., 2017, Table 3, p. 13; p. 16; footnote 4; 2017-era ctsem `summary.ctsemFit.R`; JSS PDF re-opened 2026-08-23T17:20Z; the affected variance is free `T0VAR`, not `asymDIFFUSION`; unstandardised `t0_b` is defined for a zero coefficient and for zero predictor variance and is not `T0TIPREDEFFECTstd`; `TIPREDEFFECTstd` `B · √v / √(-q / (2 a))` is the continuous coefficient and is not `T0TIPREDEFFECTstd`; `asymTIPREDEFFECTstd` `(-B / a) · √v / √p` is the total change and is not `T0TIPREDEFFECTstd`; `t0_b · √v / √(trait + p_0 + added)` uses `TRAITVAR` and is not `T0TIPREDEFFECTstd`; `TRAITVAR` is not the standardisation variance); 2017-era `addedT0TIPREDVAR` is `t0_b² v` (Driver et al., 2017, Table 3, p. 13; p. 16; §7.2; 2017-era ctsem `summary.ctsemFit.R`; JSS PDF re-opened 2026-08-23T18:20Z; `T0TIPREDEFFECT %*% TIPREDVAR %*% t(T0TIPREDEFFECT)` immediately after `T0TIPREDEFFECTstd`; form `t0_b` first, then square, then multiply by `v`; a zero coefficient or zero predictor variance is exactly zero; free `T0TIPREDEFFECT` does not require `a < 0`; `(B / a)² v` is `addedTIPREDVAR` and is not this first-occasion map; `t0_b · √v / √p_0` is `T0TIPREDEFFECTstd` and is not this variance; free `T0VAR` is not this extra TI variance; `TRAITVAR` is not this extra TI variance; Equation 5 of 2017-era `addedT0TIPREDVAR` is `λ² t0_b² v` (Driver et al., 2017, Eq. 5, p. 5; Table 3, p. 13; Table 2, p. 12; 2017-era ctsem `summary.ctsemFit.R`; JSS PDF re-opened 2026-08-23T19:10Z; form `t0_b² v` first, then `(λ extra) λ` with `θ = 0`; a zero loading or zero extra is exactly zero; `t0_b² v` is the latent extra, not the observed extra; `λ² p_0 + θ` is first-occasion observed variance, not this extra; `λ² (B / a)² v` is Eq. 5 of `addedTIPREDVAR`, not this first-occasion observed extra; `MANIFESTVAR` `θ` is not this extra; Equation 5 of §7.2 `addedTIPREDVAR` is `λ² (B / a)² v`; form `(B / a)² v` first, then `(λ extra) λ` with `θ = 0`; a zero loading or zero extra is exactly zero; lasting asymptotic extra requires `a < 0`; `(B / a)² v` is the latent extra, not the observed extra; `λ² t0_b² v` is first-occasion extra observed TI variance, not this extra; `λ² p + θ` is stationary observed variance, not this extra; `MANIFESTVAR` `θ` is not this extra; p. 16 `TDPREDEFFECTstd` is `m · √v / √(-q / (2 a))` after strictly positive `asymDIFFUSION` and strictly positive time-dependent predictor variance; unstandardised `M` is not `TDPREDEFFECTstd`; `TIPREDEFFECTstd` is not `TDPREDEFFECTstd` even when `M = B`; intercept-style `A^{-1}[e^{A Δt} − I] M · √v / √p` is not `TDPREDEFFECTstd`; `m · √v / √(trait + p + added)` uses `TRAITVAR` and is not `TDPREDEFFECTstd`; Table 3 / p. 16 `T0TDPREDEFFECTstd` is `t0_m · √v / √p_0` after strictly positive free `T0VAR` and strictly positive TD predictor variance; unstandardised `t0_m` is not `T0TDPREDEFFECTstd`; `TDPREDEFFECTstd` uses `asymDIFFUSION` and is not `T0TDPREDEFFECTstd`; `T0TIPREDEFFECTstd` is not `T0TDPREDEFFECTstd` even when `t0_m = t0_b`; `t0_m · √v / √(trait + p_0 + added)` uses `TRAITVAR` and is not `T0TDPREDEFFECTstd`; free `T0VAR` does not require `a < 0`; p. 16 `T0VARstd` is `p_0 / p_0 = 1` after strictly positive free `T0VAR` (`solve(sqrt(diag(T0VAR))) %&% T0VAR`; OpenMx `%&%` is `t(A) %*% B %*% A`; default ridge is 0); unstandardised `T0VAR` is not `T0VARstd`; `T0TDPREDEFFECTstd` is not `T0VARstd`; `addedT0TIPREDVAR` is not `T0VARstd`; p. 16 `TRAITVARstd` is `trait / trait = 1` after strictly positive `TRAITVAR` (`solve(sqrt(diag(TRAITVAR))) %&% TRAITVAR`; OpenMx `%&%` is `t(A) %*% B %*% A`; no ridge addend); unstandardised `TRAITVAR` is not `TRAITVARstd`; `T0VARstd` is not `TRAITVARstd` even when both equal 1; `addedT0TIPREDVAR` is not `TRAITVARstd`; p. 16 `MANIFESTTRAITVARstd` is `ψ / ψ = 1` after strictly positive `MANIFESTTRAITVAR` (`solve(sqrt(diag(MANIFESTTRAITVAR))) %&% MANIFESTTRAITVAR`; OpenMx `%&%` is `t(A) %*% B %*% A`; 2017-era source adds ridging; default ridge is 0); unstandardised `MANIFESTTRAITVAR` is not `MANIFESTTRAITVARstd`; `TRAITVARstd` is not `MANIFESTTRAITVARstd` even when both equal 1; `MANIFESTVAR` is not `MANIFESTTRAITVARstd`; p. 16 `MANIFESTVARstd` is `θ / θ = 1` after strictly positive `MANIFESTVAR` (`solve(sqrt(diag(MANIFESTVAR))) %&% MANIFESTVAR`; OpenMx `%&%` is `t(A) %*% B %*% A`; 2017-era source adds ridging; default ridge is 0; 2017-era `dimnames` assignment to `latentNames` is a source bug); unstandardised `MANIFESTVAR` is not `MANIFESTVARstd`; `MANIFESTTRAITVARstd` is not `MANIFESTVARstd` even when both equal 1; Equation 5 `Var(y)` is not `MANIFESTVARstd`; p. 16 `TIPREDVARstd` is `v / v = 1` after strictly positive `TIPREDVAR` (`solve(sqrt(diag(TIPREDVAR))) %&% TIPREDVAR`; OpenMx `%&%` is `t(A) %*% B %*% A`; 2017-era source adds ridging; default ridge is 0; `dimnames` are `TIpredNames`); unstandardised `TIPREDVAR` is not `TIPREDVARstd`; `MANIFESTVARstd` is not `TIPREDVARstd` even when both equal 1; §7.2 `addedTIPREDVAR` is not `TIPREDVARstd`; p. 16 `asymDIFFUSIONstd` is `p / p = 1` after strictly positive `asymDIFFUSION` (`solve(sqrt(diag(asymDIFFUSION))) %&% asymDIFFUSION`; OpenMx `%&%` is `t(A) %*% B %*% A`; 2017-era source adds ridging; default ridge is 0; `dimnames` are `latentNames`); unstandardised `asymDIFFUSION` is not `asymDIFFUSIONstd`; `TIPREDVARstd` is not `asymDIFFUSIONstd` even when both equal 1; `DIFFUSIONstd` `−2 a` is not `asymDIFFUSIONstd`; p. 16 `discreteCINTstd` is `A^{-1}[e^{A Δt} − I] κ / √p` after strictly positive `asymDIFFUSION`; unstandardised `discreteCINT` is not `discreteCINTstd`; `κ / √p` is not `discreteCINTstd`; `(-κ / a) / √p` is not `discreteCINTstd`; `asymCINTstd` is `(-κ / a) / √p` after strictly positive `asymDIFFUSION`; unstandardised `asymCINT` is not `asymCINTstd`; `κ / √p` is not `asymCINTstd`; `discreteCINTstd` is not `asymCINTstd`; `T0MEANSstd` is `μ_0 / √p_0` after strictly positive free `T0VAR`; unstandardised `T0MEANS` is not `T0MEANSstd`; `T0VARstd` is not `T0MEANSstd`; `μ_0 / √asymDIFFUSION` is not `T0MEANSstd`; `MANIFESTMEANSstd` is `τ / √θ` after strictly positive `MANIFESTVAR`; unstandardised `MANIFESTMEANS` is not `MANIFESTMEANSstd`; `MANIFESTVARstd` is not `MANIFESTMEANSstd`; `τ / √(λ² Var(η) + θ)` is not `MANIFESTMEANSstd`; p. 16 `CINTstd` is `κ / √p` after strictly positive `asymDIFFUSION`; unstandardised `CINT` is not `CINTstd`; `asymCINTstd` is not `CINTstd`; `discreteCINTstd` is not `CINTstd`; `κ / √(trait + p + added)` is not `CINTstd`;))))), irregular already-centered residual lag, and strong/strict-gated latent means on the stacked psychometric PR (two-observation residual variance is identically `0` and caps at strong/scalar; Putnick & Bornstein, 2016, PMC5145197 opened 2026-08-19T22:15Z); full ESEM/DSEM remaining | partial | | CPU bounded multithreading + GPU/VRAM streaming/parity | ADR 0001/0006 | future `compute_backend` | accepted-target | | CPU bounded multithreading + GPU/VRAM streaming/parity | ADR 0001/0006 | future `compute_backend` | accepted-target | diff --git a/docs/research/compositional-cluster-pair-gates.md b/docs/research/compositional-cluster-pair-gates.md index 815eb9453..41b3a19a8 100644 --- a/docs/research/compositional-cluster-pair-gates.md +++ b/docs/research/compositional-cluster-pair-gates.md @@ -6,8 +6,10 @@ Euclidean coordinates and scores recovered cluster assignments with label-invariant pair precision and recall against known truth. -This slice does not fit a graphical lasso, run Leiden clustering, or claim -that a topic cluster is a causal construct. +This slice does not fit a graphical lasso or claim that a topic cluster is +a causal construct. Repeated Leiden partitions now live in +`docs/research/leiden-consensus.md`; they do **not** authorize Euclidean +distances on raw topic proportions. ## Authority diff --git a/docs/research/leiden-consensus.md b/docs/research/leiden-consensus.md new file mode 100644 index 000000000..61737c903 --- /dev/null +++ b/docs/research/leiden-consensus.md @@ -0,0 +1,70 @@ +# Leiden consensus partitions (GAP-009 remainder) + +## Scope + +This slice replaces the greedy union-find stand-in inside +`network_analysis::consensus` with Traag, Waltman, and van Eck (2019) +Leiden community detection on admitted positive topic–topic edges: + +1. fast local moving of nodes under Newman–Girvan modularity with γ = 1; +2. a refinement step that starts from singletons inside each community + and forbids internally disconnected communities; +3. aggregation of well-connected communities before the next local-move + pass; +4. co-assignment consensus across independently perturbed replicates + (Monti, Tamayo, Mesirov, & Golub, 2003; Hennig, 2007) with an explicit caller-supplied + `edge_drop_probability`. + +Isolated topics stay unclustered. Negative-effect edges never enter the +partition graph. This slice does not fit a graphical lasso, does not +claim that a topic community is a causal construct, and does not add an +export workflow. + +## Authoritative sources + +Traag, V. A., Waltman, L., & van Eck, N. J. (2019). From Louvain to +Leiden: Guaranteeing well-connected communities. *Scientific Reports, 9*, +Article 5233. https://doi.org/10.1038/s41598-019-41695-z + +Monti, S., Tamayo, P., Mesirov, J., & Golub, T. (2003). Consensus clustering: A resampling-based method for +class discovery and visualization of gene expression microarray data. +*Machine Learning, 52*(1–2), 91–118. +https://doi.org/10.1023/A:1023949509487 + +Hennig, C. (2007). Cluster-wise assessment of cluster stability. +*Computational Statistics & Data Analysis, 52*(1), 258–271. +https://doi.org/10.1016/j.csda.2006.11.025 + +Newman, M. E. J., & Girvan, M. (2004). Finding and evaluating community +structure in networks. *Physical Review E, 69*(2), 026113. +https://doi.org/10.1103/PhysRevE.69.026113 + +## Formula notes + +- Undirected edge weight `w_ij` is the admitted positive effect. Self-loops + and non-positive effects are dropped. Duplicate pairs are summed. +- Strength `k_i = Σ_j w_ij`. Total weight `m = (Σ_i k_i) / 2`. +- Modularity gain of moving node `i` from community `c` to `d` at γ = 1: + `ΔQ = (k_{i→d} − k_{i→c}) / m − k_i (Σ_d − Σ_c + k_i) / (2 m²)`, + where `Σ_c` still includes `k_i`. +- Louvain is rejected: it can emit internally disconnected communities + (Traag et al., 2019). Union-find is rejected: it merges every surviving + edge into one component and makes no modularity claim. + +## Verification + +- empty and edgeless graphs stay singletons; +- two nodes with a positive edge share a community; +- two triangles joined by a weak bridge stay two communities (union-find + would glue them) both at the Leiden partition and at the consensus + wrapper with zero drop probability; +- out-of-range, self, and non-positive edges are ignored; +- identical seeds reproduce the partition; +- isolated topics keep a distinct unclustered assignment at consensus; +- aggregation sums parallel supernode edges and keeps distinct neighbours; +- the connectedness oracle rejects a disconnected community and accepts + singletons; +- aggregation continues when nodes moved even if the aggregate did not + shrink; +- out-of-range and self-loop endpoints stay unclustered at consensus; +- permuting equal edges does not change the consensus assignment. diff --git a/docs/research/standards-and-literature.md b/docs/research/standards-and-literature.md index 8ce0ecc77..e61b73579 100644 --- a/docs/research/standards-and-literature.md +++ b/docs/research/standards-and-literature.md @@ -158,9 +158,11 @@ Fisher, R. A. (1921). On the "probable error" of a coefficient of correlation de Friedman, J., Hastie, T., & Tibshirani, R. (2008). Sparse inverse covariance estimation with the graphical lasso. *Biostatistics, 9*(3), 432–441. https://doi.org/10.1093/biostatistics/kxm045 -Hennig, C. (2007). Cluster-wise assessment of cluster stability. *Computational Statistics & Data Analysis, 52*(1), 258–281. https://doi.org/10.1016/j.csda.2006.11.025 +Hennig, C. (2007). Cluster-wise assessment of cluster stability. *Computational Statistics & Data Analysis, 52*(1), 258–271. https://doi.org/10.1016/j.csda.2006.11.025 -Monti, S. (2003). Consensus clustering: A resampling-based method for class discovery and visualization of gene expression microarray data. *Machine Learning, 52*(1–2), 91–118. https://doi.org/10.1023/A:1023949509487 +Monti, S., Tamayo, P., Mesirov, J., & Golub, T. (2003). Consensus clustering: A resampling-based method for class discovery and visualization of gene expression microarray data. *Machine Learning, 52*(1–2), 91–118. https://doi.org/10.1023/A:1023949509487 + +Newman, M. E. J., & Girvan, M. (2004). Finding and evaluating community structure in networks. *Physical Review E, 69*(2), 026113. https://doi.org/10.1103/PhysRevE.69.026113 Traag, V. A., Waltman, L., & van Eck, N. J. (2019). From Louvain to Leiden: Guaranteeing well-connected communities. *Scientific Reports, 9*, Article 5233. https://doi.org/10.1038/s41598-019-41695-z @@ -180,9 +182,15 @@ credible intervals and selection fractions are percentile bootstrap quantities over posterior draws (Efron, 1979). Consensus clusters come from resampled co-assignment matrices with an explicit, caller-supplied perturbation probability rather than an implicit constant (Monti, 2003; -Hennig, 2007). No numeric threshold or weight in this path is chosen by -heuristic; each is either an explicit parameter with stated provenance -or a value derived from these primary sources. +Hennig, 2007). Each co-assignment replicate is a Traag, Waltman, and van +Eck (2019) Leiden partition of the surviving positive edges under +Newman–Girvan modularity with γ = 1; Louvain is rejected because it can +emit internally disconnected communities, and union-find is rejected +because it merges every surviving edge into one component and makes no +modularity claim. Isolated topics remain unclustered. See +`docs/research/leiden-consensus.md`. No numeric threshold or weight in +this path is chosen by heuristic; each is either an explicit parameter +with stated provenance or a value derived from these primary sources. ## Time, events, and topic detection and tracking