facett_graphview/analysis/

centrality.rs

//! **Centrality** — who matters in the graph, five ways. All return a `Vec<f32>`
//! indexed by node (`0..n`). Deterministic; no RNG.

use std::collections::VecDeque;

use super::Adjacency;

/// **Degree centrality** — undirected degree normalised by `n - 1` (the fraction of
/// other nodes a node is adjacent to). `0.0` for `n ≤ 1`.
#[must_use]
pub fn degree(g: &Adjacency) -> Vec<f32> {
    let denom = (g.n.saturating_sub(1)).max(1) as f32;
    (0..g.n).map(|i| g.degree(i) as f32 / denom).collect()
}

/// **In-degree** and **out-degree** centrality (directed), normalised by `n - 1`.
#[must_use]
pub fn in_out_degree(g: &Adjacency) -> (Vec<f32>, Vec<f32>) {
    let denom = (g.n.saturating_sub(1)).max(1) as f32;
    let ins = (0..g.n).map(|i| g.in_degree(i) as f32 / denom).collect();
    let outs = (0..g.n).map(|i| g.out_degree(i) as f32 / denom).collect();
    (ins, outs)
}

/// **Closeness centrality** (undirected) — `(reachable-1) / Σ dist`, scaled by the
/// reachable fraction so nodes in small components aren't over-rated (Wasserman–Faust).
/// `0.0` for isolated nodes. BFS from each node (unweighted).
#[must_use]
pub fn closeness(g: &Adjacency) -> Vec<f32> {
    let n = g.n;
    (0..n)
        .map(|s| {
            let dist = bfs_dist(g, s);
            let mut sum = 0.0f32;
            let mut reach = 0usize;
            for (i, &d) in dist.iter().enumerate() {
                if i != s && d >= 0 {
                    sum += d as f32;
                    reach += 1;
                }
            }
            if sum <= 0.0 || n <= 1 {
                0.0
            } else {
                // (reach / sum) · (reach / (n-1)) — harmonic-ish, component-scaled.
                (reach as f32 / sum) * (reach as f32 / (n - 1) as f32)
            }
        })
        .collect()
}

/// **Betweenness centrality** (Brandes' algorithm, undirected, unweighted). The
/// fraction of shortest paths through each node; normalised by `(n-1)(n-2)` so it
/// lands in `[0, 1]`. `O(n·m)`.
#[must_use]
pub fn betweenness(g: &Adjacency) -> Vec<f32> {
    let n = g.n;
    let mut bc = vec![0.0f32; n];
    for s in 0..n {
        // Single-source shortest-path counting (Brandes).
        let mut stack = Vec::new();
        let mut pred: Vec<Vec<usize>> = vec![Vec::new(); n];
        let mut sigma = vec![0.0f32; n];
        let mut dist = vec![-1i32; n];
        sigma[s] = 1.0;
        dist[s] = 0;
        let mut q = VecDeque::new();
        q.push_back(s);
        while let Some(v) = q.pop_front() {
            stack.push(v);
            for &(w, _) in &g.und[v] {
                if dist[w] < 0 {
                    dist[w] = dist[v] + 1;
                    q.push_back(w);
                }
                if dist[w] == dist[v] + 1 {
                    sigma[w] += sigma[v];
                    pred[w].push(v);
                }
            }
        }
        // Accumulate dependencies back-to-front.
        let mut delta = vec![0.0f32; n];
        while let Some(w) = stack.pop() {
            for &v in &pred[w] {
                if sigma[w] > 0.0 {
                    delta[v] += (sigma[v] / sigma[w]) * (1.0 + delta[w]);
                }
            }
            if w != s {
                bc[w] += delta[w];
            }
        }
    }
    // Undirected: each pair counted twice → halve. Normalise by (n-1)(n-2).
    let norm = if n > 2 { ((n - 1) * (n - 2)) as f32 } else { 1.0 };
    for x in &mut bc {
        *x = (*x * 0.5) / norm;
    }
    bc
}

/// **PageRank** (directed) with damping `d` (typically `0.85`) for `iters` power
/// iterations. Dangling nodes (no out-edges) redistribute their mass uniformly.
/// Returns a probability distribution summing to ≈1.
#[must_use]
pub fn pagerank(g: &Adjacency, d: f32, iters: usize) -> Vec<f32> {
    let n = g.n;
    if n == 0 {
        return Vec::new();
    }
    let base = (1.0 - d) / n as f32;
    let mut rank = vec![1.0f32 / n as f32; n];
    // Out-weight sums for normalising contributions.
    let out_sum: Vec<f32> = (0..n).map(|i| g.out[i].iter().map(|(_, w)| *w).sum::<f32>()).collect();
    for _ in 0..iters {
        let mut next = vec![base; n];
        // Dangling mass (nodes with no out-edges) spread uniformly.
        let dangling: f32 = (0..n).filter(|&i| g.out[i].is_empty()).map(|i| rank[i]).sum();
        let dangle_share = d * dangling / n as f32;
        for v in &mut next {
            *v += dangle_share;
        }
        for i in 0..n {
            if out_sum[i] <= 0.0 {
                continue;
            }
            for &(j, w) in &g.out[i] {
                next[j] += d * rank[i] * (w / out_sum[i]);
            }
        }
        rank = next;
    }
    // Renormalise against tiny drift.
    let total: f32 = rank.iter().sum();
    if total > 0.0 {
        for v in &mut rank {
            *v /= total;
        }
    }
    rank
}

