xgraph/core/algo/

centrality.rs

//! Centrality measures, written once against [`GraphView`].
//!
//! Unifying the two previous copies also settles where they disagreed. The heterogeneous
//! versions of Brandes' algorithm and of closeness mixed edge weights into a breadth-first
//! sweep: nodes came off the queue in hop order while distances accumulated weights, so
//! with non-uniform weights the shortest-path counts were simply wrong. Both measures are
//! defined here on hop counts, which is what "betweenness" and "closeness" mean when no
//! distance metric is named, and matches the homogeneous implementations.
//!
//! The iterative measures — eigenvector, PageRank and Katz — do use weights, and validate
//! their parameters instead of quietly returning a diverged vector.

use crate::core::algo::connectivity::Connectivity;
use crate::core::algo::traversal::TraversalError;
use crate::core::view::{GraphView, NodeId};
use std::collections::{HashMap, VecDeque};

/// Error type for centrality computation failures.
#[derive(Debug, Clone, PartialEq, thiserror::Error)]
pub enum CentralityError {
    /// A node referenced during the computation does not exist.
    #[error("Invalid node reference: node ID {0} not found")]
    InvalidNode(NodeId),
    /// A parameter fell outside the range the measure is defined for.
    #[error("Parameter `{name}` is {value}, expected {expected}")]
    InvalidParameter {
        /// The parameter's name.
        name: &'static str,
        /// The value supplied.
        value: f64,
        /// What the measure requires.
        expected: &'static str,
    },
}

impl From<TraversalError> for CentralityError {
    fn from(error: TraversalError) -> Self {
        match error {
            TraversalError::InvalidNode(id) => CentralityError::InvalidNode(id),
        }
    }
}

/// Centrality measures available on every [`GraphView`].
pub trait Centrality: GraphView {
    /// Degree centrality: the number of edges incident to each node.
    ///
    /// Parallel edges count separately, so a node joined to one neighbour by three edges
    /// has degree three.
    fn degree_centrality(&self) -> HashMap<NodeId, usize> {
        self.node_id_vec()
            .into_iter()
            .map(|node| (node, self.degree(node)))
            .collect()
    }

    /// Betweenness centrality: how often each node lies on a shortest path between others.
    ///
    /// Computed with Brandes' algorithm over hop counts, normalized so that the maximum
    /// possible value is 1. Graphs of fewer than three nodes have no intermediate positions
    /// to occupy, so every score is zero.
    fn betweenness_centrality(&self) -> Result<HashMap<NodeId, f64>, CentralityError> {
        let nodes = self.node_id_vec();
        let count = nodes.len();

        // With fewer than three nodes there is no path with an interior, and the usual
        // normalisation constant would be a division by zero.
        if count < 3 {
            return Ok(nodes.into_iter().map(|node| (node, 0.0)).collect());
        }

        // Node IDs may be sparse, so the per-source arrays are indexed by position.
        let index_of: HashMap<NodeId, usize> = nodes
            .iter()
            .enumerate()
            .map(|(index, &id)| (id, index))
            .collect();

        let contributions = nodes
            .iter()
            .map(|&source| brandes_from(self, source, &nodes, &index_of))
            .collect::<Result<Vec<_>, _>>()?;

        let betweenness = accumulate_betweenness(&nodes, contributions);

        // On an undirected graph each pair is walked from both ends.
        let pairs = ((count - 1) * (count - 2)) as f64;
        let scale = if self.is_directed() {
            1.0 / pairs
        } else {
            2.0 / pairs
        };
        let halve = if self.is_directed() { 1.0 } else { 0.5 };

        Ok(normalize_betweenness(betweenness, halve * scale))
    }

