oxicuda_graphalg/shortest_path/
bidijkstra.rs

1//! Bidirectional Dijkstra (meet in the middle).
2
3use std::cmp::Ordering;
4use std::collections::BinaryHeap;
5
6use crate::error::{GraphalgError, GraphalgResult};
7use crate::repr::weighted_graph::WeightedGraph;
8
9#[derive(Debug, Clone)]
10pub struct BiDijkstraOutput {
11    pub dist: f64,
12    pub path: Vec<usize>,
13}
14
15#[derive(Debug, Clone, Copy, PartialEq)]
16struct Item {
17    dist: f64,
18    node: usize,
19}
20impl Eq for Item {}
21impl PartialOrd for Item {
22    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
23        Some(self.cmp(other))
24    }
25}
26impl Ord for Item {
27    fn cmp(&self, other: &Self) -> Ordering {
28        other
29            .dist
30            .partial_cmp(&self.dist)
31            .unwrap_or(Ordering::Equal)
32            .then_with(|| other.node.cmp(&self.node))
33    }
34}
35
36pub fn bidirectional_dijkstra(
37    g: &WeightedGraph,
38    source: usize,
39    target: usize,
40) -> GraphalgResult<BiDijkstraOutput> {
41    if source >= g.n || target >= g.n {
42        return Err(GraphalgError::SourceOutOfRange {
43            node: source.max(target),
44            n: g.n,
45        });
46    }
47    if source == target {
48        return Ok(BiDijkstraOutput {
49            dist: 0.0,
50            path: vec![source],
51        });
52    }
53    let rev = g.reverse();
54    let mut df = vec![f64::INFINITY; g.n];
55    let mut db = vec![f64::INFINITY; g.n];
56    let mut pf = vec![usize::MAX; g.n];
57    let mut pb = vec![usize::MAX; g.n];
58    df[source] = 0.0;
59    db[target] = 0.0;
60    let mut hf: BinaryHeap<Item> = BinaryHeap::new();
61    let mut hb: BinaryHeap<Item> = BinaryHeap::new();
62    hf.push(Item {
63        dist: 0.0,
64        node: source,
65    });
66    hb.push(Item {
67        dist: 0.0,
68        node: target,
69    });
70    let mut best = f64::INFINITY;
71    let mut meet = usize::MAX;
72    while !hf.is_empty() || !hb.is_empty() {
73        // expand smaller frontier
74        let f_top = hf.peek().map(|i| i.dist).unwrap_or(f64::INFINITY);
75        let b_top = hb.peek().map(|i| i.dist).unwrap_or(f64::INFINITY);
76        if f_top + b_top >= best {
77            break;
78        }
79        if f_top <= b_top {
80            if let Some(Item { dist: d, node: u }) = hf.pop() {
81                if d > df[u] {
82                    continue;
83                }
84                for &(v, w) in g.neighbors(u)? {
85                    if w < 0.0 {
86                        return Err(GraphalgError::NegativeWeight {
87                            edge: (u, v),
88                            weight: w,
89                        });
90                    }
91                    let nd = d + w;
92                    if nd < df[v] {
93                        df[v] = nd;
94                        pf[v] = u;
95                        hf.push(Item { dist: nd, node: v });
96                        if db[v].is_finite() && df[v] + db[v] < best {
97                            best = df[v] + db[v];
98                            meet = v;
99                        }
100                    }
101                }
102            }
103        } else if let Some(Item { dist: d, node: u }) = hb.pop() {
104            if d > db[u] {
105                continue;
106            }
107            for &(v, w) in rev.neighbors(u)? {
108                if w < 0.0 {
109                    return Err(GraphalgError::NegativeWeight {
110                        edge: (v, u),
111                        weight: w,
112                    });
113                }
114                let nd = d + w;
115                if nd < db[v] {
116                    db[v] = nd;
117                    pb[v] = u;
118                    hb.push(Item { dist: nd, node: v });
119                    if df[v].is_finite() && df[v] + db[v] < best {
120                        best = df[v] + db[v];
121                        meet = v;
122                    }
123                }
124            }
125        }
126    }
127    if !best.is_finite() {
128        return Err(GraphalgError::NoSolution(format!(
129            "no path from {source} to {target}"
130        )));
131    }
132    let mut left = Vec::new();
133    let mut cur = meet;
134    while cur != source {
135        left.push(cur);
136        let p = pf[cur];
137        if p == usize::MAX {
138            return Err(GraphalgError::NumericalInstability(
139                "forward parent missing".to_string(),
140            ));
141        }
142        cur = p;
143    }
144    left.push(source);
145    left.reverse();
146    let mut cur = meet;
147    while cur != target {
148        let p = pb[cur];
149        if p == usize::MAX {
150            return Err(GraphalgError::NumericalInstability(
151                "backward parent missing".to_string(),
152            ));
153        }
154        cur = p;
155        left.push(cur);
156    }
157    Ok(BiDijkstraOutput {
158        dist: best,
159        path: left,
160    })
161}
162
163#[cfg(test)]
164mod tests {
165    use super::*;
166
167    #[test]
168    fn bidir_dijkstra_path() {
169        let mut g = WeightedGraph::new(4);
170        g.add_edge(0, 1, 1.0).expect("ok");
171        g.add_edge(0, 2, 4.0).expect("ok");
172        g.add_edge(1, 2, 2.0).expect("ok");
173        g.add_edge(1, 3, 5.0).expect("ok");
174        g.add_edge(2, 3, 1.0).expect("ok");
175        let out = bidirectional_dijkstra(&g, 0, 3).expect("ok");
176        assert!((out.dist - 4.0).abs() < 1e-12);
177    }
178}
oxicuda_graphalg/shortest_path/bidijkstra.rs

oxicuda_graphalg/shortest_path/
bidijkstra.rs
