sim_lib_discrete_graph/
path.rs

1//! Shortest paths: single-source Dijkstra and Bellman-Ford, plus all-pairs and
2//! reachability as thin wrappers over the algebra spine's semiring closure.
3
4use crate::error::GraphError;
5use crate::graph::Graph;
6use core::cmp::Reverse;
7use sim_lib_discrete_algebra::{AlgebraLimits, BoolRing, Matrix, MinPlus};
8use std::collections::BinaryHeap;
9
10/// Single-source distances and predecessor forest.
11#[derive(Debug, Clone, PartialEq, Eq)]
12pub struct PathResult<W> {
13    /// `distances[v]` is the shortest distance to `v`, or `None` if unreachable.
14    pub distances: Vec<Option<W>>,
15    /// `predecessors[v]` is the node `v` was reached from on a shortest path.
16    pub predecessors: Vec<Option<usize>>,
17}
18
19/// Directed out-arcs `(target, weight)` of `node`, honoring directedness.
20fn out_arcs<N, W: Clone>(graph: &Graph<N, W>, node: usize) -> Vec<(usize, W)> {
21    let undirected = !graph.is_directed();
22    let mut arcs = Vec::new();
23    for e in &graph.edges {
24        if e.source == node {
25            arcs.push((e.target, e.weight.clone()));
26        } else if undirected && e.target == node {
27            arcs.push((e.source, e.weight.clone()));
28        }
29    }
30    arcs
31}
32
33/// Dijkstra's algorithm over non-negative `u64` weights.
34///
35/// # Examples
36///
37/// On a directed graph where the two-hop route `0 -> 1 -> 2` (1 + 2 = 3) beats
38/// the direct edge `0 -> 2` (5), the shortest distance to node `2` is `3`:
39///
40/// ```
41/// use sim_lib_discrete_graph::{dijkstra, Directedness, Graph};
42///
43/// let mut g: Graph<(), u64> = Graph::with_nodes(vec![(), (), ()], Directedness::Directed);
44/// g.add_edge(0, 1, 1).unwrap();
45/// g.add_edge(1, 2, 2).unwrap();
46/// g.add_edge(0, 2, 5).unwrap();
47///
48/// let r = dijkstra(&g, 0).unwrap();
49/// assert_eq!(r.distances, vec![Some(0), Some(1), Some(3)]);
50/// assert_eq!(r.predecessors[2], Some(1)); // reached via node 1
51/// ```
52pub fn dijkstra<N>(graph: &Graph<N, u64>, source: usize) -> Result<PathResult<u64>, GraphError> {
53    graph.validate()?;
54    let n = graph.node_count();
55    if source >= n {
56        return Err(GraphError::NodeOutOfRange {
57            node: source,
58            count: n,
59        });
60    }
61    let mut dist = vec![None; n];
62    let mut pred = vec![None; n];
63    let mut heap: BinaryHeap<Reverse<(u64, usize)>> = BinaryHeap::new();
64    dist[source] = Some(0);
65    heap.push(Reverse((0, source)));
66    while let Some(Reverse((d, u))) = heap.pop() {
67        if dist[u].is_some_and(|best| d > best) {
68            continue;
69        }
70        for (v, w) in out_arcs(graph, u) {
71            // Saturating adversarial weights must not wrap into a spuriously
72            // short distance; an overflowing relaxation is simply no shorter path.
73            let Some(nd) = d.checked_add(w) else {
74                continue;
75            };
76            if dist[v].is_none_or(|best| nd < best) {
77                dist[v] = Some(nd);
78                pred[v] = Some(u);
79                heap.push(Reverse((nd, v)));
80            }
81        }
82    }
83    Ok(PathResult {
84        distances: dist,
85        predecessors: pred,
86    })
87}
88
89/// Bellman-Ford over `i64` weights. Returns the result and whether a
90/// negative-weight cycle is reachable from the source.
91pub fn bellman_ford<N>(
92    graph: &Graph<N, i64>,
93    source: usize,
94) -> Result<(PathResult<i64>, bool), GraphError> {
95    graph.validate()?;
96    let n = graph.node_count();
97    if source >= n {
98        return Err(GraphError::NodeOutOfRange {
99            node: source,
100            count: n,
101        });
102    }
103    let undirected = !graph.is_directed();
104    // Collect all directed arcs once.
105    let mut arcs: Vec<(usize, usize, i64)> = Vec::new();
106    for e in &graph.edges {
107        arcs.push((e.source, e.target, e.weight));
108        if undirected {
109            arcs.push((e.target, e.source, e.weight));
110        }
111    }
112    let mut dist: Vec<Option<i64>> = vec![None; n];
113    let mut pred = vec![None; n];
114    dist[source] = Some(0);
115    for _ in 0..n.saturating_sub(1) {
116        let mut changed = false;
117        for &(a, b, w) in &arcs {
118            if let Some(da) = dist[a] {
119                // An overflowing relaxation is treated as no shorter path, never
120                // a wrapped (and falsely shorter) distance.
121                let Some(nd) = da.checked_add(w) else {
122                    continue;
123                };
124                if dist[b].is_none_or(|best| nd < best) {
125                    dist[b] = Some(nd);
126                    pred[b] = Some(a);
127                    changed = true;
128                }
129            }
130        }
131        if !changed {
132            break;
133        }
134    }
135    let mut negative_cycle = false;
136    for &(a, b, w) in &arcs {
137        if let Some(da) = dist[a]
138            && let Some(nd) = da.checked_add(w)
139            && dist[b].is_none_or(|best| nd < best)
140        {
141            negative_cycle = true;
142            break;
143        }
144    }
145    Ok((
146        PathResult {
147            distances: dist,
148            predecessors: pred,
149        },
150        negative_cycle,
151    ))
152}
153
154/// All-pairs shortest paths as the min-plus closure of the adjacency matrix.
155/// This is a thin wrapper over the spine; it does not re-implement Floyd-Warshall.
156pub fn all_pairs_shortest_paths<N>(graph: &Graph<N, i64>) -> Result<Matrix<MinPlus>, GraphError> {
157    graph.validate()?;
158    let n = graph.node_count();
159    let undirected = !graph.is_directed();
160    let mut m = Matrix::filled(n, n, MinPlus::Inf);
161    for e in &graph.edges {
162        // Accumulate by semiring add (min) so parallel edges keep the cheapest.
163        m.data[e.source * n + e.target] = min_plus_add(m.data[e.source * n + e.target], e.weight);
164        if undirected {
165            m.data[e.target * n + e.source] =
166                min_plus_add(m.data[e.target * n + e.source], e.weight);
167        }
168    }
169    Ok(m.closure(AlgebraLimits::default())?)
170}
171
172fn min_plus_add(cur: MinPlus, w: i64) -> MinPlus {
173    use sim_lib_discrete_algebra::Semiring;
174    cur.add(&MinPlus::Fin(w))
175}
176
177/// Reachability as the boolean closure of the adjacency matrix. Thin wrapper.
178pub fn reachability<N, W>(graph: &Graph<N, W>) -> Result<Matrix<BoolRing>, GraphError> {
179    graph.validate()?;
180    let n = graph.node_count();
181    let undirected = !graph.is_directed();
182    let mut m = Matrix::filled(n, n, BoolRing(false));
183    for e in &graph.edges {
184        m.data[e.source * n + e.target] = BoolRing(true);
185        if undirected {
186            m.data[e.target * n + e.source] = BoolRing(true);
187        }
188    }
189    Ok(m.closure(AlgebraLimits::default())?)
190}
191
192#[cfg(test)]
193mod tests {
194    use super::*;
195    use crate::edge::Directedness;
196
197    #[test]
198    fn dijkstra_row_equals_all_pairs_row() {
199        // Same structure over u64 (Dijkstra) and i64 (all-pairs closure).
200        let edges = [(0usize, 1usize, 1u64), (1, 2, 2), (0, 2, 5), (2, 3, 1)];
201        let mut gu: Graph<u8, u64> = Graph::with_nodes(vec![0, 1, 2, 3], Directedness::Directed);
202        let mut gi: Graph<u8, i64> = Graph::with_nodes(vec![0, 1, 2, 3], Directedness::Directed);
203        for &(s, t, w) in &edges {
204            gu.add_edge(s, t, w).unwrap();
205            gi.add_edge(s, t, w as i64).unwrap();
206        }
207        let dj = dijkstra(&gu, 0).unwrap();
208        let ap = all_pairs_shortest_paths(&gi).unwrap();
209        for j in 0..4 {
210            let from_closure = match ap.data[j] {
211                MinPlus::Fin(d) => Some(d as u64),
212                MinPlus::Inf => None,
213            };
214            assert_eq!(dj.distances[j], from_closure, "node {j}");
215        }
216    }
217
218    #[test]
219    fn bellman_ford_handles_negative_edge_without_cycle() {
220        let mut g: Graph<u8, i64> = Graph::with_nodes(vec![0, 1, 2], Directedness::Directed);
221        g.add_edge(0, 1, 4).unwrap();
222        g.add_edge(0, 2, 5).unwrap();
223        g.add_edge(2, 1, -3).unwrap(); // 0->2->1 = 2 beats direct 4
224        let (res, neg) = bellman_ford(&g, 0).unwrap();
225        assert!(!neg);
226        assert_eq!(res.distances[1], Some(2));
227    }
228
229    #[test]
230    fn bellman_ford_detects_negative_cycle() {
231        let mut g: Graph<u8, i64> = Graph::with_nodes(vec![0, 1], Directedness::Directed);
232        g.add_edge(0, 1, 1).unwrap();
233        g.add_edge(1, 0, -2).unwrap(); // cycle weight -1
234        let (_res, neg) = bellman_ford(&g, 0).unwrap();
235        assert!(neg);
236    }
237
238    #[test]
239    fn near_max_weights_do_not_wrap_distance() {
240        // Dijkstra: two near-u64::MAX hops must not wrap to a tiny distance.
241        let mut gu: Graph<u8, u64> = Graph::with_nodes(vec![0, 1, 2], Directedness::Directed);
242        gu.add_edge(0, 1, u64::MAX - 1).unwrap();
243        gu.add_edge(1, 2, u64::MAX - 1).unwrap();
244        let dj = dijkstra(&gu, 0).unwrap();
245        assert_eq!(dj.distances[1], Some(u64::MAX - 1));
246        // 2 is only reachable via an overflowing relaxation, so it stays unreached.
247        assert_eq!(dj.distances[2], None);
248
249        // Bellman-Ford: two near-i64::MAX hops likewise must not wrap.
250        let mut gi: Graph<u8, i64> = Graph::with_nodes(vec![0, 1, 2], Directedness::Directed);
251        gi.add_edge(0, 1, i64::MAX - 1).unwrap();
252        gi.add_edge(1, 2, i64::MAX - 1).unwrap();
253        let (res, neg) = bellman_ford(&gi, 0).unwrap();
254        assert!(!neg);
255        assert_eq!(res.distances[1], Some(i64::MAX - 1));
256        assert_eq!(res.distances[2], None);
257    }
258
259    #[test]
260    fn reachability_is_transitive() {
261        let mut g: Graph<u8, u64> = Graph::with_nodes(vec![0, 1, 2], Directedness::Directed);
262        g.add_edge(0, 1, 1).unwrap();
263        g.add_edge(1, 2, 1).unwrap();
264        let r = reachability(&g).unwrap();
265        assert_eq!(r.data[2], BoolRing(true)); // 0 reaches 2
266        assert_eq!(r.data[6], BoolRing(false)); // 2 does not reach 0
267    }
268}
sim_lib_discrete_graph/path.rs

sim_lib_discrete_graph/
path.rs
