ratio_bus/
lib.rs

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

ratio_bus/
lib.rs