    /// Betweenness centrality, with the per-source sweeps shared across cores.
    ///
    /// Identical in result to [`betweenness_centrality`](Centrality::betweenness_centrality)
    /// — each source contributes independently, so splitting the work changes only how long
    /// it takes. Worth reaching for on large graphs, where the serial version is
    /// O(|V| · |E|).
    ///
    /// Offered as a separate method rather than switching the existing one on a feature:
    /// cargo unifies features across a dependency graph, so a feature that changed this
    /// method's bounds would break unrelated crates that happen to share the build.
    #[cfg(feature = "parallel")]
    fn par_betweenness_centrality(&self) -> Result<HashMap<NodeId, f64>, CentralityError>
    where
        Self: Sync,
    {
        use rayon::prelude::*;

        let nodes = self.node_id_vec();
        let count = nodes.len();
        if count < 3 {
            return Ok(nodes.into_iter().map(|node| (node, 0.0)).collect());
        }

        let index_of: HashMap<NodeId, usize> = nodes
            .iter()
            .enumerate()
            .map(|(index, &id)| (id, index))
            .collect();

        let contributions = nodes
            .par_iter()
            .map(|&source| brandes_from(self, source, &nodes, &index_of))
            .collect::<Result<Vec<_>, _>>()?;

        let betweenness = accumulate_betweenness(&nodes, contributions);
        let pairs = ((count - 1) * (count - 2)) as f64;
        let scale = if self.is_directed() {
            1.0 / pairs
        } else {
            2.0 / pairs
        };
        let halve = if self.is_directed() { 1.0 } else { 0.5 };
        Ok(normalize_betweenness(betweenness, halve * scale))
    }

    /// Closeness centrality: how near each node is to all the nodes it can reach.
    ///
    /// Measured in hops, as the number of reachable nodes divided by the total distance to
    /// them. A node that reaches nothing scores zero.
    fn closeness_centrality(&self) -> Result<HashMap<NodeId, f64>, CentralityError> {
        Ok(self
            .node_id_vec()
            .into_iter()
            .map(|source| (source, closeness_from(self, source)))
            .collect())
    }

    /// Closeness centrality, with the per-source sweeps shared across cores.
    ///
    /// Identical in result to [`closeness_centrality`](Centrality::closeness_centrality).
    #[cfg(feature = "parallel")]
    fn par_closeness_centrality(&self) -> Result<HashMap<NodeId, f64>, CentralityError>
    where
        Self: Sync,
    {
        use rayon::prelude::*;

        Ok(self
            .node_id_vec()
            .par_iter()
            .map(|&source| (source, closeness_from(self, source)))
            .collect())
    }

    /// Eigenvector centrality: a node is important when its neighbours are important.
    ///
    /// Computed by power iteration, stopping once the largest change between rounds falls
    /// below `tolerance` or `max_iter` rounds have run. The result is scaled to sum to 1.
    fn eigenvector_centrality(
        &self,
        max_iter: usize,
        tolerance: f64,
    ) -> Result<HashMap<NodeId, f64>, CentralityError>
    where
        Self::Weight: Into<f64>,
    {
        check_non_negative("tolerance", tolerance)?;

        let nodes = self.node_id_vec();
        let mut centrality: HashMap<NodeId, f64> = nodes.iter().map(|&node| (node, 1.0)).collect();
        if nodes.is_empty() {
            return Ok(centrality);
        }

        for _ in 0..max_iter {
            let mut next: HashMap<NodeId, f64> = HashMap::with_capacity(nodes.len());
            let mut largest_change = 0.0_f64;

            for &node in &nodes {
                let mut score = 1.0;
                for link in self.links(node) {
                    let weight: f64 = link.weight.into();
                    score += weight * centrality.get(&link.target).copied().unwrap_or(0.0);
                }
                largest_change = largest_change.max((score - centrality[&node]).abs());
                next.insert(node, score);
            }

            let norm = next.values().map(|value| value * value).sum::<f64>().sqrt();
            if norm > 0.0 {
                for value in next.values_mut() {
                    *value /= norm;
                }
            }
            centrality = next;

            if largest_change < tolerance {
                break;
            }
        }

        normalize_to_sum(&mut centrality);
        Ok(centrality)
    }

