rust_igraph/algorithms/properties/
stochastic.rs

1//! Stochastic adjacency matrix (ALGO-PR-043).
2//!
3//! Returns the row-wise or column-wise normalized adjacency matrix,
4//! representing transition probabilities of a random walk.
5//!
6//! Counterpart of `igraph_get_stochastic` from
7//! `references/igraph/src/misc/conversion.c`.
8
9use crate::core::{Graph, IgraphError, IgraphResult};
10
11/// Compute the stochastic adjacency matrix of a graph.
12///
13/// The stochastic matrix is the adjacency matrix normalized so that each
14/// row (or column) sums to 1. Row-normalized = right-stochastic =
15/// random walk following edge directions. Column-normalized =
16/// left-stochastic = random walk against edge directions.
17///
18/// When a vertex has zero out-degree (or in-degree for column-wise),
19/// the corresponding row (or column) will be all zeros.
20///
21/// - `column_wise`: if false, row-normalize (right-stochastic);
22///   if true, column-normalize (left-stochastic).
23/// - `weights`: optional edge weights. If `None`, each edge counts as 1.
24///
25/// Returns an n×n matrix as `Vec<Vec<f64>>`.
26///
27/// # Examples
28///
29/// ```
30/// use rust_igraph::{Graph, get_stochastic};
31///
32/// // Path 0-1-2: row-stochastic
33/// let mut g = Graph::with_vertices(3);
34/// g.add_edge(0, 1).unwrap();
35/// g.add_edge(1, 2).unwrap();
36/// let s = get_stochastic(&g, false, None).unwrap();
37/// // Row 0 sums to 1: only edge to vertex 1
38/// assert!((s[0][1] - 1.0).abs() < 1e-10);
39/// // Row 1: edges to 0 and 2, each 0.5
40/// assert!((s[1][0] - 0.5).abs() < 1e-10);
41/// assert!((s[1][2] - 0.5).abs() < 1e-10);
42/// ```
43pub fn get_stochastic(
44    graph: &Graph,
45    column_wise: bool,
46    weights: Option<&[f64]>,
47) -> IgraphResult<Vec<Vec<f64>>> {
48    let n = graph.vcount() as usize;
49    let ecount = graph.ecount();
50
51    if let Some(w) = weights {
52        if w.len() != ecount {
53            return Err(IgraphError::InvalidArgument(format!(
54                "get_stochastic: weight vector length ({}) does not match edge count ({ecount})",
55                w.len()
56            )));
57        }
58    }
59
60    let directed = graph.is_directed();
61
62    // Compute strength (weighted degree) for normalization
63    let mut sums = vec![0.0_f64; n];
64    for eid in 0..ecount {
65        #[allow(clippy::cast_possible_truncation)]
66        let (from, to) = graph.edge(eid as u32)?;
67        let w = weights.map_or(1.0, |ws| ws[eid]);
68
69        if directed {
70            if column_wise {
71                sums[to as usize] += w;
72            } else {
73                sums[from as usize] += w;
74            }
75        } else {
76            // Undirected: both endpoints get the weight
77            sums[from as usize] += w;
78            if from == to {
79                sums[from as usize] += w;
80            } else {
81                sums[to as usize] += w;
82            }
83        }
84    }
85
86    let mut res = vec![vec![0.0_f64; n]; n];
87
88    if directed {
89        for eid in 0..ecount {
90            #[allow(clippy::cast_possible_truncation)]
91            let (from, to) = graph.edge(eid as u32)?;
92            let f = from as usize;
93            let t = to as usize;
94            let w = weights.map_or(1.0, |ws| ws[eid]);
95
96            let divisor = if column_wise { sums[t] } else { sums[f] };
97            if divisor > 0.0 {
98                res[f][t] += w / divisor;
99            }
100        }
101    } else {
102        for eid in 0..ecount {
103            #[allow(clippy::cast_possible_truncation)]
104            let (from, to) = graph.edge(eid as u32)?;
105            let f = from as usize;
106            let t = to as usize;
107            let w = weights.map_or(1.0, |ws| ws[eid]);
108
109            let div_from = if column_wise { sums[t] } else { sums[f] };
110            if div_from > 0.0 {
111                res[f][t] += w / div_from;
112            }
113
114            let div_to = if column_wise { sums[f] } else { sums[t] };
115            if div_to > 0.0 {
116                res[t][f] += w / div_to;
117            }
118        }
119    }
120
121    Ok(res)
122}
123
124#[cfg(test)]
125mod tests {
126    use super::*;
127
128    fn approx_eq(a: f64, b: f64) -> bool {
129        (a - b).abs() < 1e-10
130    }
131
132    fn row_sum(matrix: &[Vec<f64>], row: usize) -> f64 {
133        matrix[row].iter().sum()
134    }
135
136    fn col_sum(matrix: &[Vec<f64>], col: usize) -> f64 {
137        matrix.iter().map(|row| row[col]).sum()
138    }
139
140    #[test]
141    fn empty_graph() {
142        let g = Graph::with_vertices(0);
143        let s = get_stochastic(&g, false, None).unwrap();
144        assert!(s.is_empty());
145    }
146
147    #[test]
148    fn isolated_vertex() {
149        let g = Graph::with_vertices(2);
150        let s = get_stochastic(&g, false, None).unwrap();
151        // All zeros
152        for row in &s {
153            for &val in row {
154                assert!(approx_eq(val, 0.0));
155            }
156        }
157    }
158
159    #[test]
160    fn single_edge_row() {
161        let mut g = Graph::with_vertices(2);
162        g.add_edge(0, 1).unwrap();
163        let s = get_stochastic(&g, false, None).unwrap();
164        // Row 0: only connects to 1, sum = 1
165        assert!(approx_eq(s[0][1], 1.0));
166        assert!(approx_eq(s[1][0], 1.0));
167        assert!(approx_eq(row_sum(&s, 0), 1.0));
168        assert!(approx_eq(row_sum(&s, 1), 1.0));
169    }
170
171    #[test]
172    fn path_3_row() {
173        let mut g = Graph::with_vertices(3);
174        g.add_edge(0, 1).unwrap();
175        g.add_edge(1, 2).unwrap();
176        let s = get_stochastic(&g, false, None).unwrap();
177        // vertex 0: deg=1, only to 1
178        assert!(approx_eq(s[0][1], 1.0));
179        // vertex 1: deg=2, to 0 and 2
180        assert!(approx_eq(s[1][0], 0.5));
181        assert!(approx_eq(s[1][2], 0.5));
182        // vertex 2: deg=1, only to 1
183        assert!(approx_eq(s[2][1], 1.0));
184        // All rows sum to 1
185        for i in 0..3 {
186            assert!(approx_eq(row_sum(&s, i), 1.0));
187        }
188    }
189
190    #[test]
191    fn path_3_column() {
192        let mut g = Graph::with_vertices(3);
193        g.add_edge(0, 1).unwrap();
194        g.add_edge(1, 2).unwrap();
195        let s = get_stochastic(&g, true, None).unwrap();
196        // Columns should sum to 1 for non-isolated vertices
197        assert!(approx_eq(col_sum(&s, 0), 1.0));
198        assert!(approx_eq(col_sum(&s, 1), 1.0));
199        assert!(approx_eq(col_sum(&s, 2), 1.0));
200    }
201
202    #[test]
203    fn directed_row() {
204        let mut g = Graph::new(3, true).unwrap();
205        g.add_edge(0, 1).unwrap();
206        g.add_edge(0, 2).unwrap();
207        g.add_edge(1, 2).unwrap();
208        let s = get_stochastic(&g, false, None).unwrap();
209        // out-deg(0)=2: s[0][1]=0.5, s[0][2]=0.5
210        assert!(approx_eq(s[0][1], 0.5));
211        assert!(approx_eq(s[0][2], 0.5));
212        // out-deg(1)=1: s[1][2]=1.0
213        assert!(approx_eq(s[1][2], 1.0));
214        // out-deg(2)=0: row all zeros
215        assert!(approx_eq(row_sum(&s, 2), 0.0));
216    }
217
218    #[test]
219    fn directed_column() {
220        let mut g = Graph::new(3, true).unwrap();
221        g.add_edge(0, 1).unwrap();
222        g.add_edge(0, 2).unwrap();
223        g.add_edge(1, 2).unwrap();
224        let s = get_stochastic(&g, true, None).unwrap();
225        // in-deg(0)=0: column all zeros
226        assert!(approx_eq(col_sum(&s, 0), 0.0));
227        // in-deg(1)=1: s[0][1]=1.0
228        assert!(approx_eq(s[0][1], 1.0));
229        // in-deg(2)=2: s[0][2]+s[1][2]=1.0
230        assert!(approx_eq(s[0][2], 0.5));
231        assert!(approx_eq(s[1][2], 0.5));
232    }
233
234    #[test]
235    fn weighted_row() {
236        let mut g = Graph::with_vertices(3);
237        g.add_edge(0, 1).unwrap();
238        g.add_edge(0, 2).unwrap();
239        let weights = vec![2.0, 3.0];
240        let s = get_stochastic(&g, false, Some(&weights)).unwrap();
241        // strength(0)=5: s[0][1]=2/5, s[0][2]=3/5
242        assert!(approx_eq(s[0][1], 2.0 / 5.0));
243        assert!(approx_eq(s[0][2], 3.0 / 5.0));
244        assert!(approx_eq(row_sum(&s, 0), 1.0));
245    }
246
247    #[test]
248    fn self_loop() {
249        let mut g = Graph::with_vertices(2);
250        g.add_edge(0, 0).unwrap();
251        g.add_edge(0, 1).unwrap();
252        let s = get_stochastic(&g, false, None).unwrap();
253        // deg(0)=3 (self-loop counts 2 + edge to 1)
254        // s[0][0] = 1/3 (self-loop once in matrix) — wait, the C code
255        // adds to res[from][to], so self-loop adds w/sum once.
256        // Actually for undirected, both from→to and to→from are added.
257        // Self-loop: from=to=0, so res[0][0] gets w/sums[0] twice
258        // sums[0] = 2 + 1 = 3 (self-loop=2, edge=1)
259        // res[0][0] += 1/3 (first add) + 1/3 (second add) = 2/3
260        // res[0][1] += 1/3
261        // res[1][0] += 1/1 = 1.0
262        assert!(approx_eq(s[0][0], 2.0 / 3.0));
263        assert!(approx_eq(s[0][1], 1.0 / 3.0));
264        assert!(approx_eq(row_sum(&s, 0), 1.0));
265    }
266
267    #[test]
268    fn complete_k3_row() {
269        let mut g = Graph::with_vertices(3);
270        g.add_edge(0, 1).unwrap();
271        g.add_edge(0, 2).unwrap();
272        g.add_edge(1, 2).unwrap();
273        let s = get_stochastic(&g, false, None).unwrap();
274        // All degrees = 2, uniform distribution to neighbors
275        for i in 0..3 {
276            assert!(approx_eq(row_sum(&s, i), 1.0));
277        }
278        assert!(approx_eq(s[0][1], 0.5));
279        assert!(approx_eq(s[0][2], 0.5));
280    }
281
282    #[test]
283    fn weight_mismatch() {
284        let mut g = Graph::with_vertices(2);
285        g.add_edge(0, 1).unwrap();
286        let result = get_stochastic(&g, false, Some(&[1.0, 2.0]));
287        assert!(result.is_err());
288    }
289}
rust_igraph/algorithms/properties/stochastic.rs

rust_igraph/algorithms/properties/
stochastic.rs
