rust_igraph/algorithms/paths/
johnson.rs

1//! Johnson's all-pairs shortest paths (ALGO-SP-003).
2//!
3//! Counterpart of `igraph_distances_johnson()` from
4//! `references/igraph/src/paths/johnson.c:83`. Computes the
5//! shortest-path distance matrix between every pair of vertices,
6//! supporting **negative edge weights** in directed graphs.
7//!
8//! Algorithm:
9//!
10//! Fast path: if no edge weight is negative, fall through to
11//! [`dijkstra_distances`] from every source — that is asymptotically
12//! the same (`O(V (V+E) log V)`) and avoids the reweighting overhead.
13//!
14//! Slow path (negative weights present on a *directed* graph): first
15//! compute potentials `h[v]` via a "virtual source" Bellman-Ford.
16//! The standard trick — initialising every vertex's distance to 0
17//! in SPFA is equivalent to attaching a virtual vertex with
18//! zero-weight outgoing edges to every real vertex and running BF
19//! from it (matches upstream's `bfres` row computation). Reweight
20//! each edge `(u, v, w)` to `w' = w + h[u] - h[v]`. By the triangle
21//! inequality on `h`, every `w'` is ≥ 0 (modulo float roundoff,
22//! which we snap to 0). Run [`dijkstra_distances`] from every source
23//! on the reweighted graph to get `d'`. Recover the true distance:
24//! `d[u][v] = d'[u][v] - h[u] + h[v]`.
25//!
26//! Time complexity: `O(V·E + V·(V+E)·log V)` — V Dijkstras plus one
27//! BF.
28//!
29//! Constraints (matching upstream):
30//! - Undirected graphs with negative weights are rejected
31//!   ([`IgraphError::InvalidArgument`]). This isn't a Johnson
32//!   limitation per se — an undirected negative edge is itself a
33//!   2-cycle of negative weight.
34//! - Negative cycles reachable from anywhere are detected by the
35//!   inner BF and surfaced as [`IgraphError::InvalidArgument`].
36
37use std::collections::VecDeque;
38
39use crate::algorithms::paths::dijkstra::dijkstra_distances;
40use crate::core::{Graph, IgraphError, IgraphResult, VertexId};
41
42fn validate_weights(graph: &Graph, weights: &[f64]) -> IgraphResult<()> {
43    let m = graph.ecount();
44    if weights.len() != m {
45        return Err(IgraphError::InvalidArgument(format!(
46            "weights vector size ({}) differs from edge count ({})",
47            weights.len(),
48            m
49        )));
50    }
51    for (e, &w) in weights.iter().enumerate() {
52        if w.is_nan() {
53            return Err(IgraphError::InvalidArgument(format!(
54                "weight at edge {e} is NaN"
55            )));
56        }
57    }
58    Ok(())
59}
60
61/// SPFA with every vertex initialised to distance 0. Returns the
62/// Johnson potential vector `h[v]`. Equivalent to running
63/// Bellman-Ford from a virtual source with zero-weight outgoing
64/// edges to every real vertex.
65///
66/// Surfaces a negative cycle (any vertex popped more than `vcount`
67/// times) as [`IgraphError::InvalidArgument`].
68fn compute_potentials(graph: &Graph, weights: &[f64]) -> IgraphResult<Vec<f64>> {
69    let n = graph.vcount();
70    let n_usize = n as usize;
71    let mut dist: Vec<f64> = vec![0.0; n_usize];
72    let mut queue: VecDeque<VertexId> = VecDeque::with_capacity(n_usize);
73    let mut in_queue: Vec<bool> = vec![true; n_usize];
74    let mut num_queued: Vec<u32> = vec![0; n_usize];
75    for v in 0..n {
76        queue.push_back(v);
77    }
78
79    while let Some(j) = queue.pop_front() {
80        let j_idx = j as usize;
81        in_queue[j_idx] = false;
82        num_queued[j_idx] = num_queued[j_idx]
83            .checked_add(1)
84            .ok_or(IgraphError::Internal("num_queued overflow"))?;
85        if num_queued[j_idx] > n {
86            return Err(IgraphError::InvalidArgument(
87                "negative cycle detected while computing Johnson potentials".to_string(),
88            ));
89        }
90
91        // OUT-mode SPFA: traverse outgoing incident edges.
92        let incidents = graph.incident(j)?;
93        for eid in incidents {
94            let w = weights[eid as usize];
95            if w == f64::INFINITY {
96                continue;
97            }
98            let target = graph.edge_other(eid, j)?;
99            let target_idx = target as usize;
100            let altdist = dist[j_idx] + w;
101            if altdist < dist[target_idx] {
102                dist[target_idx] = altdist;
103                if !in_queue[target_idx] {
104                    in_queue[target_idx] = true;
105                    queue.push_back(target);
106                }
107            }
108        }
109    }
110
111    Ok(dist)
112}
113
114/// All-pairs shortest-path distances via Johnson's algorithm.
115///
116/// Returns a `vcount × vcount` matrix where `result[u][v]` is the
117/// distance from `u` to `v`: `Some(d)` if reachable, `None`
118/// otherwise. `result[u][u]` is `Some(0.0)` for every valid `u`.
119///
120/// Use this when (a) you need many shortest paths and (b) at least
121/// one edge weight is negative. For non-negative weights this
122/// function transparently falls through to V independent Dijkstra
123/// runs — slightly wasteful but always correct.
124///
125/// `weights[e]` is the weight of edge `e`; length must equal
126/// `graph.ecount()`. NaN weights are rejected. Positive-infinite
127/// weights are treated as "edge ignored" (matches upstream).
128///
129/// **Constraint**: undirected graphs with any negative weight are
130/// rejected ([`IgraphError::InvalidArgument`]) — an undirected
131/// negative edge is itself a length-2 negative cycle.
132///
133/// Counterpart of `igraph_distances_johnson(_, _, vss_all(),
134/// vss_all(), weights, IGRAPH_OUT)`.
135///
136/// # Examples
137///
138/// ```
139/// use rust_igraph::{Graph, johnson_distances};
140///
141/// // Directed diamond 0→1 (3), 0→2 (1), 1→3 (-2), 2→3 (4).
142/// // The negative edge 1→3 makes 0→1→3 cheaper than 0→2→3.
143/// let mut g = Graph::new(4, true).unwrap();
144/// g.add_edge(0, 1).unwrap();
145/// g.add_edge(0, 2).unwrap();
146/// g.add_edge(1, 3).unwrap();
147/// g.add_edge(2, 3).unwrap();
148/// let d = johnson_distances(&g, &[3.0, 1.0, -2.0, 4.0]).unwrap();
149/// // Row 0 (distances from 0): [0, 3, 1, 1]
150/// assert_eq!(d[0], vec![Some(0.0), Some(3.0), Some(1.0), Some(1.0)]);
151/// ```
152pub fn johnson_distances(graph: &Graph, weights: &[f64]) -> IgraphResult<Vec<Vec<Option<f64>>>> {
153    validate_weights(graph, weights)?;
154    let vcount = graph.vcount();
155
156    let has_negative = weights.iter().any(|&w| w < 0.0);
157
158    // Fast path: no negative weights ⇒ Dijkstra from every source.
159    if !has_negative {
160        let mut out = Vec::with_capacity(vcount as usize);
161        for v in 0..vcount {
162            out.push(dijkstra_distances(graph, v, weights)?);
163        }
164        return Ok(out);
165    }
166
167    // Negative weights but undirected: an undirected negative edge
168    // is a length-2 cycle of negative weight.
169    if !graph.is_directed() {
170        return Err(IgraphError::InvalidArgument(
171            "Johnson's algorithm cannot be applied to undirected graphs with negative weights (would form a negative cycle)".to_string(),
172        ));
173    }
174
175    // 1. Compute potentials h[v] via virtual-source SPFA.
176    let potentials = compute_potentials(graph, weights)?;
177
178    // 2. Reweight: w'(u, v) = w(u, v) + h[u] - h[v]. Snap roundoff
179    //    negatives to 0 (matches upstream johnson.c:194).
180    let ecount = graph.ecount();
181    let mut reweighted: Vec<f64> = Vec::with_capacity(ecount);
182    for (e, &orig_w) in weights.iter().enumerate() {
183        let eid = u32::try_from(e)
184            .map_err(|_| IgraphError::Internal("edge id exceeds u32::MAX in johnson"))?;
185        let (from, to) = graph.edge(eid)?;
186        let mut new_w = orig_w + potentials[from as usize] - potentials[to as usize];
187        if new_w < 0.0 {
188            new_w = 0.0;
189        }
190        reweighted.push(new_w);
191    }
192
193    // 3. Run Dijkstra from every source on the reweighted graph.
194    // 4. Recover: d[u][v] = d'[u][v] - h[u] + h[v].
195    let mut out = Vec::with_capacity(vcount as usize);
196    for src in 0..vcount {
197        let dprime = dijkstra_distances(graph, src, &reweighted)?;
198        let mut row: Vec<Option<f64>> = Vec::with_capacity(vcount as usize);
199        for (target, dp) in dprime.into_iter().enumerate() {
200            let recovered = dp.map(|d| d - potentials[src as usize] + potentials[target]);
201            row.push(recovered);
202        }
203        out.push(row);
204    }
205    Ok(out)
206}
207
208#[cfg(test)]
209mod tests {
210    use super::*;
211
212    #[test]
213    fn empty_graph_returns_empty_matrix() {
214        let g = Graph::with_vertices(0);
215        let d = johnson_distances(&g, &[]).unwrap();
216        assert!(d.is_empty());
217    }
218
219    #[test]
220    fn no_edges_diagonal_zero_off_diagonal_none() {
221        let g = Graph::with_vertices(3);
222        let d = johnson_distances(&g, &[]).unwrap();
223        assert_eq!(d.len(), 3);
224        for (u, row) in d.iter().enumerate() {
225            for (v, entry) in row.iter().enumerate() {
226                if u == v {
227                    assert_eq!(entry, &Some(0.0));
228                } else {
229                    assert_eq!(entry, &None);
230                }
231            }
232        }
233    }
234
235    #[test]
236    fn positive_weights_match_pairwise_dijkstra() {
237        // Triangle 0-1-2 with positive weights — both fast paths.
238        let mut g = Graph::with_vertices(3);
239        g.add_edge(0, 1).unwrap();
240        g.add_edge(0, 2).unwrap();
241        g.add_edge(1, 2).unwrap();
242        let weights = [1.0, 4.0, 2.0];
243        let d = johnson_distances(&g, &weights).unwrap();
244        // Row 0: 0, 1, 3 (via shortcut)
245        assert_eq!(d[0], vec![Some(0.0), Some(1.0), Some(3.0)]);
246        assert_eq!(d[1], vec![Some(1.0), Some(0.0), Some(2.0)]);
247        assert_eq!(d[2], vec![Some(3.0), Some(2.0), Some(0.0)]);
248    }
249
250    #[test]
251    fn directed_with_negative_edge_diamond() {
252        let mut g = Graph::new(4, true).unwrap();
253        g.add_edge(0, 1).unwrap();
254        g.add_edge(0, 2).unwrap();
255        g.add_edge(1, 3).unwrap();
256        g.add_edge(2, 3).unwrap();
257        let weights = [3.0, 1.0, -2.0, 4.0];
258        let d = johnson_distances(&g, &weights).unwrap();
259        // From 0: 0→1=3, 0→2=1, 0→3=min(3-2, 1+4)=1
260        assert_eq!(d[0], vec![Some(0.0), Some(3.0), Some(1.0), Some(1.0)]);
261        // From 1: only out-edge is 1→3=-2; 0,2 unreachable
262        assert_eq!(d[1], vec![None, Some(0.0), None, Some(-2.0)]);
263        // From 2: only out-edge is 2→3=4
264        assert_eq!(d[2], vec![None, None, Some(0.0), Some(4.0)]);
265        // From 3: sink, only itself
266        assert_eq!(d[3], vec![None, None, None, Some(0.0)]);
267    }
268
269    #[test]
270    fn unreachable_components_yield_none() {
271        let mut g = Graph::with_vertices(4);
272        g.add_edge(0, 1).unwrap();
273        g.add_edge(2, 3).unwrap();
274        let d = johnson_distances(&g, &[1.0, 2.0]).unwrap();
275        // Two components: 0,1 reach each other; 2,3 reach each other; cross None.
276        assert_eq!(d[0][0], Some(0.0));
277        assert_eq!(d[0][1], Some(1.0));
278        assert_eq!(d[0][2], None);
279        assert_eq!(d[0][3], None);
280        assert_eq!(d[2][3], Some(2.0));
281    }
282
283    #[test]
284    fn undirected_with_negative_weight_errors() {
285        let mut g = Graph::with_vertices(2);
286        g.add_edge(0, 1).unwrap();
287        let err = johnson_distances(&g, &[-1.0]).unwrap_err();
288        assert!(matches!(err, IgraphError::InvalidArgument(_)));
289    }
290
291    #[test]
292    fn negative_cycle_in_directed_errors() {
293        let mut g = Graph::new(3, true).unwrap();
294        g.add_edge(0, 1).unwrap();
295        g.add_edge(1, 2).unwrap();
296        g.add_edge(2, 0).unwrap();
297        // 1 + 1 - 3 = -1 < 0 ⇒ negative cycle
298        let err = johnson_distances(&g, &[1.0, 1.0, -3.0]).unwrap_err();
299        assert!(matches!(err, IgraphError::InvalidArgument(_)));
300    }
301
302    #[test]
303    fn weights_size_mismatch_errors() {
304        let mut g = Graph::with_vertices(2);
305        g.add_edge(0, 1).unwrap();
306        let err = johnson_distances(&g, &[1.0, 2.0]).unwrap_err();
307        assert!(matches!(err, IgraphError::InvalidArgument(_)));
308    }
309
310    #[test]
311    fn nan_weight_errors() {
312        let mut g = Graph::with_vertices(2);
313        g.add_edge(0, 1).unwrap();
314        let err = johnson_distances(&g, &[f64::NAN]).unwrap_err();
315        assert!(matches!(err, IgraphError::InvalidArgument(_)));
316    }
317
318    #[test]
319    fn zero_weights_ok() {
320        let mut g = Graph::new(3, true).unwrap();
321        g.add_edge(0, 1).unwrap();
322        g.add_edge(1, 2).unwrap();
323        let d = johnson_distances(&g, &[0.0, 0.0]).unwrap();
324        assert_eq!(d[0], vec![Some(0.0), Some(0.0), Some(0.0)]);
325    }
326
327    #[test]
328    fn directed_negative_edges_no_cycle_complex() {
329        // 4 vertices: 0→1 (5), 0→2 (-3), 1→3 (1), 2→1 (4), 2→3 (2).
330        // From 0: 0→2→3 = -3 + 2 = -1; 0→2→1→3 = -3+4+1 = 2; so 0→3 = -1.
331        // From 0→1: direct 5, or 0→2→1 = -3+4 = 1, so 0→1 = 1.
332        let mut g = Graph::new(4, true).unwrap();
333        g.add_edge(0, 1).unwrap();
334        g.add_edge(0, 2).unwrap();
335        g.add_edge(1, 3).unwrap();
336        g.add_edge(2, 1).unwrap();
337        g.add_edge(2, 3).unwrap();
338        let weights = [5.0, -3.0, 1.0, 4.0, 2.0];
339        let d = johnson_distances(&g, &weights).unwrap();
340        assert_eq!(d[0][0], Some(0.0));
341        assert_eq!(d[0][1], Some(1.0));
342        assert_eq!(d[0][2], Some(-3.0));
343        assert_eq!(d[0][3], Some(-1.0));
344    }
345}
rust_igraph/algorithms/paths/johnson.rs

rust_igraph/algorithms/paths/
johnson.rs