    /// PageRank: the stationary distribution of a random walk that teleports.
    ///
    /// `damping` is the probability of following an edge rather than teleporting; 0.85 is
    /// the conventional value. Iteration stops once the total change between rounds falls
    /// below `tolerance`.
    ///
    /// The rank a node passes on is split among its outgoing edges in proportion to their
    /// weights. A node with no outgoing edges keeps nothing: its rank leaks out of the
    /// system, so the vector is renormalized at the end.
    fn pagerank(
        &self,
        damping: f64,
        max_iter: usize,
        tolerance: f64,
    ) -> Result<HashMap<NodeId, f64>, CentralityError>
    where
        Self::Weight: Into<f64>,
    {
        check_probability("damping", damping)?;
        check_non_negative("tolerance", tolerance)?;

        let nodes = self.node_id_vec();
        let count = nodes.len();
        if count == 0 {
            return Ok(HashMap::new());
        }

        let mut rank: HashMap<NodeId, f64> = nodes
            .iter()
            .map(|&node| (node, 1.0 / count as f64))
            .collect();

        // Total outgoing weight per node, computed once. The previous implementation
        // rediscovered a node's inbound links by scanning every other node on every
        // iteration, which made each round cost O(|V|^2) instead of O(|V| + |E|).
        let out_weight: HashMap<NodeId, f64> = nodes
            .iter()
            .map(|&node| {
                let total: f64 = self.links(node).map(|link| link.weight.into()).sum();
                (node, total)
            })
            .collect();

        let teleport = (1.0 - damping) / count as f64;

        for _ in 0..max_iter {
            let mut next: HashMap<NodeId, f64> =
                nodes.iter().map(|&node| (node, teleport)).collect();

            for &node in &nodes {
                let outgoing = out_weight[&node];
                if outgoing <= 0.0 {
                    continue;
                }
                let share = damping * rank[&node] / outgoing;
                for link in self.links(node) {
                    let weight: f64 = link.weight.into();
                    if let Some(slot) = next.get_mut(&link.target) {
                        *slot += share * weight;
                    }
                }
            }

            let change: f64 = nodes
                .iter()
                .map(|node| (next[node] - rank[node]).abs())
                .sum();
            rank = next;

            if change < tolerance {
                break;
            }
        }

        normalize_to_sum(&mut rank);
        Ok(rank)
    }

    /// Katz centrality: influence that decays with distance.
    ///
    /// Each node starts at `beta` and accumulates `alpha` times the score of each
    /// neighbour, weighted by the connecting edge. `alpha` must be small enough for the
    /// series to converge — below the reciprocal of the largest eigenvalue, for which the
    /// reciprocal of the maximum weighted degree is a safe bound.
    fn katz_centrality(
        &self,
        alpha: f64,
        beta: f64,
        max_iter: usize,
        tolerance: f64,
    ) -> Result<HashMap<NodeId, f64>, CentralityError>
    where
        Self::Weight: Into<f64>,
    {
        check_non_negative("tolerance", tolerance)?;
        if alpha <= 0.0 {
            return Err(CentralityError::InvalidParameter {
                name: "alpha",
                value: alpha,
                expected: "a positive attenuation factor",
            });
        }

        let nodes = self.node_id_vec();
        let mut centrality: HashMap<NodeId, f64> = nodes.iter().map(|&node| (node, beta)).collect();
        if nodes.is_empty() {
            return Ok(centrality);
        }

        for _ in 0..max_iter {
            let mut next: HashMap<NodeId, f64> = HashMap::with_capacity(nodes.len());
            let mut largest_change = 0.0_f64;

            for &node in &nodes {
                let mut score = beta;
                for link in self.links(node) {
                    let weight: f64 = link.weight.into();
                    score += alpha * weight * centrality.get(&link.target).copied().unwrap_or(0.0);
                }
                largest_change = largest_change.max((score - centrality[&node]).abs());
                next.insert(node, score);
            }

            centrality = next;
            if largest_change < tolerance {
                break;
            }
        }
        Ok(centrality)
    }

    /// The in-degree of each node: the number of edges pointing at it.
    ///
    /// On an undirected view this coincides with [`degree_centrality`](Centrality::degree_centrality).
    fn in_degree_centrality(&self) -> HashMap<NodeId, usize> {
        self.predecessors()
            .into_iter()
            .map(|(node, sources)| (node, sources.len()))
            .collect()
    }
}

impl<G: GraphView + ?Sized> Centrality for G {}

