rust_igraph/algorithms/layout/
gem.rs

1//! GEM force-directed layout (ALGO-LO-007).
2//!
3//! Graph Embedder algorithm by Arne Frick, Andreas Ludwig, Heiko Mehldau,
4//! "A Fast Adaptive Layout Algorithm for Undirected Graphs",
5//! Proc. Graph Drawing 1994, LNCS 894, pp. 388-403, 1995.
6//!
7//! Key features: per-vertex adaptive temperature, oscillation/rotation
8//! detection via impulse history, gravitational barycenter pull.
9
10use crate::core::{Graph, IgraphError, IgraphResult, VertexId};
11
12/// Parameters for the GEM layout algorithm.
13#[derive(Debug, Clone)]
14pub struct GemParams {
15    /// Maximum number of single-vertex iterations.
16    /// A reasonable default is `40 * n * n`.
17    pub maxiter: u32,
18    /// Maximum allowed local temperature. Default: `n` (vertex count).
19    pub temp_max: f64,
20    /// Global temperature threshold for termination. Default: `0.1`.
21    pub temp_min: f64,
22    /// Initial local temperature per vertex. Default: `sqrt(n)`.
23    pub temp_init: f64,
24}
25
26impl GemParams {
27    /// Create parameters with defaults scaled to the given vertex count.
28    ///
29    /// # Examples
30    ///
31    /// ```
32    /// use rust_igraph::GemParams;
33    ///
34    /// let params = GemParams::for_graph(20);
35    /// assert!(params.temp_max > 0.0);
36    /// assert!(params.temp_init > 0.0);
37    /// ```
38    pub fn for_graph(n: u32) -> Self {
39        let nf = f64::from(n);
40        Self {
41            maxiter: 40u32.saturating_mul(n).saturating_mul(n),
42            temp_max: nf,
43            temp_min: 0.1,
44            temp_init: nf.sqrt(),
45        }
46    }
47}
48
49/// Compute the GEM force-directed layout.
50///
51/// Places vertices using adaptive per-vertex temperatures with
52/// gravitational pull toward the barycenter, repulsive forces between
53/// all vertex pairs, and attractive spring forces along edges.
54///
55/// Edge directions are ignored (treated as undirected).
56///
57/// # Arguments
58///
59/// * `graph` — input graph.
60/// * `seed` — optional initial positions. If `None`, random positions
61///   are generated.
62/// * `params` — algorithm parameters (use `GemParams::for_graph(n)`
63///   for sensible defaults).
64///
65/// Returns `[x, y]` positions for each vertex.
66///
67/// # Errors
68///
69/// Returns `InvalidArgument` if temperatures violate
70/// `temp_min <= temp_init <= temp_max`, any temperature is non-positive,
71/// or `seed` length doesn't match vertex count.
72///
73/// # Examples
74///
75/// ```
76/// use rust_igraph::{Graph, layout_gem, GemParams};
77///
78/// let mut g = Graph::with_vertices(6);
79/// g.add_edge(0, 1).unwrap();
80/// g.add_edge(1, 2).unwrap();
81/// g.add_edge(2, 3).unwrap();
82/// g.add_edge(3, 4).unwrap();
83/// g.add_edge(4, 5).unwrap();
84/// g.add_edge(5, 0).unwrap();
85///
86/// let params = GemParams::for_graph(6);
87/// let pos = layout_gem(&g, None, &params).unwrap();
88/// assert_eq!(pos.len(), 6);
89/// assert!(pos.iter().all(|p| p[0].is_finite() && p[1].is_finite()));
90/// ```
91pub fn layout_gem(
92    graph: &Graph,
93    seed: Option<&[[f64; 2]]>,
94    params: &GemParams,
95) -> IgraphResult<Vec<[f64; 2]>> {
96    let n = graph.vcount() as usize;
97    if n == 0 {
98        return Ok(Vec::new());
99    }
100
101    if params.temp_max <= 0.0 {
102        return Err(IgraphError::InvalidArgument(
103            "temp_max must be positive".into(),
104        ));
105    }
106    if params.temp_min <= 0.0 {
107        return Err(IgraphError::InvalidArgument(
108            "temp_min must be positive".into(),
109        ));
110    }
111    if params.temp_init <= 0.0 {
112        return Err(IgraphError::InvalidArgument(
113            "temp_init must be positive".into(),
114        ));
115    }
116    if params.temp_max < params.temp_init || params.temp_init < params.temp_min {
117        return Err(IgraphError::InvalidArgument(
118            "requires temp_min <= temp_init <= temp_max".into(),
119        ));
120    }
121
122    if let Some(s) = seed {
123        if s.len() != n {
124            return Err(IgraphError::InvalidArgument(format!(
125                "seed length ({}) must equal vertex count ({})",
126                s.len(),
127                n
128            )));
129        }
130    }
131
132    // Constants from the paper
133    let elen_des2: f64 = 128.0 * 128.0;
134    let gamma: f64 = 1.0 / 16.0;
135    let alpha_o: f64 = std::f64::consts::PI;
136    let alpha_r: f64 = std::f64::consts::PI / 3.0;
137    let sigma_o: f64 = 1.0 / 3.0;
138    let sigma_r: f64 = 1.0 / (2.0 * n as f64);
139
140    // Compute vertex "mass" (degree * (degree/2 + 1))
141    let mut phi = vec![0.0_f64; n];
142    for v in 0..n {
143        let deg = graph.degree(v as VertexId).unwrap_or(0) as f64;
144        phi[v] = deg * (deg / 2.0 + 1.0);
145    }
146
147    // Initialize positions
148    let mut pos = if let Some(s) = seed {
149        s.to_vec()
150    } else {
151        let width_half = n as f64 * 100.0;
152        let mut rng = SplitMix64::new(42);
153        (0..n)
154            .map(|_| {
155                [
156                    rng.next_uniform() * 2.0 * width_half - width_half,
157                    rng.next_uniform() * 2.0 * width_half - width_half,
158                ]
159            })
160            .collect()
161    };
162
163    // Barycenter
164    let mut barycenter_x: f64 = pos.iter().map(|p| p[0]).sum();
165    let mut barycenter_y: f64 = pos.iter().map(|p| p[1]).sum();
166
167    // Per-vertex state
168    let mut impulse_x = vec![0.0_f64; n];
169    let mut impulse_y = vec![0.0_f64; n];
170    let mut temp = vec![params.temp_init; n];
171    let mut skew_gauge = vec![0.0_f64; n];
172
173    // Permutation for random vertex selection
174    let mut perm: Vec<usize> = (0..n).collect();
175    let mut perm_pointer: usize = 0;
176    let mut rng = SplitMix64::new(123);
177
178    // Build adjacency list (undirected)
179    let mut adj: Vec<Vec<VertexId>> = vec![Vec::new(); n];
180    let ecount = graph.ecount();
181    for eid in 0..ecount as u32 {
182        if let Ok((src, tgt)) = graph.edge(eid) {
183            if src != tgt {
184                adj[src as usize].push(tgt);
185                adj[tgt as usize].push(src);
186            }
187        }
188    }
189
190    let mut temp_global = params.temp_init * n as f64;
191    let mut maxiter = params.maxiter;
192
193    while temp_global > params.temp_min * n as f64 && maxiter > 0 {
194        // Choose vertex v
195        if perm_pointer == 0 {
196            fisher_yates_shuffle(&mut perm, &mut rng);
197            perm_pointer = n;
198        }
199        perm_pointer -= 1;
200        let v = perm[perm_pointer];
201
202        // Gravitational pull toward barycenter
203        let nf = n as f64;
204        let mut px = (barycenter_x / nf - pos[v][0]) * gamma * phi[v];
205        let mut py = (barycenter_y / nf - pos[v][1]) * gamma * phi[v];
206
207        // Random jitter
208        px += rng.next_uniform() * 64.0 - 32.0;
209        py += rng.next_uniform() * 64.0 - 32.0;
210
211        // Repulsive forces from all other vertices
212        for u in 0..n {
213            if u == v {
214                continue;
215            }
216            let dx = pos[v][0] - pos[u][0];
217            let dy = pos[v][1] - pos[u][1];
218            let dist2 = dx * dx + dy * dy;
219            if dist2 != 0.0 {
220                px += dx * elen_des2 / dist2;
221                py += dy * elen_des2 / dist2;
222            }
223        }
224
225        // Attractive forces from neighbors
226        for &u in &adj[v] {
227            let ui = u as usize;
228            let dx = pos[v][0] - pos[ui][0];
229            let dy = pos[v][1] - pos[ui][1];
230            let dist2 = dx * dx + dy * dy;
231            if phi[v] != 0.0 {
232                px -= dx * dist2 / (elen_des2 * phi[v]);
233                py -= dy * dist2 / (elen_des2 * phi[v]);
234            }
235        }
236
237        // Update position
238        if px != 0.0 || py != 0.0 {
239            let plen = (px * px + py * py).sqrt();
240            px *= temp[v] / plen;
241            py *= temp[v] / plen;
242            pos[v][0] += px;
243            pos[v][1] += py;
244            barycenter_x += px;
245            barycenter_y += py;
246        }
247
248        // Temperature adjustment via oscillation/rotation detection
249        let pvx = impulse_x[v];
250        let pvy = impulse_y[v];
251        if pvx != 0.0 || pvy != 0.0 {
252            let beta = (pvy - py).atan2(pvx - px);
253            let sin_beta = beta.sin();
254            let sign_sin_beta = if sin_beta > 0.0 {
255                1.0
256            } else if sin_beta < 0.0 {
257                -1.0
258            } else {
259                0.0
260            };
261            let cos_beta = beta.cos();
262            let abs_cos_beta = cos_beta.abs();
263            let old_temp = temp[v];
264
265            if sin_beta >= (std::f64::consts::FRAC_PI_2 + alpha_r / 2.0).sin() {
266                skew_gauge[v] += sigma_r * sign_sin_beta;
267            }
268            if abs_cos_beta >= (alpha_o / 2.0).cos() {
269                temp[v] *= sigma_o * cos_beta;
270            }
271            temp[v] *= 1.0 - skew_gauge[v].abs();
272            if temp[v] > params.temp_max {
273                temp[v] = params.temp_max;
274            }
275            if temp[v] < 0.0 {
276                temp[v] = 0.0;
277            }
278            impulse_x[v] = px;
279            impulse_y[v] = py;
280            temp_global += temp[v] - old_temp;
281        }
282
283        maxiter -= 1;
284    }
285
286    Ok(pos)
287}
288
289// ═══════════════════════════════════════════════════════════════════
290// Internal RNG (deterministic, no external deps)
291// ═══════════════════════════════════════════════════════════════════
292
293struct SplitMix64 {
294    state: u64,
295}
296
297impl SplitMix64 {
298    fn new(seed: u64) -> Self {
299        Self { state: seed }
300    }
301
302    fn next_u64(&mut self) -> u64 {
303        self.state = self.state.wrapping_add(0x9E37_79B9_7F4A_7C15);
304        let mut z = self.state;
305        z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
306        z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
307        z ^ (z >> 31)
308    }
309
310    fn next_uniform(&mut self) -> f64 {
311        (self.next_u64() >> 11) as f64 / ((1u64 << 53) as f64)
312    }
313}
314
315fn fisher_yates_shuffle(perm: &mut [usize], rng: &mut SplitMix64) {
316    let n = perm.len();
317    for i in (1..n).rev() {
318        let j = (rng.next_u64() as usize) % (i + 1);
319        perm.swap(i, j);
320    }
321}
322
323// ═══════════════════════════════════════════════════════════════════
324// Tests
325// ═══════════════════════════════════════════════════════════════════
326
327#[cfg(test)]
328mod tests {
329    use super::*;
330
331    #[test]
332    fn test_gem_empty() {
333        let g = Graph::with_vertices(0);
334        let params = GemParams::for_graph(0);
335        let pos = layout_gem(&g, None, &params).unwrap();
336        assert!(pos.is_empty());
337    }
338
339    #[test]
340    fn test_gem_single_vertex() {
341        let g = Graph::with_vertices(1);
342        let params = GemParams::for_graph(1);
343        let pos = layout_gem(&g, None, &params).unwrap();
344        assert_eq!(pos.len(), 1);
345        assert!(pos[0][0].is_finite());
346        assert!(pos[0][1].is_finite());
347    }
348
349    #[test]
350    fn test_gem_triangle() {
351        let mut g = Graph::with_vertices(3);
352        g.add_edge(0, 1).unwrap();
353        g.add_edge(1, 2).unwrap();
354        g.add_edge(2, 0).unwrap();
355        let params = GemParams::for_graph(3);
356        let pos = layout_gem(&g, None, &params).unwrap();
357        assert_eq!(pos.len(), 3);
358        for p in &pos {
359            assert!(p[0].is_finite());
360            assert!(p[1].is_finite());
361        }
362    }
363
364    #[test]
365    fn test_gem_path() {
366        let mut g = Graph::with_vertices(5);
367        for i in 0..4 {
368            g.add_edge(i, i + 1).unwrap();
369        }
370        let params = GemParams::for_graph(5);
371        let pos = layout_gem(&g, None, &params).unwrap();
372        assert_eq!(pos.len(), 5);
373        for p in &pos {
374            assert!(p[0].is_finite());
375            assert!(p[1].is_finite());
376        }
377    }
378
379    #[test]
380    fn test_gem_no_overlap() {
381        let mut g = Graph::with_vertices(4);
382        g.add_edge(0, 1).unwrap();
383        g.add_edge(1, 2).unwrap();
384        g.add_edge(2, 3).unwrap();
385        g.add_edge(3, 0).unwrap();
386        let params = GemParams::for_graph(4);
387        let pos = layout_gem(&g, None, &params).unwrap();
388        // Vertices shouldn't all collapse to the same point
389        let mut all_same = true;
390        for i in 1..4 {
391            if (pos[i][0] - pos[0][0]).abs() > 1e-6 || (pos[i][1] - pos[0][1]).abs() > 1e-6 {
392                all_same = false;
393                break;
394            }
395        }
396        assert!(!all_same, "all vertices collapsed to the same point");
397    }
398
399    #[test]
400    fn test_gem_with_seed() {
401        let mut g = Graph::with_vertices(3);
402        g.add_edge(0, 1).unwrap();
403        g.add_edge(1, 2).unwrap();
404        let seed = vec![[0.0, 0.0], [100.0, 0.0], [50.0, 86.6]];
405        let params = GemParams::for_graph(3);
406        let pos = layout_gem(&g, Some(&seed), &params).unwrap();
407        assert_eq!(pos.len(), 3);
408        for p in &pos {
409            assert!(p[0].is_finite());
410            assert!(p[1].is_finite());
411        }
412    }
413
414    #[test]
415    fn test_gem_seed_wrong_length() {
416        let g = Graph::with_vertices(3);
417        let seed = vec![[0.0, 0.0], [1.0, 0.0]];
418        let params = GemParams::for_graph(3);
419        let result = layout_gem(&g, Some(&seed), &params);
420        assert!(result.is_err());
421    }
422
423    #[test]
424    fn test_gem_invalid_temp() {
425        let g = Graph::with_vertices(3);
426        let params = GemParams {
427            maxiter: 100,
428            temp_max: -1.0,
429            temp_min: 0.1,
430            temp_init: 1.0,
431        };
432        assert!(layout_gem(&g, None, &params).is_err());
433
434        let params2 = GemParams {
435            maxiter: 100,
436            temp_max: 10.0,
437            temp_min: 5.0,
438            temp_init: 2.0,
439        };
440        assert!(layout_gem(&g, None, &params2).is_err());
441    }
442
443    #[test]
444    fn test_gem_deterministic() {
445        let mut g = Graph::with_vertices(4);
446        g.add_edge(0, 1).unwrap();
447        g.add_edge(1, 2).unwrap();
448        g.add_edge(2, 3).unwrap();
449        let params = GemParams::for_graph(4);
450        let pos1 = layout_gem(&g, None, &params).unwrap();
451        let pos2 = layout_gem(&g, None, &params).unwrap();
452        for i in 0..4 {
453            assert!((pos1[i][0] - pos2[i][0]).abs() < 1e-10);
454            assert!((pos1[i][1] - pos2[i][1]).abs() < 1e-10);
455        }
456    }
457
458    #[test]
459    fn test_gem_disconnected() {
460        let mut g = Graph::with_vertices(4);
461        g.add_edge(0, 1).unwrap();
462        g.add_edge(2, 3).unwrap();
463        let params = GemParams::for_graph(4);
464        let pos = layout_gem(&g, None, &params).unwrap();
465        assert_eq!(pos.len(), 4);
466        for p in &pos {
467            assert!(p[0].is_finite());
468            assert!(p[1].is_finite());
469        }
470    }
471
472    #[test]
473    fn test_gem_star() {
474        let mut g = Graph::with_vertices(6);
475        for i in 1..6 {
476            g.add_edge(0, i).unwrap();
477        }
478        let params = GemParams {
479            maxiter: 1000,
480            temp_max: 6.0,
481            temp_min: 0.1,
482            temp_init: 2.4,
483        };
484        let pos = layout_gem(&g, None, &params).unwrap();
485        assert_eq!(pos.len(), 6);
486        for p in &pos {
487            assert!(p[0].is_finite());
488            assert!(p[1].is_finite());
489        }
490    }
491}
rust_igraph/algorithms/layout/gem.rs

rust_igraph/algorithms/layout/
gem.rs
