rust_igraph/algorithms/properties/
joint_degree_matrix.rs

1//! Joint degree matrix (ALGO-PR-050).
2//!
3//! Counterpart of `igraph_joint_degree_matrix()` from
4//! `references/igraph/src/misc/mixing.c:528`.
5//!
6//! Entry `(i-1, j-1)` of the matrix counts edges (or sums edge weights)
7//! between degree-`i` and degree-`j` vertices.
8//! For undirected graphs the matrix is symmetric.
9//! For directed graphs rows correspond to out-degree, columns to in-degree.
10
11use crate::core::{Graph, IgraphError, IgraphResult};
12
13/// Compute the joint degree matrix of a graph.
14///
15/// Returns a matrix where entry `(d1-1, d2-1)` counts (or sums weights of)
16/// edges between vertices of degree `d1` and degree `d2`.
17///
18/// - For undirected graphs: both endpoints contribute their total degree.
19///   The matrix is symmetric. Each edge is counted once (but contributes
20///   to both `[d_u-1][d_v-1]` and `[d_v-1][d_u-1]` unless `d_u == d_v`).
21/// - For directed graphs: row index is source out-degree, column index is
22///   target in-degree. Each edge is counted exactly once.
23///
24/// `weights`: optional edge weights (length must equal `ecount()`).
25///
26/// Returns an empty matrix for graphs with no edges.
27///
28/// # Examples
29///
30/// ```
31/// use rust_igraph::{Graph, joint_degree_matrix};
32///
33/// // Triangle: all vertices have degree 2. JDM[1][1] = 3 (three edges).
34/// let mut g = Graph::with_vertices(3);
35/// g.add_edge(0, 1).unwrap();
36/// g.add_edge(1, 2).unwrap();
37/// g.add_edge(2, 0).unwrap();
38/// let jdm = joint_degree_matrix(&g, None).unwrap();
39/// // max_degree = 2, so matrix is 2×2. Only [1][1] is nonzero.
40/// assert_eq!(jdm.len(), 2);
41/// assert!((jdm[1][1] - 3.0).abs() < 1e-10);
42/// assert!((jdm[0][0]).abs() < 1e-10);
43/// ```
44pub fn joint_degree_matrix(graph: &Graph, weights: Option<&[f64]>) -> IgraphResult<Vec<Vec<f64>>> {
45    let n = graph.vcount();
46    let ecount = graph.ecount();
47
48    if let Some(w) = weights {
49        if w.len() != ecount {
50            return Err(IgraphError::InvalidArgument(format!(
51                "joint_degree_matrix: weights length ({}) does not match edge count ({ecount})",
52                w.len()
53            )));
54        }
55    }
56
57    if ecount == 0 {
58        return Ok(Vec::new());
59    }
60
61    let n_usize = n as usize;
62    let directed = graph.is_directed();
63
64    let m_u32 =
65        u32::try_from(ecount).map_err(|_| IgraphError::Internal("ecount exceeds u32::MAX"))?;
66
67    if directed {
68        // Compute out-degree and in-degree for each vertex.
69        let mut out_deg = vec![0u32; n_usize];
70        let mut in_deg = vec![0u32; n_usize];
71
72        for eid in 0..m_u32 {
73            let (from, to) = graph.edge(eid)?;
74            out_deg[from as usize] += 1;
75            in_deg[to as usize] += 1;
76        }
77
78        let max_out = out_deg.iter().copied().max().unwrap_or(0) as usize;
79        let max_in = in_deg.iter().copied().max().unwrap_or(0) as usize;
80
81        if max_out == 0 || max_in == 0 {
82            return Ok(Vec::new());
83        }
84
85        let mut jdm = vec![vec![0.0_f64; max_in]; max_out];
86
87        for eid in 0..m_u32 {
88            let (from, to) = graph.edge(eid)?;
89            let d_out = out_deg[from as usize] as usize;
90            let d_in = in_deg[to as usize] as usize;
91            let w = weights.map_or(1.0, |ws| ws[eid as usize]);
92            jdm[d_out - 1][d_in - 1] += w;
93        }
94
95        Ok(jdm)
96    } else {
97        // Compute degree for each vertex.
98        let mut deg = vec![0u32; n_usize];
99
100        for eid in 0..m_u32 {
101            let (from, to) = graph.edge(eid)?;
102            deg[from as usize] += 1;
103            deg[to as usize] += 1;
104        }
105
106        let max_deg = deg.iter().copied().max().unwrap_or(0) as usize;
107
108        if max_deg == 0 {
109            return Ok(Vec::new());
110        }
111
112        let mut jdm = vec![vec![0.0_f64; max_deg]; max_deg];
113
114        for eid in 0..m_u32 {
115            let (from, to) = graph.edge(eid)?;
116            let d1 = deg[from as usize] as usize;
117            let d2 = deg[to as usize] as usize;
118            let w = weights.map_or(1.0, |ws| ws[eid as usize]);
119            jdm[d1 - 1][d2 - 1] += w;
120            if d1 != d2 {
121                jdm[d2 - 1][d1 - 1] += w;
122            }
123        }
124
125        Ok(jdm)
126    }
127}
128
129#[cfg(test)]
130mod tests {
131    use super::*;
132
133    fn close(a: f64, b: f64) -> bool {
134        (a - b).abs() < 1e-10
135    }
136
137    #[test]
138    fn empty_graph() {
139        let g = Graph::with_vertices(0);
140        let jdm = joint_degree_matrix(&g, None).unwrap();
141        assert!(jdm.is_empty());
142    }
143
144    #[test]
145    fn no_edges() {
146        let g = Graph::with_vertices(5);
147        let jdm = joint_degree_matrix(&g, None).unwrap();
148        assert!(jdm.is_empty());
149    }
150
151    #[test]
152    fn triangle_undirected() {
153        // All degree-2. JDM is 2×2, only [1][1] = 3 (three edges between d=2 vertices).
154        let mut g = Graph::with_vertices(3);
155        g.add_edge(0, 1).unwrap();
156        g.add_edge(1, 2).unwrap();
157        g.add_edge(2, 0).unwrap();
158        let jdm = joint_degree_matrix(&g, None).unwrap();
159        assert_eq!(jdm.len(), 2);
160        assert_eq!(jdm[0].len(), 2);
161        assert!(close(jdm[0][0], 0.0));
162        assert!(close(jdm[0][1], 0.0));
163        assert!(close(jdm[1][0], 0.0));
164        assert!(close(jdm[1][1], 3.0));
165    }
166
167    #[test]
168    fn path_3_undirected() {
169        // 0-1-2: deg = [1, 2, 1].
170        // Edge (0,1): d=1 and d=2 → jdm[0][1] += 1, jdm[1][0] += 1.
171        // Edge (1,2): d=2 and d=1 → jdm[1][0] += 1, jdm[0][1] += 1.
172        // Result: jdm[0][1] = 2, jdm[1][0] = 2. Max degree = 2.
173        let mut g = Graph::with_vertices(3);
174        g.add_edge(0, 1).unwrap();
175        g.add_edge(1, 2).unwrap();
176        let jdm = joint_degree_matrix(&g, None).unwrap();
177        assert_eq!(jdm.len(), 2);
178        assert!(close(jdm[0][0], 0.0));
179        assert!(close(jdm[0][1], 2.0));
180        assert!(close(jdm[1][0], 2.0));
181        assert!(close(jdm[1][1], 0.0));
182    }
183
184    #[test]
185    fn star_undirected() {
186        // Star: 0 connected to 1,2,3. deg(0)=3, deg(1)=deg(2)=deg(3)=1.
187        // Edges: (0,1), (0,2), (0,3). Each connects d=3 with d=1.
188        // jdm[2][0] += 3, jdm[0][2] += 3. Matrix is 3×3.
189        let mut g = Graph::with_vertices(4);
190        g.add_edge(0, 1).unwrap();
191        g.add_edge(0, 2).unwrap();
192        g.add_edge(0, 3).unwrap();
193        let jdm = joint_degree_matrix(&g, None).unwrap();
194        assert_eq!(jdm.len(), 3);
195        assert!(close(jdm[0][2], 3.0));
196        assert!(close(jdm[2][0], 3.0));
197        assert!(close(jdm[0][0], 0.0));
198        assert!(close(jdm[1][1], 0.0));
199        assert!(close(jdm[2][2], 0.0));
200    }
201
202    #[test]
203    fn directed_simple() {
204        // 0→1→2. out_deg = [1,1,0], in_deg = [0,1,1].
205        // Edge (0→1): out_deg[0]=1, in_deg[1]=1 → jdm[0][0] += 1.
206        // Edge (1→2): out_deg[1]=1, in_deg[2]=1 → jdm[0][0] += 1.
207        // max_out = 1, max_in = 1. Matrix is 1×1: jdm[0][0] = 2.
208        let mut g = Graph::new(3, true).unwrap();
209        g.add_edge(0, 1).unwrap();
210        g.add_edge(1, 2).unwrap();
211        let jdm = joint_degree_matrix(&g, None).unwrap();
212        assert_eq!(jdm.len(), 1);
213        assert_eq!(jdm[0].len(), 1);
214        assert!(close(jdm[0][0], 2.0));
215    }
216
217    #[test]
218    fn directed_star() {
219        // 0→1, 0→2, 0→3. out_deg = [3,0,0,0], in_deg = [0,1,1,1].
220        // Edge (0→1): out_deg[0]=3, in_deg[1]=1 → jdm[2][0] += 1.
221        // Edge (0→2): out_deg[0]=3, in_deg[2]=1 → jdm[2][0] += 1.
222        // Edge (0→3): out_deg[0]=3, in_deg[3]=1 → jdm[2][0] += 1.
223        // max_out = 3, max_in = 1. Matrix is 3×1: jdm[2][0] = 3.
224        let mut g = Graph::new(4, true).unwrap();
225        g.add_edge(0, 1).unwrap();
226        g.add_edge(0, 2).unwrap();
227        g.add_edge(0, 3).unwrap();
228        let jdm = joint_degree_matrix(&g, None).unwrap();
229        assert_eq!(jdm.len(), 3);
230        assert_eq!(jdm[0].len(), 1);
231        assert!(close(jdm[2][0], 3.0));
232        assert!(close(jdm[0][0], 0.0));
233        assert!(close(jdm[1][0], 0.0));
234    }
235
236    #[test]
237    fn weighted_undirected() {
238        // 0-1 (w=3), 1-2 (w=2). deg = [1, 2, 1].
239        // Edge 0-1: d1=1, d2=2 → jdm[0][1] += 3, jdm[1][0] += 3.
240        // Edge 1-2: d1=2, d2=1 → jdm[1][0] += 2, jdm[0][1] += 2.
241        // jdm[0][1] = 5, jdm[1][0] = 5.
242        let mut g = Graph::with_vertices(3);
243        g.add_edge(0, 1).unwrap();
244        g.add_edge(1, 2).unwrap();
245        let jdm = joint_degree_matrix(&g, Some(&[3.0, 2.0])).unwrap();
246        assert!(close(jdm[0][1], 5.0));
247        assert!(close(jdm[1][0], 5.0));
248        assert!(close(jdm[0][0], 0.0));
249        assert!(close(jdm[1][1], 0.0));
250    }
251
252    #[test]
253    fn self_loop_undirected() {
254        // 0-0 (self-loop), 0-1. deg = [3, 1] (self-loop adds 2 to degree of 0).
255        // Edge (0,0): d1=3, d2=3 (same vertex) → jdm[2][2] += 1. d1==d2 so no double.
256        // Edge (0,1): d1=3, d2=1 → jdm[2][0] += 1, jdm[0][2] += 1.
257        let mut g = Graph::with_vertices(2);
258        g.add_edge(0, 0).unwrap();
259        g.add_edge(0, 1).unwrap();
260        let jdm = joint_degree_matrix(&g, None).unwrap();
261        assert_eq!(jdm.len(), 3);
262        assert!(close(jdm[2][2], 1.0));
263        assert!(close(jdm[2][0], 1.0));
264        assert!(close(jdm[0][2], 1.0));
265    }
266
267    #[test]
268    fn total_equals_edge_count_undirected() {
269        // For unweighted undirected: sum of JDM diagonal = edges within same-degree,
270        // sum of upper triangle = edges between different-degree.
271        // Total sum should equal ecount (each edge counted once in lower+diagonal,
272        // but symmetric so total = 2*off_diag + diag... actually let's just
273        // check sum of upper-triangle + diagonal = ecount.
274        let mut g = Graph::with_vertices(5);
275        g.add_edge(0, 1).unwrap();
276        g.add_edge(1, 2).unwrap();
277        g.add_edge(2, 3).unwrap();
278        g.add_edge(3, 4).unwrap();
279        g.add_edge(4, 0).unwrap(); // cycle
280        let jdm = joint_degree_matrix(&g, None).unwrap();
281        // All vertices have degree 2. So jdm[1][1] = 5.
282        assert!(close(jdm[1][1], 5.0));
283        // Total of matrix (symmetric, same-degree so just diagonal): 5.
284        let total: f64 = jdm.iter().flat_map(|r| r.iter()).sum();
285        // For same-degree case, total = ecount (no double counting).
286        assert!(close(total, 5.0));
287    }
288
289    #[test]
290    fn weights_mismatch_errors() {
291        let mut g = Graph::with_vertices(2);
292        g.add_edge(0, 1).unwrap();
293        assert!(joint_degree_matrix(&g, Some(&[1.0, 2.0])).is_err());
294    }
295
296    #[test]
297    fn k4_undirected() {
298        // K4: all degree 3. 6 edges. JDM[2][2] = 6.
299        let mut g = Graph::with_vertices(4);
300        for u in 0..4u32 {
301            for v in (u + 1)..4 {
302                g.add_edge(u, v).unwrap();
303            }
304        }
305        let jdm = joint_degree_matrix(&g, None).unwrap();
306        assert_eq!(jdm.len(), 3);
307        assert!(close(jdm[2][2], 6.0));
308        assert!(close(jdm[0][0], 0.0));
309    }
310}
rust_igraph/algorithms/properties/joint_degree_matrix.rs

rust_igraph/algorithms/properties/
joint_degree_matrix.rs
