Skip to main content

weavatrix_graph/algo/network/centrality/
spectral.rs

1#![allow(clippy::cast_precision_loss)]
2
3use super::super::adjacency::{SlotAdjacency, adjacency};
4use super::math::square_root;
5use crate::algo::traversal::Direction;
6use crate::{GraphError, IndexGraphView, Result, String, Vec};
7
8#[derive(Debug, Clone, PartialEq)]
9pub struct IterativeCentrality<Node> {
10    scores: Vec<(Node, f64)>,
11    iterations: usize,
12    converged: bool,
13}
14
15impl<Node> IterativeCentrality<Node> {
16    #[must_use]
17    pub fn scores(&self) -> &[(Node, f64)] {
18        &self.scores
19    }
20
21    #[must_use]
22    pub const fn iterations(&self) -> usize {
23        self.iterations
24    }
25
26    #[must_use]
27    pub const fn converged(&self) -> bool {
28        self.converged
29    }
30
31    #[must_use]
32    pub fn into_scores(self) -> Vec<(Node, f64)> {
33        self.scores
34    }
35}
36
37/// Computes eigenvector centrality by shifted power iteration.
38///
39/// # Errors
40///
41/// Returns an error for zero iterations or an invalid tolerance.
42pub fn eigenvector_centrality<G>(
43    graph: &G,
44    max_iterations: usize,
45    tolerance: f64,
46) -> Result<IterativeCentrality<G::Node>>
47where
48    G: IndexGraphView,
49{
50    validate(max_iterations, tolerance)?;
51    let adjacent = adjacency(graph, Direction::Outgoing);
52    let initial = if adjacent.nodes.is_empty() {
53        0.0
54    } else {
55        1.0 / square_root(adjacent.nodes.len() as f64)
56    };
57    let mut current = vec![initial; graph.node_bound()];
58    let mut next = vec![0.0; graph.node_bound()];
59    for iteration in 1..=max_iterations {
60        next.fill(0.0);
61        for &node in &adjacent.nodes {
62            let source = G::node_slot(node);
63            next[source] += current[source];
64            for &target in &adjacent.neighbors[source] {
65                next[target] += current[source];
66            }
67        }
68        normalize::<G>(&mut next, &adjacent);
69        if converged::<G>(&current, &next, &adjacent, tolerance) {
70            return Ok(result::<G>(&adjacent, &next, iteration, true));
71        }
72        core::mem::swap(&mut current, &mut next);
73    }
74    Ok(result::<G>(&adjacent, &current, max_iterations, false))
75}
76
77/// Computes Katz centrality by fixed-point iteration.
78///
79/// # Errors
80///
81/// Returns an error for invalid iteration, tolerance, alpha, or beta values.
82pub fn katz_centrality<G>(
83    graph: &G,
84    alpha: f64,
85    beta: f64,
86    max_iterations: usize,
87    tolerance: f64,
88) -> Result<IterativeCentrality<G::Node>>
89where
90    G: IndexGraphView,
91{
92    validate(max_iterations, tolerance)?;
93    if !alpha.is_finite() || alpha < 0.0 || !beta.is_finite() {
94        return Err(invalid_parameter("alpha/beta"));
95    }
96    let adjacent = adjacency(graph, Direction::Outgoing);
97    let mut current = vec![1.0; graph.node_bound()];
98    let mut next = vec![beta; graph.node_bound()];
99    for iteration in 1..=max_iterations {
100        next.fill(beta);
101        for &node in &adjacent.nodes {
102            let source = G::node_slot(node);
103            for &target in &adjacent.neighbors[source] {
104                next[target] += alpha * current[source];
105            }
106        }
107        if next.iter().any(|score| !score.is_finite()) {
108            return Err(invalid_parameter("alpha"));
109        }
110        if converged::<G>(&current, &next, &adjacent, tolerance) {
111            return Ok(result::<G>(&adjacent, &next, iteration, true));
112        }
113        core::mem::swap(&mut current, &mut next);
114    }
115    Ok(result::<G>(&adjacent, &current, max_iterations, false))
116}
117
118fn normalize<G>(scores: &mut [f64], adjacent: &SlotAdjacency<G::Node>)
119where
120    G: IndexGraphView,
121{
122    let norm = adjacent
123        .nodes
124        .iter()
125        .map(|node| {
126            let value = scores[G::node_slot(*node)];
127            value * value
128        })
129        .sum::<f64>();
130    let norm = square_root(norm);
131    if norm > 0.0 {
132        for node in &adjacent.nodes {
133            scores[G::node_slot(*node)] /= norm;
134        }
135    }
136}
137
138fn converged<G>(
139    current: &[f64],
140    next: &[f64],
141    adjacent: &SlotAdjacency<G::Node>,
142    tolerance: f64,
143) -> bool
144where
145    G: IndexGraphView,
146{
147    adjacent
148        .nodes
149        .iter()
150        .map(|node| (next[G::node_slot(*node)] - current[G::node_slot(*node)]).abs())
151        .sum::<f64>()
152        <= tolerance * adjacent.nodes.len() as f64
153}
154
155fn result<G>(
156    adjacent: &SlotAdjacency<G::Node>,
157    scores: &[f64],
158    iterations: usize,
159    converged: bool,
160) -> IterativeCentrality<G::Node>
161where
162    G: IndexGraphView,
163{
164    IterativeCentrality {
165        scores: adjacent
166            .nodes
167            .iter()
168            .copied()
169            .map(|node| (node, scores[G::node_slot(node)]))
170            .collect(),
171        iterations,
172        converged,
173    }
174}
175
176fn validate(max_iterations: usize, tolerance: f64) -> Result<()> {
177    if max_iterations == 0 || !tolerance.is_finite() || tolerance <= 0.0 {
178        return Err(invalid_parameter("max_iterations/tolerance"));
179    }
180    Ok(())
181}
182
183fn invalid_parameter(parameter: &'static str) -> GraphError {
184    GraphError::InvalidAlgorithmParameter {
185        algorithm: "centrality",
186        parameter,
187        value: String::from("must be finite, positive, and non-zero where required"),
188    }
189}