/// How near one node is to everything it can reach, in hops.
///
/// Kept as a free function so the serial and parallel drivers share one implementation.
fn closeness_from<G>(graph: &G, source: NodeId) -> f64
where
    G: GraphView + ?Sized,
{
    let mut distance: HashMap<NodeId, usize> = HashMap::from([(source, 0)]);
    let mut queue = VecDeque::from([source]);

    while let Some(node) = queue.pop_front() {
        let step = distance[&node] + 1;
        for link in graph.links(node) {
            if let std::collections::hash_map::Entry::Vacant(slot) = distance.entry(link.target) {
                slot.insert(step);
                queue.push_back(link.target);
            }
        }
    }

    let total: usize = distance.values().sum();
    let reachable = distance.len().saturating_sub(1);
    if total > 0 && reachable > 0 {
        reachable as f64 / total as f64
    } else {
        0.0
    }
}

/// Sums the per-source dependencies into one score per node.
fn accumulate_betweenness(nodes: &[NodeId], contributions: Vec<Vec<f64>>) -> HashMap<NodeId, f64> {
    let mut betweenness: HashMap<NodeId, f64> = nodes.iter().map(|&node| (node, 0.0)).collect();
    for dependencies in contributions {
        for (index, dependency) in dependencies.into_iter().enumerate() {
            *betweenness
                .get_mut(&nodes[index])
                .expect("every node has an entry") += dependency;
        }
    }
    betweenness
}

/// Scales the raw scores so that the maximum possible value is 1.
fn normalize_betweenness(
    mut betweenness: HashMap<NodeId, f64>,
    scale: f64,
) -> HashMap<NodeId, f64> {
    for value in betweenness.values_mut() {
        *value *= scale;
    }
    betweenness
}

/// One source's pass of Brandes' algorithm.
///
/// Returns, per node position, how much of the shortest-path traffic out of `source` runs
/// through that node. Summing these over every source gives betweenness. Kept as a free
/// function so that the serial and parallel drivers share one implementation.
fn brandes_from<G>(
    graph: &G,
    source: NodeId,
    nodes: &[NodeId],
    index_of: &HashMap<NodeId, usize>,
) -> Result<Vec<f64>, CentralityError>
where
    G: GraphView + ?Sized,
{
    let count = nodes.len();
    let source_index = index_of[&source];

    let mut order = Vec::with_capacity(count);
    let mut predecessors: Vec<Vec<usize>> = vec![Vec::new(); count];
    let mut path_count = vec![0.0_f64; count];
    let mut distance: Vec<Option<usize>> = vec![None; count];

    path_count[source_index] = 1.0;
    distance[source_index] = Some(0);

    let mut queue = VecDeque::from([source]);
    while let Some(node) = queue.pop_front() {
        let node_index = index_of[&node];
        order.push(node_index);

        for link in graph.links(node) {
            let target_index = *index_of
                .get(&link.target)
                .ok_or(CentralityError::InvalidNode(link.target))?;
            let step = distance[node_index].expect("dequeued nodes carry a distance") + 1;

            if distance[target_index].is_none() {
                distance[target_index] = Some(step);
                queue.push_back(link.target);
            }
            // A parallel edge is a second shortest path through the same pair, so it
            // contributes to the count exactly as a distinct route would.
            if distance[target_index] == Some(step) {
                path_count[target_index] += path_count[node_index];
                predecessors[target_index].push(node_index);
            }
        }
    }

    let mut dependency = vec![0.0_f64; count];
    for &node_index in order.iter().rev() {
        for &predecessor in &predecessors[node_index] {
            dependency[predecessor] +=
                (path_count[predecessor] / path_count[node_index]) * (1.0 + dependency[node_index]);
        }
    }
    // The source itself is not an intermediate node on its own paths.
    dependency[source_index] = 0.0;
    Ok(dependency)
}

/// Rejects a parameter that should be a probability.
fn check_probability(name: &'static str, value: f64) -> Result<(), CentralityError> {
    if !(0.0..=1.0).contains(&value) || value.is_nan() {
        return Err(CentralityError::InvalidParameter {
            name,
            value,
            expected: "a value between 0 and 1",
        });
    }
    Ok(())
}

/// Rejects a parameter that should not be negative.
fn check_non_negative(name: &'static str, value: f64) -> Result<(), CentralityError> {
    if value < 0.0 || value.is_nan() {
        return Err(CentralityError::InvalidParameter {
            name,
            value,
            expected: "a non-negative value",
        });
    }
    Ok(())
}

