rust_igraph/algorithms/properties/
laplacian.rs

1//! Graph Laplacian matrix (ALGO-PR-041).
2//!
3//! Returns the Laplacian matrix L = D - A in dense form, with optional
4//! normalization (unnormalized, symmetric, left-stochastic, right-stochastic).
5//!
6//! Counterpart of `igraph_get_laplacian` from
7//! `references/igraph/src/properties/spectral.c`.
8
9use crate::algorithms::properties::degree::DegreeMode;
10use crate::core::{Graph, IgraphError, IgraphResult};
11
12/// Laplacian normalization mode.
13#[derive(Debug, Clone, Copy, PartialEq, Eq)]
14pub enum LaplacianNormalization {
15    /// L = D - A (unnormalized).
16    Unnormalized,
17    /// Symmetric normalization: D^{-1/2} L D^{-1/2}.
18    Symmetric,
19    /// Left stochastic: D^{-1} L (row-normalized).
20    Left,
21    /// Right stochastic: L D^{-1} (column-normalized).
22    Right,
23}
24
25/// Compute the Laplacian matrix of a graph.
26///
27/// Returns an n×n matrix stored in row-major order as `Vec<Vec<f64>>`.
28///
29/// - `mode`: degree mode (Out/In/All). For undirected graphs, forced to All.
30/// - `normalization`: normalization method.
31/// - `weights`: optional edge weights. If `None`, unweighted.
32///
33/// # Examples
34///
35/// ```
36/// use rust_igraph::{Graph, DegreeMode, get_laplacian, LaplacianNormalization};
37///
38/// // Path 0-1-2: unnormalized Laplacian
39/// let mut g = Graph::with_vertices(3);
40/// g.add_edge(0, 1).unwrap();
41/// g.add_edge(1, 2).unwrap();
42/// let lap = get_laplacian(&g, DegreeMode::All, LaplacianNormalization::Unnormalized, None).unwrap();
43/// assert!((lap[0][0] - 1.0).abs() < 1e-10);  // deg(0) = 1
44/// assert!((lap[1][1] - 2.0).abs() < 1e-10);  // deg(1) = 2
45/// assert!((lap[0][1] + 1.0).abs() < 1e-10);  // -A[0][1] = -1
46/// ```
47pub fn get_laplacian(
48    graph: &Graph,
49    mode: DegreeMode,
50    normalization: LaplacianNormalization,
51    weights: Option<&[f64]>,
52) -> IgraphResult<Vec<Vec<f64>>> {
53    let n = graph.vcount() as usize;
54    let ecount = graph.ecount();
55
56    if let Some(w) = weights {
57        if w.len() != ecount {
58            return Err(IgraphError::InvalidArgument(format!(
59                "get_laplacian: weight vector length ({}) does not match edge count ({ecount})",
60                w.len()
61            )));
62        }
63    }
64
65    let directed = graph.is_directed();
66    let mode_eff = if directed { mode } else { DegreeMode::All };
67    let treat_as_undirected = !directed || mode == DegreeMode::All;
68
69    // Compute degree/strength
70    let mut degree = vec![0.0_f64; n];
71    for eid in 0..ecount {
72        #[allow(clippy::cast_possible_truncation)]
73        let (from, to) = graph.edge(eid as u32)?;
74        let w = weights.map_or(1.0, |ws| ws[eid]);
75
76        match mode_eff {
77            DegreeMode::Out => {
78                degree[from as usize] += w;
79                if !directed && from != to {
80                    degree[to as usize] += w;
81                }
82            }
83            DegreeMode::In => {
84                degree[to as usize] += w;
85                if !directed && from != to {
86                    degree[from as usize] += w;
87                }
88            }
89            DegreeMode::All => {
90                degree[from as usize] += w;
91                if from == to {
92                    degree[from as usize] += w;
93                } else {
94                    degree[to as usize] += w;
95                }
96            }
97        }
98    }
99
100    let mut res = vec![vec![0.0_f64; n]; n];
101
102    // Set diagonal based on normalization
103    let mut norm_factor = vec![0.0_f64; n];
104    for i in 0..n {
105        match normalization {
106            LaplacianNormalization::Unnormalized => {
107                res[i][i] = degree[i];
108            }
109            LaplacianNormalization::Symmetric => {
110                if degree[i] > 0.0 {
111                    res[i][i] = 1.0;
112                    norm_factor[i] = 1.0 / degree[i].sqrt();
113                }
114            }
115            LaplacianNormalization::Left | LaplacianNormalization::Right => {
116                if degree[i] > 0.0 {
117                    res[i][i] = 1.0;
118                    norm_factor[i] = 1.0 / degree[i];
119                }
120            }
121        }
122    }
123
124    // Set off-diagonal
125    for eid in 0..ecount {
126        #[allow(clippy::cast_possible_truncation)]
127        let (from, to) = graph.edge(eid as u32)?;
128        let w = weights.map_or(1.0, |ws| ws[eid]);
129        let f = from as usize;
130        let t = to as usize;
131
132        match normalization {
133            LaplacianNormalization::Unnormalized => {
134                res[f][t] -= w;
135                if treat_as_undirected {
136                    res[t][f] -= w;
137                }
138            }
139            LaplacianNormalization::Symmetric => {
140                let nf = norm_factor[f] * norm_factor[t];
141                let wn = w * nf;
142                res[f][t] -= wn;
143                if treat_as_undirected {
144                    res[t][f] -= wn;
145                }
146            }
147            LaplacianNormalization::Left => {
148                res[f][t] -= w * norm_factor[f];
149                if treat_as_undirected {
150                    res[t][f] -= w * norm_factor[t];
151                }
152            }
153            LaplacianNormalization::Right => {
154                res[f][t] -= w * norm_factor[t];
155                if treat_as_undirected {
156                    res[t][f] -= w * norm_factor[f];
157                }
158            }
159        }
160    }
161
162    Ok(res)
163}
164
165#[cfg(test)]
166mod tests {
167    use super::*;
168
169    fn approx_eq(a: f64, b: f64) -> bool {
170        (a - b).abs() < 1e-10
171    }
172
173    #[test]
174    fn empty_graph() {
175        let g = Graph::with_vertices(0);
176        let lap = get_laplacian(
177            &g,
178            DegreeMode::All,
179            LaplacianNormalization::Unnormalized,
180            None,
181        )
182        .unwrap();
183        assert!(lap.is_empty());
184    }
185
186    #[test]
187    fn single_vertex() {
188        let g = Graph::with_vertices(1);
189        let lap = get_laplacian(
190            &g,
191            DegreeMode::All,
192            LaplacianNormalization::Unnormalized,
193            None,
194        )
195        .unwrap();
196        assert_eq!(lap.len(), 1);
197        assert!(approx_eq(lap[0][0], 0.0));
198    }
199
200    #[test]
201    fn path_3_unnormalized() {
202        let mut g = Graph::with_vertices(3);
203        g.add_edge(0, 1).unwrap();
204        g.add_edge(1, 2).unwrap();
205        let lap = get_laplacian(
206            &g,
207            DegreeMode::All,
208            LaplacianNormalization::Unnormalized,
209            None,
210        )
211        .unwrap();
212        // L = [[ 1, -1,  0],
213        //      [-1,  2, -1],
214        //      [ 0, -1,  1]]
215        assert!(approx_eq(lap[0][0], 1.0));
216        assert!(approx_eq(lap[0][1], -1.0));
217        assert!(approx_eq(lap[0][2], 0.0));
218        assert!(approx_eq(lap[1][0], -1.0));
219        assert!(approx_eq(lap[1][1], 2.0));
220        assert!(approx_eq(lap[1][2], -1.0));
221        assert!(approx_eq(lap[2][0], 0.0));
222        assert!(approx_eq(lap[2][1], -1.0));
223        assert!(approx_eq(lap[2][2], 1.0));
224    }
225
226    #[test]
227    fn row_sum_zero_unnormalized() {
228        let mut g = Graph::with_vertices(4);
229        g.add_edge(0, 1).unwrap();
230        g.add_edge(0, 2).unwrap();
231        g.add_edge(1, 3).unwrap();
232        g.add_edge(2, 3).unwrap();
233        let lap = get_laplacian(
234            &g,
235            DegreeMode::All,
236            LaplacianNormalization::Unnormalized,
237            None,
238        )
239        .unwrap();
240        for row in &lap {
241            let sum: f64 = row.iter().sum();
242            assert!(approx_eq(sum, 0.0));
243        }
244    }
245
246    #[test]
247    fn symmetric_normalization() {
248        let mut g = Graph::with_vertices(3);
249        g.add_edge(0, 1).unwrap();
250        g.add_edge(1, 2).unwrap();
251        let lap =
252            get_laplacian(&g, DegreeMode::All, LaplacianNormalization::Symmetric, None).unwrap();
253        // Diagonal should be 1 for non-isolated vertices
254        assert!(approx_eq(lap[0][0], 1.0));
255        assert!(approx_eq(lap[1][1], 1.0));
256        assert!(approx_eq(lap[2][2], 1.0));
257        // L_sym[0][1] = -1/(sqrt(1)*sqrt(2)) = -1/sqrt(2)
258        assert!(approx_eq(lap[0][1], -1.0 / 2.0_f64.sqrt()));
259    }
260
261    #[test]
262    fn left_normalization() {
263        let mut g = Graph::with_vertices(3);
264        g.add_edge(0, 1).unwrap();
265        g.add_edge(1, 2).unwrap();
266        let lap = get_laplacian(&g, DegreeMode::All, LaplacianNormalization::Left, None).unwrap();
267        // Rows should sum to 0
268        for row in &lap {
269            let sum: f64 = row.iter().sum();
270            assert!(approx_eq(sum, 0.0));
271        }
272        // Diagonal is 1
273        assert!(approx_eq(lap[0][0], 1.0));
274        assert!(approx_eq(lap[1][1], 1.0));
275    }
276
277    #[test]
278    fn weighted() {
279        let mut g = Graph::with_vertices(3);
280        g.add_edge(0, 1).unwrap();
281        g.add_edge(1, 2).unwrap();
282        let weights = vec![2.0, 3.0];
283        let lap = get_laplacian(
284            &g,
285            DegreeMode::All,
286            LaplacianNormalization::Unnormalized,
287            Some(&weights),
288        )
289        .unwrap();
290        // deg(0) = 2, deg(1) = 5, deg(2) = 3
291        assert!(approx_eq(lap[0][0], 2.0));
292        assert!(approx_eq(lap[1][1], 5.0));
293        assert!(approx_eq(lap[2][2], 3.0));
294        assert!(approx_eq(lap[0][1], -2.0));
295        assert!(approx_eq(lap[1][2], -3.0));
296    }
297
298    #[test]
299    fn directed_out_mode() {
300        let mut g = Graph::new(3, true).unwrap();
301        g.add_edge(0, 1).unwrap();
302        g.add_edge(0, 2).unwrap();
303        g.add_edge(1, 2).unwrap();
304        let lap = get_laplacian(
305            &g,
306            DegreeMode::Out,
307            LaplacianNormalization::Unnormalized,
308            None,
309        )
310        .unwrap();
311        // out-deg: 0→2, 1→1, 2→0
312        assert!(approx_eq(lap[0][0], 2.0));
313        assert!(approx_eq(lap[1][1], 1.0));
314        assert!(approx_eq(lap[2][2], 0.0));
315        // Off-diagonal: only from→to
316        assert!(approx_eq(lap[0][1], -1.0));
317        assert!(approx_eq(lap[0][2], -1.0));
318        assert!(approx_eq(lap[1][0], 0.0)); // no edge 1→0
319    }
320
321    #[test]
322    fn directed_in_mode() {
323        let mut g = Graph::new(3, true).unwrap();
324        g.add_edge(0, 1).unwrap();
325        g.add_edge(0, 2).unwrap();
326        g.add_edge(1, 2).unwrap();
327        let lap = get_laplacian(
328            &g,
329            DegreeMode::In,
330            LaplacianNormalization::Unnormalized,
331            None,
332        )
333        .unwrap();
334        // in-deg: 0→0, 1→1, 2→2
335        assert!(approx_eq(lap[0][0], 0.0));
336        assert!(approx_eq(lap[1][1], 1.0));
337        assert!(approx_eq(lap[2][2], 2.0));
338    }
339
340    #[test]
341    fn self_loop_undirected() {
342        let mut g = Graph::with_vertices(2);
343        g.add_edge(0, 0).unwrap();
344        g.add_edge(0, 1).unwrap();
345        let lap = get_laplacian(
346            &g,
347            DegreeMode::All,
348            LaplacianNormalization::Unnormalized,
349            None,
350        )
351        .unwrap();
352        // L_ii = d_i - A_ii. deg(0)=3, A[0][0]=2 → L[0][0]=1
353        assert!(approx_eq(lap[0][0], 1.0));
354        assert!(approx_eq(lap[1][1], 1.0));
355        assert!(approx_eq(lap[0][1], -1.0));
356    }
357
358    #[test]
359    fn weight_mismatch() {
360        let mut g = Graph::with_vertices(2);
361        g.add_edge(0, 1).unwrap();
362        let result = get_laplacian(
363            &g,
364            DegreeMode::All,
365            LaplacianNormalization::Unnormalized,
366            Some(&[1.0, 2.0]),
367        );
368        assert!(result.is_err());
369    }
370
371    #[test]
372    fn isolated_vertex_symmetric() {
373        let mut g = Graph::with_vertices(3);
374        g.add_edge(0, 1).unwrap();
375        // vertex 2 is isolated
376        let lap =
377            get_laplacian(&g, DegreeMode::All, LaplacianNormalization::Symmetric, None).unwrap();
378        // Isolated vertex: diagonal = 0, norm_factor = 0
379        assert!(approx_eq(lap[2][2], 0.0));
380    }
381}
rust_igraph/algorithms/properties/laplacian.rs

rust_igraph/algorithms/properties/
laplacian.rs
