modewise_stats/
nbs.rs

1//! Network-based statistic: connected components of suprathreshold edges (Zalesky et al. 2010).
2
3/// Where an edge feature sits: which network (component k of the connectivity) and which nodes.
4#[derive(Debug, Clone, Copy, PartialEq, Eq)]
5pub struct Edge {
6    pub network: usize,
7    pub a: usize,
8    pub b: usize,
9}
10
11/// A connected set of suprathreshold edges within one network.
12#[derive(Debug, Clone, PartialEq)]
13pub struct Component {
14    pub network: usize,
15    /// Indices of the edge features in the component.
16    pub edges: Vec<usize>,
17    /// Permutation p of the component's size against the maximum size under the null.
18    pub p: f64,
19}
20
21/// Components (as lists of edge feature indices) among edges whose statistic reaches `threshold`.
22pub(crate) fn components(
23    stat: &[f64],
24    edges: &[Option<Edge>],
25    threshold: f64,
26) -> Vec<(usize, Vec<usize>)> {
27    let nodes = edges
28        .iter()
29        .flatten()
30        .map(|e| e.a.max(e.b) + 1)
31        .max()
32        .unwrap_or(0);
33    let networks = edges
34        .iter()
35        .flatten()
36        .map(|e| e.network + 1)
37        .max()
38        .unwrap_or(0);
39    let mut parent: Vec<usize> = (0..nodes * networks).collect();
40    fn find(parent: &mut [usize], mut x: usize) -> usize {
41        while parent[x] != x {
42            parent[x] = parent[parent[x]];
43            x = parent[x];
44        }
45        x
46    }
47    let supra: Vec<(usize, Edge)> = edges
48        .iter()
49        .enumerate()
50        .filter_map(|(f, e)| e.filter(|_| stat[f] >= threshold).map(|e| (f, e)))
51        .collect();
52    for (_, e) in &supra {
53        let (x, y) = (
54            find(&mut parent, e.network * nodes + e.a),
55            find(&mut parent, e.network * nodes + e.b),
56        );
57        parent[x] = y;
58    }
59    let mut groups: std::collections::BTreeMap<usize, (usize, Vec<usize>)> =
60        std::collections::BTreeMap::new();
61    for (f, e) in supra {
62        let root = find(&mut parent, e.network * nodes + e.a);
63        groups
64            .entry(root)
65            .or_insert((e.network, Vec::new()))
66            .1
67            .push(f);
68    }
69    groups.into_values().collect()
70}
71
72/// The size (in edges) of the largest component.
73pub(crate) fn largest(stat: &[f64], edges: &[Option<Edge>], threshold: f64) -> usize {
74    components(stat, edges, threshold)
75        .iter()
76        .map(|(_, e)| e.len())
77        .max()
78        .unwrap_or(0)
79}
80
81#[cfg(test)]
82mod tests {
83    use super::*;
84
85    #[test]
86    fn components_follow_shared_nodes_within_a_network() {
87        let e = |network, a, b| Some(Edge { network, a, b });
88        // Network 0: 0-1, 1-2 (joined), 3-4 (separate); network 1: 0-1 (separate network).
89        let edges = [
90            e(0, 0, 1),
91            e(0, 1, 2),
92            e(0, 3, 4),
93            e(1, 0, 1),
94            None,
95            e(0, 2, 3),
96        ];
97        let stat = [3.5, 4.0, 3.1, 5.0, 9.0, 1.0];
98        let mut found = components(&stat, &edges, 3.0);
99        found.sort_by_key(|(_, e)| e.clone());
100        assert_eq!(found, [(0, vec![0, 1]), (0, vec![2]), (1, vec![3])]);
101        assert_eq!(largest(&stat, &edges, 3.0), 2);
102        let bridged = [3.5, 4.0, 3.1, 5.0, 9.0, 3.0];
103        assert_eq!(
104            largest(&bridged, &edges, 3.0),
105            4,
106            "edge 2-3 joins both into one"
107        );
108        assert_eq!(largest(&stat, &edges, 10.0), 0);
109        assert_eq!(largest(&[], &[], 1.0), 0);
110    }
111}
modewise_stats/nbs.rs

modewise_stats/
nbs.rs