/// Scales a score map so its values sum to 1, leaving an all-zero map untouched.
fn normalize_to_sum(scores: &mut HashMap<NodeId, f64>) {
    let total: f64 = scores.values().sum();
    if total > 0.0 {
        for value in scores.values_mut() {
            *value /= total;
        }
    }
}

#[cfg(test)]
#[cfg(feature = "graph")]
mod tests {
    use super::*;
    use crate::core::EdgeFiltered;
    use crate::graph::Graph;
    use float_cmp::approx_eq;

    /// A path a - b - c, in which only `b` sits between anything.
    fn path_of_three() -> (Graph<u32, u32, ()>, NodeId, NodeId, NodeId) {
        let mut graph = Graph::<u32, u32, ()>::new(false);
        let a = graph.add_node(0);
        let b = graph.add_node(1);
        let c = graph.add_node(2);
        graph.add_edge(a, b, 1, ()).unwrap();
        graph.add_edge(b, c, 1, ()).unwrap();
        (graph, a, b, c)
    }

    #[test]
    fn test_degree_counts_parallel_edges_separately() {
        let mut graph = Graph::<u32, u32, ()>::new(false);
        let a = graph.add_node(0);
        let b = graph.add_node(1);
        graph.add_edge(a, b, 1, ()).unwrap();
        graph.add_edge(a, b, 2, ()).unwrap();

        let degree = graph.degree_centrality();
        assert_eq!(degree[&a], 2);
        assert_eq!(degree[&b], 2);
    }

    #[test]
    fn test_betweenness_finds_the_middle_of_a_path() {
        let (graph, a, b, c) = path_of_three();
        let betweenness = graph.betweenness_centrality().unwrap();

        assert!(approx_eq!(f64, betweenness[&b], 1.0, epsilon = 1e-9));
        assert!(approx_eq!(f64, betweenness[&a], 0.0, epsilon = 1e-9));
        assert!(approx_eq!(f64, betweenness[&c], 0.0, epsilon = 1e-9));
    }

    /// The normalisation constant divides by `(n - 1)(n - 2)`, so graphs of one or two
    /// nodes have to be answered before it is applied rather than yielding NaN.
    #[test]
    fn test_betweenness_of_tiny_graphs_is_zero() {
        let mut graph = Graph::<u32, u32, ()>::new(false);
        let a = graph.add_node(0);
        let b = graph.add_node(1);
        graph.add_edge(a, b, 1, ()).unwrap();

        let betweenness = graph.betweenness_centrality().unwrap();
        assert_eq!(betweenness[&a], 0.0);
        assert_eq!(betweenness[&b], 0.0);
        assert!(betweenness.values().all(|value| value.is_finite()));
    }

    #[test]
    fn test_betweenness_survives_sparse_ids() {
        let mut graph = Graph::<u32, u32, ()>::new(false);
        let a = graph.add_node(0);
        let gap = graph.add_node(1);
        let b = graph.add_node(2);
        let c = graph.add_node(3);
        graph.remove_node(gap).unwrap();
        graph.add_edge(a, b, 1, ()).unwrap();
        graph.add_edge(b, c, 1, ()).unwrap();

        let betweenness = graph.betweenness_centrality().unwrap();
        assert_eq!(betweenness.len(), 3);
        assert!(approx_eq!(f64, betweenness[&b], 1.0, epsilon = 1e-9));
    }

    #[test]
    fn test_closeness_prefers_the_centre() {
        let (graph, a, b, _) = path_of_three();
        let closeness = graph.closeness_centrality().unwrap();

        // `b` reaches both others in one hop; `a` needs one and two.
        assert!(approx_eq!(f64, closeness[&b], 1.0, epsilon = 1e-9));
        assert!(approx_eq!(f64, closeness[&a], 2.0 / 3.0, epsilon = 1e-9));
    }

    #[test]
    fn test_isolated_node_has_zero_closeness() {
        let mut graph = Graph::<u32, u32, ()>::new(false);
        let alone = graph.add_node(0);
        assert_eq!(graph.closeness_centrality().unwrap()[&alone], 0.0);
    }

