sim_lib_discrete_graph/
path.rs

1//! Shortest paths: single-source Dijkstra and Bellman-Ford, checked all-pairs
2//! shortest paths, and reachability over the algebra spine's semiring closure.
3
4use crate::certificate::{ShortestPathCertificate, verify_shortest_paths};
5use crate::error::GraphError;
6use crate::graph::Graph;
7use core::cmp::Reverse;
8use sim_lib_discrete_algebra::{AlgebraLimits, BoolRing, Matrix, MinPlus};
9use std::collections::BinaryHeap;
10
11/// Single-source distances and predecessor forest.
12#[derive(Debug, Clone, PartialEq, Eq)]
13pub struct PathResult<W> {
14    /// `distances[v]` is the shortest distance to `v`, or `None` if unreachable.
15    pub distances: Vec<Option<W>>,
16    /// `predecessors[v]` is the node `v` was reached from on a shortest path.
17    pub predecessors: Vec<Option<usize>>,
18}
19
20/// One shortest path between two nodes, with a verifiable predecessor-tree
21/// certificate for the source.
22#[derive(Debug, Clone, PartialEq, Eq)]
23pub struct ShortestPath<N> {
24    /// Source node.
25    pub source: usize,
26    /// Goal node.
27    pub goal: usize,
28    /// Node labels along the selected shortest path, including endpoints.
29    pub nodes: Vec<N>,
30    /// Total path weight, or `None` when the goal is unreachable.
31    pub distance: Option<i64>,
32    /// The shortest-path tree certificate produced by Bellman-Ford.
33    pub certificate: ShortestPathCertificate,
34}
35
36/// Directed out-arcs `(target, weight)` of `node`, honoring directedness.
37fn out_arcs<N, W: Clone>(graph: &Graph<N, W>, node: usize) -> Vec<(usize, W)> {
38    let undirected = !graph.is_directed();
39    let mut arcs = Vec::new();
40    for e in &graph.edges {
41        if e.source == node {
42            arcs.push((e.target, e.weight.clone()));
43        } else if undirected && e.target == node {
44            arcs.push((e.source, e.weight.clone()));
45        }
46    }
47    arcs
48}
49
50/// Dijkstra's algorithm over non-negative `u64` weights.
51///
52/// # Examples
53///
54/// On a directed graph where the two-hop route `0 -> 1 -> 2` (1 + 2 = 3) beats
55/// the direct edge `0 -> 2` (5), the shortest distance to node `2` is `3`:
56///
57/// ```
58/// use sim_lib_discrete_graph::{dijkstra, Directedness, Graph};
59///
60/// let mut g: Graph<(), u64> = Graph::with_nodes(vec![(), (), ()], Directedness::Directed);
61/// g.add_edge(0, 1, 1).unwrap();
62/// g.add_edge(1, 2, 2).unwrap();
63/// g.add_edge(0, 2, 5).unwrap();
64///
65/// let r = dijkstra(&g, 0).unwrap();
66/// assert_eq!(r.distances, vec![Some(0), Some(1), Some(3)]);
67/// assert_eq!(r.predecessors[2], Some(1)); // reached via node 1
68/// ```
69pub fn dijkstra<N>(graph: &Graph<N, u64>, source: usize) -> Result<PathResult<u64>, GraphError> {
70    graph.validate()?;
71    let n = graph.node_count();
72    if source >= n {
73        return Err(GraphError::NodeOutOfRange {
74            node: source,
75            count: n,
76        });
77    }
78    let mut dist = vec![None; n];
79    let mut pred = vec![None; n];
80    let mut heap: BinaryHeap<Reverse<(u64, usize)>> = BinaryHeap::new();
81    dist[source] = Some(0);
82    heap.push(Reverse((0, source)));
83    while let Some(Reverse((d, u))) = heap.pop() {
84        if dist[u].is_some_and(|best| d > best) {
85            continue;
86        }
87        for (v, w) in out_arcs(graph, u) {
88            // Saturating adversarial weights must not wrap into a spuriously
89            // short distance; an overflowing relaxation is simply no shorter path.
90            let Some(nd) = d.checked_add(w) else {
91                continue;
92            };
93            if dist[v].is_none_or(|best| nd < best) {
94                dist[v] = Some(nd);
95                pred[v] = Some(u);
96                heap.push(Reverse((nd, v)));
97            }
98        }
99    }
100    Ok(PathResult {
101        distances: dist,
102        predecessors: pred,
103    })
104}
105
106/// Bellman-Ford over `i64` weights. Returns the result and whether a
107/// negative-weight cycle is reachable from the source.
108pub fn bellman_ford<N>(
109    graph: &Graph<N, i64>,
110    source: usize,
111) -> Result<(PathResult<i64>, bool), GraphError> {
112    graph.validate()?;
113    let n = graph.node_count();
114    if source >= n {
115        return Err(GraphError::NodeOutOfRange {
116            node: source,
117            count: n,
118        });
119    }
120    let undirected = !graph.is_directed();
121    // Collect all directed arcs once.
122    let mut arcs: Vec<(usize, usize, i64)> = Vec::new();
123    for e in &graph.edges {
124        arcs.push((e.source, e.target, e.weight));
125        if undirected {
126            arcs.push((e.target, e.source, e.weight));
127        }
128    }
129    let mut dist: Vec<Option<i64>> = vec![None; n];
130    let mut pred = vec![None; n];
131    dist[source] = Some(0);
132    for _ in 0..n.saturating_sub(1) {
133        let mut changed = false;
134        for &(a, b, w) in &arcs {
135            if let Some(da) = dist[a] {
136                let nd = da.checked_add(w).ok_or_else(|| {
137                    GraphError::WeightOverflow("Bellman-Ford relaxation".to_string())
138                })?;
139                if dist[b].is_none_or(|best| nd < best) {
140                    dist[b] = Some(nd);
141                    pred[b] = Some(a);
142                    changed = true;
143                }
144            }
145        }
146        if !changed {
147            break;
148        }
149    }
150    let mut negative_cycle = false;
151    for &(a, b, w) in &arcs {
152        if let Some(da) = dist[a] {
153            let nd = da.checked_add(w).ok_or_else(|| {
154                GraphError::WeightOverflow("Bellman-Ford cycle check".to_string())
155            })?;
156            if dist[b].is_none_or(|best| nd < best) {
157                negative_cycle = true;
158                break;
159            }
160        }
161    }
162    Ok((
163        PathResult {
164            distances: dist,
165            predecessors: pred,
166        },
167        negative_cycle,
168    ))
169}
170
171/// Return one shortest path and its reusable certificate.
172///
173/// The helper delegates search to Bellman-Ford and verifies the produced
174/// [`ShortestPathCertificate`] before returning. The graph may be directed or
175/// undirected and may contain negative edges, but negative cycles are rejected.
176///
177/// ```
178/// use sim_lib_discrete_graph::{Directedness, Graph, shortest_path};
179///
180/// let mut g = Graph::with_nodes(vec!["start", "via", "goal"], Directedness::Directed);
181/// g.add_edge(0, 1, 1).unwrap();
182/// g.add_edge(1, 2, 1).unwrap();
183/// g.add_edge(0, 2, 5).unwrap();
184///
185/// let path = shortest_path(&g, 0, 2).unwrap();
186/// assert_eq!(path.nodes, vec!["start", "via", "goal"]);
187/// assert_eq!(path.distance, Some(2));
188/// assert_eq!(path.certificate.predecessors[2], Some(1));
189/// ```
190pub fn shortest_path<N: Clone>(
191    graph: &Graph<N, i64>,
192    source: usize,
193    goal: usize,
194) -> Result<ShortestPath<N>, GraphError> {
195    graph.validate()?;
196    let n = graph.node_count();
197    for node in [source, goal] {
198        if node >= n {
199            return Err(GraphError::NodeOutOfRange { node, count: n });
200        }
201    }
202
203    let (paths, negative_cycle) = bellman_ford(graph, source)?;
204    if negative_cycle {
205        return Err(GraphError::NegativeCycle);
206    }
207    let certificate = ShortestPathCertificate {
208        source,
209        predecessors: paths.predecessors,
210    };
211    verify_shortest_paths(graph, &certificate)?;
212
213    let nodes = if paths.distances[goal].is_some() {
214        let mut reversed = Vec::new();
215        let mut current = goal;
216        loop {
217            reversed.push(graph.nodes[current].clone());
218            if current == source {
219                break;
220            }
221            current = certificate.predecessors[current].ok_or_else(|| {
222                GraphError::CertificateInvalid("path predecessor gap".to_string())
223            })?;
224        }
225        reversed.reverse();
226        reversed
227    } else {
228        Vec::new()
229    };
230
231    Ok(ShortestPath {
232        source,
233        goal,
234        nodes,
235        distance: paths.distances[goal],
236        certificate,
237    })
238}
239
240/// Checked all-pairs shortest paths.
241///
242/// The algebra crate's tropical semiring is a bounded saturating model. This
243/// graph-facing API instead shares Bellman-Ford's fail-closed `i64` overflow
244/// policy so single-source and all-pairs shortest paths agree at numeric
245/// extremes.
246pub fn all_pairs_shortest_paths<N>(graph: &Graph<N, i64>) -> Result<Matrix<MinPlus>, GraphError> {
247    graph.validate()?;
248    let n = graph.node_count();
249    let mut m = Matrix::try_filled_with_limits(n, n, MinPlus::Inf, AlgebraLimits::default())?;
250    for source in 0..n {
251        let (paths, negative_cycle) = bellman_ford(graph, source)?;
252        if negative_cycle {
253            return Err(GraphError::NegativeCycle);
254        }
255        for (target, distance) in paths.distances.into_iter().enumerate() {
256            if let Some(distance) = distance {
257                m.set(source, target, MinPlus::Fin(distance))?;
258            }
259        }
260    }
261    Ok(m)
262}
263
264/// Reachability as the boolean closure of the adjacency matrix. Thin wrapper.
265pub fn reachability<N, W>(graph: &Graph<N, W>) -> Result<Matrix<BoolRing>, GraphError> {
266    graph.validate()?;
267    let n = graph.node_count();
268    let undirected = !graph.is_directed();
269    let mut m = Matrix::try_filled_with_limits(n, n, BoolRing(false), AlgebraLimits::default())?;
270    for e in &graph.edges {
271        m.data[e.source * n + e.target] = BoolRing(true);
272        if undirected {
273            m.data[e.target * n + e.source] = BoolRing(true);
274        }
275    }
276    Ok(m.closure(AlgebraLimits::default())?)
277}
278
279#[cfg(test)]
280mod tests {
281    use super::*;
282    use crate::edge::Directedness;
283
284    #[test]
285    fn dijkstra_row_equals_all_pairs_row() {
286        // Same structure over u64 (Dijkstra) and i64 (all-pairs closure).
287        let edges = [(0usize, 1usize, 1u64), (1, 2, 2), (0, 2, 5), (2, 3, 1)];
288        let mut gu: Graph<u8, u64> = Graph::with_nodes(vec![0, 1, 2, 3], Directedness::Directed);
289        let mut gi: Graph<u8, i64> = Graph::with_nodes(vec![0, 1, 2, 3], Directedness::Directed);
290        for &(s, t, w) in &edges {
291            gu.add_edge(s, t, w).unwrap();
292            gi.add_edge(s, t, w as i64).unwrap();
293        }
294        let dj = dijkstra(&gu, 0).unwrap();
295        let ap = all_pairs_shortest_paths(&gi).unwrap();
296        for j in 0..4 {
297            let from_closure = match ap.data[j] {
298                MinPlus::Fin(d) => Some(d as u64),
299                MinPlus::Inf => None,
300            };
301            assert_eq!(dj.distances[j], from_closure, "node {j}");
302        }
303    }
304
305    #[test]
306    fn bellman_ford_handles_negative_edge_without_cycle() {
307        let mut g: Graph<u8, i64> = Graph::with_nodes(vec![0, 1, 2], Directedness::Directed);
308        g.add_edge(0, 1, 4).unwrap();
309        g.add_edge(0, 2, 5).unwrap();
310        g.add_edge(2, 1, -3).unwrap(); // 0->2->1 = 2 beats direct 4
311        let (res, neg) = bellman_ford(&g, 0).unwrap();
312        assert!(!neg);
313        assert_eq!(res.distances[1], Some(2));
314    }
315
316    #[test]
317    fn bellman_ford_detects_negative_cycle() {
318        let mut g: Graph<u8, i64> = Graph::with_nodes(vec![0, 1], Directedness::Directed);
319        g.add_edge(0, 1, 1).unwrap();
320        g.add_edge(1, 0, -2).unwrap(); // cycle weight -1
321        let (_res, neg) = bellman_ford(&g, 0).unwrap();
322        assert!(neg);
323    }
324
325    #[test]
326    fn near_max_weights_do_not_wrap_distance() {
327        // Dijkstra: two near-u64::MAX hops must not wrap to a tiny distance.
328        let mut gu: Graph<u8, u64> = Graph::with_nodes(vec![0, 1, 2], Directedness::Directed);
329        gu.add_edge(0, 1, u64::MAX - 1).unwrap();
330        gu.add_edge(1, 2, u64::MAX - 1).unwrap();
331        let dj = dijkstra(&gu, 0).unwrap();
332        assert_eq!(dj.distances[1], Some(u64::MAX - 1));
333        // 2 is only reachable via an overflowing relaxation, so it stays unreached.
334        assert_eq!(dj.distances[2], None);
335
336        // Bellman-Ford: two near-i64::MAX hops fail closed instead of wrapping
337        // or silently dropping the overflowing reachable relaxation.
338        let mut gi: Graph<u8, i64> = Graph::with_nodes(vec![0, 1, 2], Directedness::Directed);
339        gi.add_edge(0, 1, i64::MAX - 1).unwrap();
340        gi.add_edge(1, 2, i64::MAX - 1).unwrap();
341        assert!(matches!(
342            bellman_ford(&gi, 0),
343            Err(GraphError::WeightOverflow(_))
344        ));
345    }
346
347    #[test]
348    fn all_pairs_shortest_paths_rejects_overflow() {
349        let mut g: Graph<u8, i64> = Graph::with_nodes(vec![0, 1, 2], Directedness::Directed);
350        g.add_edge(0, 1, i64::MAX - 1).unwrap();
351        g.add_edge(1, 2, i64::MAX - 1).unwrap();
352
353        assert!(matches!(
354            all_pairs_shortest_paths(&g),
355            Err(GraphError::WeightOverflow(_))
356        ));
357    }
358
359    #[test]
360    fn all_pairs_shortest_paths_rejects_negative_cycle() {
361        let mut g: Graph<u8, i64> = Graph::with_nodes(vec![0, 1], Directedness::Directed);
362        g.add_edge(0, 1, 1).unwrap();
363        g.add_edge(1, 0, -2).unwrap();
364
365        assert_eq!(all_pairs_shortest_paths(&g), Err(GraphError::NegativeCycle));
366    }
367
368    #[test]
369    fn bellman_ford_rejects_negative_overflow() {
370        let mut g: Graph<u8, i64> = Graph::with_nodes(vec![0, 1, 2], Directedness::Directed);
371        g.add_edge(0, 1, i64::MIN + 1).unwrap();
372        g.add_edge(1, 2, -2).unwrap();
373
374        assert!(matches!(
375            bellman_ford(&g, 0),
376            Err(GraphError::WeightOverflow(_))
377        ));
378    }
379
380    #[test]
381    fn reachability_is_transitive() {
382        let mut g: Graph<u8, u64> = Graph::with_nodes(vec![0, 1, 2], Directedness::Directed);
383        g.add_edge(0, 1, 1).unwrap();
384        g.add_edge(1, 2, 1).unwrap();
385        let r = reachability(&g).unwrap();
386        assert_eq!(r.data[2], BoolRing(true)); // 0 reaches 2
387        assert_eq!(r.data[6], BoolRing(false)); // 2 does not reach 0
388    }
389
390    #[test]
391    fn shortest_path_returns_verified_certificate() {
392        let mut g = Graph::with_nodes(vec!["start", "via", "goal"], Directedness::Directed);
393        g.add_edge(0, 1, 1).unwrap();
394        g.add_edge(1, 2, 1).unwrap();
395        g.add_edge(0, 2, 5).unwrap();
396
397        let path = shortest_path(&g, 0, 2).unwrap();
398
399        assert_eq!(path.nodes, vec!["start", "via", "goal"]);
400        assert_eq!(path.distance, Some(2));
401        assert_eq!(path.certificate.predecessors, vec![None, Some(0), Some(1)]);
402        verify_shortest_paths(&g, &path.certificate).unwrap();
403    }
404
405    #[test]
406    fn shortest_path_reports_unreachable_goal_with_certificate() {
407        let g = Graph::with_nodes(vec![0, 1], Directedness::Directed);
408
409        let path = shortest_path(&g, 0, 1).unwrap();
410
411        assert_eq!(path.nodes, Vec::<i32>::new());
412        assert_eq!(path.distance, None);
413        assert_eq!(path.certificate.predecessors, vec![None, None]);
414        verify_shortest_paths(&g, &path.certificate).unwrap();
415    }
416}
sim_lib_discrete_graph/path.rs

sim_lib_discrete_graph/
path.rs