/// **Eigenvector centrality** (undirected) via power iteration for `iters` steps.
/// Iterates `(A + I)` rather than `A` — the identity shift makes every eigenvalue
/// positive, curing the **bipartite oscillation** that traps plain `A·x` (a star /
/// path never converges because `+λ` and `−λ` have equal magnitude). L2-normalised
/// each step. `0.0` everywhere for an edgeless graph (centrality is undefined there).
#[must_use]
pub fn eigenvector(g: &Adjacency, iters: usize) -> Vec<f32> {
    let n = g.n;
    if n == 0 {
        return Vec::new();
    }
    // Edgeless ⇒ no meaningful centrality.
    if (0..n).all(|i| g.und[i].is_empty()) {
        return vec![0.0; n];
    }
    let mut x = vec![1.0f32 / (n as f32).sqrt(); n];
    for _ in 0..iters {
        // next = (A + I)·x — the +I shift breaks bipartite oscillation.
        let mut next = x.clone();
        for i in 0..n {
            for &(j, w) in &g.und[i] {
                next[i] += w * x[j];
            }
        }
        let norm = next.iter().map(|v| v * v).sum::<f32>().sqrt();
        if norm <= f32::EPSILON {
            return vec![0.0; n];
        }
        for v in &mut next {
            *v /= norm;
        }
        x = next;
    }
    // Fix the sign so the vector is non-negative (Perron–Frobenius).
    if x.iter().sum::<f32>() < 0.0 {
        for v in &mut x {
            *v = -*v;
        }
    }
    x
}

/// BFS distances from `s` over the **undirected** view; `-1` for unreachable.
fn bfs_dist(g: &Adjacency, s: usize) -> Vec<i32> {
    let mut dist = vec![-1i32; g.n];
    let mut q = VecDeque::new();
    dist[s] = 0;
    q.push_back(s);
    while let Some(v) = q.pop_front() {
        for &(w, _) in &g.und[v] {
            if dist[w] < 0 {
                dist[w] = dist[v] + 1;
                q.push_back(w);
            }
        }
    }
    dist
}

#[cfg(test)]
mod tests {
    use super::*;

    /// A 5-node path `0-1-2-3-4` — the classic betweenness/closeness fixture.
    fn path5() -> Adjacency {
        Adjacency::from_edges(5, &[(0, 1), (1, 2), (2, 3), (3, 4)])
    }

    /// A star: centre 0 connected to 1,2,3,4.
    fn star5() -> Adjacency {
        Adjacency::from_edges(5, &[(0, 1), (0, 2), (0, 3), (0, 4)])
    }

    #[test]
    fn degree_peaks_at_the_star_centre() {
        let d = degree(&star5());
        assert!((d[0] - 1.0).abs() < 1e-6, "the hub is adjacent to everyone");
        for i in 1..5 {
            assert!((d[i] - 0.25).abs() < 1e-6, "a leaf touches only the hub");
        }
    }

    #[test]
    fn betweenness_peaks_at_the_path_middle() {
        let b = betweenness(&path5());
        // The middle node (2) lies on the most shortest paths; the ends on none.
        assert!(b[2] > b[1] && b[1] > b[0], "betweenness rises toward the centre: {b:?}");
        assert!((b[0]).abs() < 1e-6 && (b[4]).abs() < 1e-6, "endpoints are on no through-path");
    }

    #[test]
    fn closeness_peaks_at_the_star_centre() {
        let c = closeness(&star5());
        for i in 1..5 {
            assert!(c[0] > c[i], "the hub is closest to all: {c:?}");
        }
    }

    #[test]
    fn pagerank_is_a_distribution_and_favours_a_sink() {
        // 1→0, 2→0, 3→0: everyone points at 0, so 0 collects the rank.
        let g = Adjacency::from_edges(4, &[(1, 0), (2, 0), (3, 0)]);
        let pr = pagerank(&g, 0.85, 50);
        let sum: f32 = pr.iter().sum();
        assert!((sum - 1.0).abs() < 1e-3, "pagerank sums to 1 (got {sum})");
        assert!(pr[0] > pr[1] && pr[0] > pr[2] && pr[0] > pr[3], "the sink ranks highest: {pr:?}");
    }

    #[test]
    fn eigenvector_peaks_at_the_star_centre() {
        let e = eigenvector(&star5(), 100);
        for i in 1..5 {
            assert!(e[0] > e[i], "the hub dominates the eigenvector: {e:?}");
        }
        // Non-negative (Perron–Frobenius sign fix).
        assert!(e.iter().all(|&v| v >= -1e-6));
    }

    #[test]
    fn edgeless_graphs_do_not_panic() {
        let g = Adjacency::from_edges(3, &[]);
        assert_eq!(betweenness(&g), vec![0.0; 3]);
        assert_eq!(closeness(&g), vec![0.0; 3]);
        assert_eq!(eigenvector(&g, 10), vec![0.0; 3]);
        let pr = pagerank(&g, 0.85, 10);
        assert!((pr.iter().sum::<f32>() - 1.0).abs() < 1e-3);
    }
}