    #[test]
    fn test_pagerank_sums_to_one_and_favours_the_hub() {
        let mut graph = Graph::<f64, u32, ()>::new(true);
        let hub = graph.add_node(0);
        let spokes: Vec<_> = (1..4).map(|i| graph.add_node(i)).collect();
        for &spoke in &spokes {
            graph.add_edge(spoke, hub, 1.0, ()).unwrap();
        }

        let rank = graph.pagerank(0.85, 100, 1e-9).unwrap();
        let total: f64 = rank.values().sum();
        assert!(approx_eq!(f64, total, 1.0, epsilon = 1e-9));
        assert!(spokes.iter().all(|spoke| rank[&hub] > rank[spoke]));
    }

    #[test]
    fn test_pagerank_rejects_a_damping_factor_outside_zero_to_one() {
        let mut graph = Graph::<f64, u32, ()>::new(true);
        graph.add_node(0);

        assert!(matches!(
            graph.pagerank(1.5, 10, 1e-6),
            Err(CentralityError::InvalidParameter {
                name: "damping",
                ..
            })
        ));
        assert!(matches!(
            graph.pagerank(f64::NAN, 10, 1e-6),
            Err(CentralityError::InvalidParameter {
                name: "damping",
                ..
            })
        ));
    }

    #[test]
    fn test_katz_rejects_a_non_positive_attenuation() {
        let mut graph = Graph::<f64, u32, ()>::new(false);
        graph.add_node(0);

        assert!(matches!(
            graph.katz_centrality(0.0, 1.0, 10, 1e-6),
            Err(CentralityError::InvalidParameter { name: "alpha", .. })
        ));
    }

    #[test]
    fn test_eigenvector_and_katz_rank_the_hub_highest() {
        let mut graph = Graph::<f64, u32, ()>::new(false);
        let hub = graph.add_node(0);
        let spokes: Vec<_> = (1..5).map(|i| graph.add_node(i)).collect();
        for &spoke in &spokes {
            graph.add_edge(hub, spoke, 1.0, ()).unwrap();
        }

        let eigenvector = graph.eigenvector_centrality(200, 1e-9).unwrap();
        let katz = graph.katz_centrality(0.1, 1.0, 200, 1e-9).unwrap();
        for spoke in &spokes {
            assert!(eigenvector[&hub] > eigenvector[spoke]);
            assert!(katz[&hub] > katz[spoke]);
        }
    }

    #[test]
    fn test_centrality_over_a_filtered_view() {
        let mut graph = Graph::<u32, u32, &str>::new(false);
        let a = graph.add_node(0);
        let b = graph.add_node(1);
        let c = graph.add_node(2);
        let rail = graph.add_edge(a, b, 1, "rail").unwrap();
        graph.add_edge(b, c, 1, "road").unwrap();

        let rail_only = EdgeFiltered::new(&graph, |edge| edge == rail);
        let degree = rail_only.degree_centrality();
        assert_eq!(degree[&a], 1);
        assert_eq!(degree[&c], 0);
    }

    /// The point of the parallel variants is that they change only the wall time.
    #[cfg(feature = "parallel")]
    #[test]
    fn test_the_parallel_measures_agree_with_the_serial_ones() {
        let mut graph = Graph::<u32, u32, ()>::new(false);
        let nodes: Vec<_> = (0..40).map(|i| graph.add_node(i)).collect();
        // A ring with chords, so the shortest-path structure is not trivial.
        for i in 0..nodes.len() {
            graph
                .add_edge(nodes[i], nodes[(i + 1) % nodes.len()], 1, ())
                .unwrap();
            if i % 3 == 0 {
                graph
                    .add_edge(nodes[i], nodes[(i + 7) % nodes.len()], 1, ())
                    .unwrap();
            }
        }

        let serial = graph.betweenness_centrality().unwrap();
        let parallel = graph.par_betweenness_centrality().unwrap();
        assert_eq!(serial.len(), parallel.len());
        for (node, value) in &serial {
            assert!(
                approx_eq!(f64, *value, parallel[node], epsilon = 1e-12),
                "{:?}: {} vs {}",
                node,
                value,
                parallel[node]
            );
        }

        assert_eq!(
            graph.closeness_centrality().unwrap(),
            graph.par_closeness_centrality().unwrap()
        );
    }
}