solvr 0.2.0

Advanced computing library for real-world problem solving - optimization, differential equations, interpolation, statistics, and more
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109

//! Betweenness centrality via Brandes' algorithm.
//!
//! Sequential BFS from each source, accumulates pair dependencies.
use crate::DType;

use numr::error::Result;
use numr::runtime::{Runtime, RuntimeClient};
use numr::sparse::SparseTensor;
use numr::tensor::Tensor;

use crate::graph::traits::types::GraphData;

use super::helpers::extract_csr_arrays;

/// Betweenness centrality via Brandes' algorithm.
///
/// Time: O(V * E) for unweighted graphs.
pub fn betweenness_centrality_impl<R, C>(
    _client: &C,
    graph: &GraphData<R>,
    normalized: bool,
) -> Result<Tensor<R>>
where
    R: Runtime<DType = DType>,
    C: RuntimeClient<R>,
{
    let (row_ptrs, col_indices, _values, n) = extract_csr_arrays(graph)?;

    let device = match &graph.adjacency {
        SparseTensor::Csr(csr) => csr.values().device().clone(),
        _ => unreachable!(),
    };

    if n == 0 {
        return Ok(Tensor::<R>::from_slice(&[] as &[f64], &[0], &device));
    }

    let mut betweenness = vec![0.0f64; n];

    // Brandes' algorithm: BFS from each source
    for s in 0..n {
        let mut stack = Vec::with_capacity(n);
        let mut predecessors: Vec<Vec<usize>> = vec![Vec::new(); n];
        let mut sigma = vec![0.0f64; n]; // number of shortest paths
        let mut dist = vec![-1i64; n];
        let mut delta = vec![0.0f64; n];

        sigma[s] = 1.0;
        dist[s] = 0;

        // BFS
        let mut queue = std::collections::VecDeque::new();
        queue.push_back(s);

        while let Some(v) = queue.pop_front() {
            stack.push(v);

            let start = row_ptrs[v] as usize;
            let end = row_ptrs[v + 1] as usize;
            for &w_idx in col_indices.iter().take(end).skip(start) {
                let w = w_idx as usize;

                // First visit?
                if dist[w] < 0 {
                    dist[w] = dist[v] + 1;
                    queue.push_back(w);
                }

                // Shortest path to w through v?
                if dist[w] == dist[v] + 1 {
                    sigma[w] += sigma[v];
                    predecessors[w].push(v);
                }
            }
        }

        // Back-propagation of dependencies
        while let Some(w) = stack.pop() {
            for &v in &predecessors[w] {
                delta[v] += (sigma[v] / sigma[w]) * (1.0 + delta[w]);
            }
            if w != s {
                betweenness[w] += delta[w];
            }
            delta[w] = 0.0;
        }
    }

    // For undirected graphs, each pair is counted twice
    if !graph.directed {
        for b in betweenness.iter_mut() {
            *b /= 2.0;
        }
    }

    // Normalize
    if normalized && n > 2 {
        let norm = if graph.directed {
            ((n - 1) * (n - 2)) as f64
        } else {
            ((n - 1) * (n - 2)) as f64 / 2.0
        };
        for b in betweenness.iter_mut() {
            *b /= norm;
        }
    }

    Ok(Tensor::<R>::from_slice(&betweenness, &[n], &device))
}