ratio_bus/
lib.rs

1#![doc = include_str!("../README.md")]
2
3// ## License
4//
5// This Source Code Form is subject to the terms of the Mozilla Public License, v. 2.0.
6// If a copy of the MPL was not distributed with this file,
7// You can obtain one at <https://mozilla.org/MPL/2.0/>.
8//
9// **Code examples both in the docstrings and rendered documentation are free to use.**
10
11use nalgebra::{DMatrix, DVector};
12use snafu::prelude::*;
13
14#[derive(Debug, Snafu)]
15pub enum Error {
16    #[snafu(display("Gamma should positive, but it was '{gamma}'"))]
17    GammaTooLow { gamma: f64 },
18}
19
20/// Gamma bus detection.
21///
22/// # Arguments
23///
24/// * `adj`: Adjacency matrix corresponding to a graph.
25/// * `gamma`: Threshold factor for bus node degrees with respect to the median of the
26///   (nonzero) remaining candidate node degrees.
27///
28/// # Returns
29///
30/// A pair of vectors containing the bus node indices and remaining non-bus node
31/// indices, respectively.
32///
33/// # Errors
34///
35/// This function will return an error if gamma is too low.
36pub fn gamma(adj: &DMatrix<f64>, gamma: f64) -> Result<(Vec<usize>, Vec<usize>), Error> {
37    ensure!(gamma >= 0.0, GammaTooLowSnafu { gamma });
38
39    let (rows, cols) = adj.shape();
40    let dim = std::cmp::max(rows, cols);
41
42    if dim <= 1 {
43        return Ok((vec![], (0..dim).collect())); // Early return for trivial cases 0 and 1 node.
44    }
45
46    let (in_degree, out_degree) = get_degrees(adj);
47    let degrees = (in_degree + out_degree).data.as_vec().clone(); // Get regular vec.
48
49    let mut nonbus: Vec<usize> = vec![];
50    let mut candidates: Vec<usize> = vec![];
51    let mut bus: Vec<usize> = vec![];
52
53    for (i, d) in degrees.iter().copied().enumerate().take(dim) {
54        if d > 0 {
55            candidates.push(i);
56        } else {
57            nonbus.push(i);
58        }
59    }
60
61    candidates.sort_by(|a, b| degrees[*a].cmp(&degrees[*b]));
62    let mut degrees: Vec<f64> = candidates.iter().map(|&idx| degrees[idx] as f64).collect();
63
64    loop {
65        let mut added = false;
66        let mid = degrees.len() / 2;
67        let median = degrees[mid];
68        let threshold = gamma * median;
69        while let (Some(degree), Some(idx)) = (degrees.pop(), candidates.pop()) {
70            if degree < threshold {
71                candidates.push(idx);
72                degrees.push(degree);
73                break;
74            }
75            bus.push(idx);
76            added = true;
77        }
78        if !added {
79            break;
80        }
81    }
82    nonbus.extend(candidates);
83    nonbus.sort(); // Return nonbus in original order.
84    Ok((bus, nonbus))
85}
86
87/// Get the in and out degrees of a graph represented by its adjacency matrix.
88///
89/// # Example
90/// ```
91/// use nalgebra::{DMatrix, DVector};
92/// use ratio_bus::get_degrees;
93/// let adj = DMatrix::from_row_slice(3, 3, &[0.0, 5.0, 10.0, 1.0, 0.0, 0.0, 3.0, 10.0, 0.0]);
94/// let (ins, outs) = get_degrees(&adj);
95/// assert_eq!(ins, DVector::from_vec(vec![2, 1, 2]));
96/// assert_eq!(outs, DVector::from_vec(vec![2, 2, 1]));
97/// ```
98pub fn get_degrees(adj: &DMatrix<f64>) -> (DVector<usize>, DVector<usize>) {
99    let (rows, cols) = adj.shape();
100    let mut interfaces = DMatrix::<usize>::zeros(rows, cols);
101    let mut in_degree = DVector::<usize>::zeros(rows);
102    let mut out_degree = DVector::<usize>::zeros(cols);
103    for i in 0..rows {
104        for j in 0..cols {
105            if adj[(i, j)] <= 0.0 {
106                continue;
107            }
108            interfaces[(i, j)] = 1;
109            in_degree[i] += 1;
110            out_degree[j] += 1;
111        }
112    }
113    (in_degree, out_degree)
114}
115
116#[cfg(test)]
117mod tests {
118    #[allow(unused_imports)]
119    use pretty_assertions::{assert_eq, assert_ne, assert_str_eq};
120
121    use super::*;
122
123    #[test]
124    fn example() {
125        let adj = DMatrix::from_row_slice(
126            4,
127            4,
128            &[
129                0.0, 5.0, 10.0, 0.0, 1.0, 0.0, 0.0, 0.0, 3.0, 10.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
130            ],
131        );
132        let (ins, outs) = get_degrees(&adj);
133        assert_eq!(
134            ins,
135            DVector::from_vec(vec![2, 1, 2, 0]),
136            "In degrees should match."
137        );
138        assert_eq!(
139            outs,
140            DVector::from_vec(vec![2, 2, 1, 0]),
141            "Out degrees should match."
142        );
143
144        let (bus, nonbus) = gamma(&adj, 100.0).unwrap();
145        assert_eq!(bus, vec![], "Bus should be empty for gamma=100.");
146        assert_eq!(
147            nonbus,
148            vec![0, 1, 2, 3],
149            "Nonbus should be full for gamma=100."
150        );
151
152        let (bus, nonbus) = gamma(&adj, 1.2).unwrap();
153        assert_eq!(bus, vec![0], "Bus should contain 0 for gamma=1.2.");
154        assert_eq!(
155            nonbus,
156            vec![1, 2, 3],
157            "Nonbus should be [1,2] for gamma=100."
158        );
159    }
160
161    #[test]
162    fn edge_cases() {
163        let single = DMatrix::from_row_slice(1, 1, &[1.0]);
164        let result = gamma(&single, -1.0);
165        result.unwrap_err();
166        assert_eq!(gamma(&single, 1.0).unwrap(), (vec![], vec![0]));
167    }
168}
ratio_bus/lib.rs

ratio_bus/
lib.rs
