feat: implement bidirectional Dijkstra

examples/dijkstra_bidirectional.rs

@@ -0,0 +1,49 @@
+//! This example demonstrates the bidirectional Dijkstra algorithm, and compares it with the
+//! regular Dijkstra algorithm on a large weighted grid.
+//!
+//! Both searches return the same shortest-path cost, but the bidirectional variant usually settles
+//! fewer nodes: instead of growing a single frontier all the way from the start to the goal, it
+//! grows two smaller frontiers that meet in the middle.
+//!
+//! We search from the centre of the grid to a corner. A single Dijkstra search from the centre
+//! must grow its frontier outwards until it reaches the corner, by which point it has covered
+//! essentially the whole grid. The bidirectional search instead grows one frontier around the
+//! centre and one around the corner; they meet roughly halfway, so together they touch far fewer
+//! cells.
+
+use pathfinding::prelude::{dijkstra, dijkstra_bidirectional};
+use std::time::Instant;
+
+const SIZE: i32 = 400;
+const START: (i32, i32) = (SIZE / 2, SIZE / 2);
+const GOAL: (i32, i32) = (SIZE, SIZE);
+
+/// Moving between two orthogonally adjacent cells costs one plus a small penalty that depends on
+/// both endpoints, so the edge weights are not all identical and Dijkstra cannot be replaced by a
+/// plain BFS. The penalty is symmetric in the two cells, so the grid is undirected and the same
+/// closure describes both successors and predecessors.
+#[expect(clippy::trivially_copy_pass_by_ref)]
+fn neighbours(&(x, y): &(i32, i32)) -> Vec<((i32, i32), usize)> {
+    [(x + 1, y), (x - 1, y), (x, y + 1), (x, y - 1)]
+        .into_iter()
+        .filter(|&(nx, ny)| (0..=SIZE).contains(&nx) && (0..=SIZE).contains(&ny))
+        .map(|(nx, ny)| ((nx, ny), 1 + (x + y + nx + ny).unsigned_abs() as usize % 3))
+        .collect()
+}
+
+fn main() {
+    let instant = Instant::now();
+    let (_, cost) = dijkstra(&START, neighbours, |p| *p == GOAL).expect("no path found");
+    let duration_dijkstra = instant.elapsed();
+
+    let instant = Instant::now();
+    let (_, cost_bidirectional) =
+        dijkstra_bidirectional(&START, &GOAL, neighbours, neighbours).expect("no path found");
+    let duration_bidirectional = instant.elapsed();
+
+    assert_eq!(cost, cost_bidirectional);
+
+    println!("Shortest path cost: {cost}");
+    println!("Dijkstra took {duration_dijkstra:?}");
+    println!("Bidirectional Dijkstra took {duration_bidirectional:?}");
+}

src/directed/dijkstra.rs

@@ -104,6 +104,215 @@ where
     })
 }
 
+/// Compute a shortest path using a bidirectional variant of the [Dijkstra search
+/// algorithm](https://en.wikipedia.org/wiki/Dijkstra's_algorithm).
+///
+/// Two searches are run simultaneously: one forward from `start` following `successors`, and one
+/// backward from `end` following `predecessors`. They progress in order of increasing cost until
+/// they meet in the middle. On large graphs this often settles far fewer nodes than a single
+/// unidirectional search, since two small frontiers are usually cheaper to grow than one large
+/// one.
+///
+/// The shortest path from `start` to `end` is computed and returned along with its total cost, in
+/// a `Some`. If no path can be found, `None` is returned instead.
+///
+/// - `start` is the starting node.
+/// - `end` is the destination node.
+/// - `successors` returns a list of successors for a given node, along with the cost for moving
+///   from the node to the successor. This cost must be non-negative.
+/// - `predecessors` returns a list of predecessors for a given node, along with the cost for
+///   moving from the predecessor to the node. This cost must be non-negative, and for a given edge
+///   it must match the cost reported by `successors`. For an undirected graph, where every edge can
+///   be traversed in both directions at the same cost, the same closure can be used for both
+///   `successors` and `predecessors`.
+///
+/// A node will never be included twice in the path as determined by the `Eq` relationship.
+///
+/// The returned path comprises both the start and end node.
+///
+/// # Example
+///
+/// We search the shortest path on a chess board to go from (1, 1) to (4, 6) doing only knight
+/// moves. Knight moves are symmetrical, so the same closure describes both the successors and the
+/// predecessors of a square.
+///
+/// ```
+/// use pathfinding::prelude::dijkstra_bidirectional;
+///
+/// fn neighbours(&(x, y): &(i32, i32)) -> Vec<((i32, i32), usize)> {
+///     vec![(x+1,y+2), (x+1,y-2), (x-1,y+2), (x-1,y-2),
+///          (x+2,y+1), (x+2,y-1), (x-2,y+1), (x-2,y-1)]
+///         .into_iter().map(|p| (p, 1)).collect()
+/// }
+///
+/// let result = dijkstra_bidirectional(&(1, 1), &(4, 6), neighbours, neighbours);
+/// assert_eq!(result.expect("no path found").1, 4);
+/// ```
+#[expect(clippy::missing_panics_doc)]
+pub fn dijkstra_bidirectional<N, C, FS, IS, FP, IP>(
+    start: &N,
+    end: &N,
+    mut successors: FS,
+    mut predecessors: FP,
+) -> Option<(Vec<N>, C)>
+where
+    N: Eq + Hash + Clone,
+    C: Zero + Ord + Copy,
+    FS: FnMut(&N) -> IS,
+    IS: IntoIterator<Item = (N, C)>,
+    FP: FnMut(&N) -> IP,
+    IP: IntoIterator<Item = (N, C)>,
+{
+    if start == end {
+        return Some((vec![start.clone()], Zero::zero()));
+    }
+
+    let mut forward: FxIndexMap<N, (usize, C)> = FxIndexMap::default();
+    forward.insert(start.clone(), (usize::MAX, Zero::zero()));
+    let mut forward_queue = BinaryHeap::new();
+    forward_queue.push(SmallestHolder {
+        cost: Zero::zero(),
+        index: 0,
+    });
+    let mut forward_settled: FxHashSet<N> = FxHashSet::default();
+
+    let mut backward: FxIndexMap<N, (usize, C)> = FxIndexMap::default();
+    backward.insert(end.clone(), (usize::MAX, Zero::zero()));
+    let mut backward_queue = BinaryHeap::new();
+    backward_queue.push(SmallestHolder {
+        cost: Zero::zero(),
+        index: 0,
+    });
+    let mut backward_settled: FxHashSet<N> = FxHashSet::default();
+
+    // Best complete path found so far, as (total cost, meeting node). The meeting node is present
+    // in both the `forward` and `backward` parent maps, so the full path can be rebuilt from it.
+    let mut best: Option<(C, N)> = None;
+
+    while !forward_queue.is_empty() && !backward_queue.is_empty() {
+        if expand_bidirectional(
+            &mut forward_queue,
+            &mut forward,
+            &mut forward_settled,
+            &backward,
+            &backward_settled,
+            &mut successors,
+            &mut best,
+        ) {
+            break;
+        }
+        if backward_queue.is_empty() {
+            break;
+        }
+        if expand_bidirectional(
+            &mut backward_queue,
+            &mut backward,
+            &mut backward_settled,
+            &forward,
+            &forward_settled,
+            &mut predecessors,
+            &mut best,
+        ) {
+            break;
+        }
+    }
+
+    best.map(|(cost, meeting)| {
+        // The forward half runs from `start` up to the meeting node.
+        let meeting_index = forward.get_index_of(&meeting).unwrap();
+        let mut path = reverse_path(&forward, |&(p, _)| p, meeting_index);
+        // The backward half runs from the meeting node towards `end`, following backward parents.
+        let mut parent = backward.get(&meeting).unwrap().0;
+        while parent != usize::MAX {
+            let (node, &(next, _)) = backward.get_index(parent).unwrap();
+            path.push(node.clone());
+            parent = next;
+        }
+        (path, cost)
+    })
+}
+
+/// Perform a single expansion step of one side of a bidirectional Dijkstra search.
+///
+/// The next unsettled node with the smallest tentative cost is popped from `queue` and settled.
+/// Its neighbours (given by `neighbours`) are relaxed into `parents`, and whenever a neighbour has
+/// already been reached by the opposite search a complete path is available and recorded in `best`
+/// if it improves on the current one.
+///
+/// Returns `true` when the popped node has already been settled by the opposite search, which means
+/// the two frontiers have met and the best recorded path is optimal.
+fn expand_bidirectional<N, C, FN, IN>(
+    queue: &mut BinaryHeap<SmallestHolder<C>>,
+    parents: &mut FxIndexMap<N, (usize, C)>,
+    settled: &mut FxHashSet<N>,
+    opposite: &FxIndexMap<N, (usize, C)>,
+    opposite_settled: &FxHashSet<N>,
+    neighbours: &mut FN,
+    best: &mut Option<(C, N)>,
+) -> bool
+where
+    N: Eq + Hash + Clone,
+    C: Zero + Ord + Copy,
+    FN: FnMut(&N) -> IN,
+    IN: IntoIterator<Item = (N, C)>,
+{
+    let Some(SmallestHolder { cost, index }) = queue.pop() else {
+        return false;
+    };
+    let node = {
+        let (node, &(_, c)) = parents.get_index(index).unwrap();
+        // A cheaper path to this node was found after this entry was queued: discard it.
+        if cost > c {
+            return false;
+        }
+        node.clone()
+    };
+    if !settled.insert(node.clone()) {
+        // Already settled through an equal-cost entry.
+        return false;
+    }
+    if opposite_settled.contains(&node) {
+        // Both searches have settled this node: the best recorded path is optimal.
+        return true;
+    }
+    for (neighbour, move_cost) in neighbours(&node) {
+        let new_cost = cost + move_cost;
+        let n;
+        match parents.entry(neighbour) {
+            Vacant(e) => {
+                n = e.index();
+                e.insert((index, new_cost));
+            }
+            Occupied(mut e) => {
+                if e.get().1 > new_cost {
+                    n = e.index();
+                    e.insert((index, new_cost));
+                } else {
+                    continue;
+                }
+            }
+        }
+        queue.push(SmallestHolder {
+            cost: new_cost,
+            index: n,
+        });
+        // If the opposite search has already reached this neighbour, the two halves form a
+        // complete path; keep it if it is the cheapest one seen so far.
+        let neighbour = parents.get_index(n).unwrap().0;
+        if let Some(&(_, opposite_cost)) = opposite.get(neighbour) {
+            let total = new_cost + opposite_cost;
+            let improved = match best {
+                Some((current, _)) => total < *current,
+                None => true,
+            };
+            if improved {
+                *best = Some((total, neighbour.clone()));
+            }
+        }
+    }
+    false
+}
+
 /// Determine all reachable nodes from a starting point as well as the
 /// minimum cost to reach them and a possible optimal parent node
 /// using the [Dijkstra search

src/lib.rs

@@ -18,7 +18,7 @@
 //! - [Bidirectional search](directed/bfs/fn.bfs_bidirectional.html): simultaneously explore paths forwards from the start and backwards from the goal ([=> Wikipedia][Bidirectional search])
 //! - [Brent](directed/cycle_detection/index.html): find a cycle in an infinite sequence ([⇒ Wikipedia][Brent])
 //! - [DFS](directed/dfs/index.html): explore a graph by going as far as possible, then backtrack ([⇒ Wikipedia][DFS])
-//! - [Dijkstra](directed/dijkstra/index.html): find the shortest path in a weighted graph ([⇒ Wikipedia][Dijkstra])
+//! - [Dijkstra](directed/dijkstra/index.html): find the shortest path in a weighted graph ([⇒ Wikipedia][Dijkstra]), optionally [searching from both ends](directed/dijkstra/fn.dijkstra_bidirectional.html) at once
 //! - [Edmonds Karp](directed/edmonds_karp/index.html): find the maximum flow in a weighted graph ([⇒ Wikipedia][Edmonds Karp])
 //! - [Floyd](directed/cycle_detection/index.html): find a cycle in an infinite sequence ([⇒ Wikipedia][Floyd])
 //! - [Fringe](directed/fringe/index.html): find the shortest path in a weighted graph using an heuristic to guide the process ([⇒ Wikipedia][Fringe])

tests/pathfinding.rs

@@ -42,6 +42,36 @@ mod ex1 {
         }
     }
 
+    #[expect(clippy::trivially_copy_pass_by_ref)]
+    fn predecessors(node: &u8) -> Vec<(u8, usize)> {
+        // Reverse the (directed) successor relation.
+        (0..9u8)
+            .flat_map(|n| successors(&n).filter_map(move |(s, c)| (s == *node).then_some((n, c))))
+            .collect()
+    }
+
+    #[test]
+    fn dijkstra_bidirectional_ok() {
+        for target in 0..9 {
+            let result = dijkstra_bidirectional(&1, &target, successors, predecessors);
+            // Same reachability and same optimal cost as the unidirectional search.
+            let expected_cost = expected(target).map(|(_, c)| c);
+            assert_eq!(result.as_ref().map(|(_, c)| *c), expected_cost);
+            if let Some((path, cost)) = result {
+                // The returned path is a genuine walk from 1 to `target` whose edges sum to `cost`.
+                assert_eq!(path.first(), Some(&1));
+                assert_eq!(path.last(), Some(&target));
+                let walked = path.windows(2).map(|w| {
+                    successors(&w[0])
+                        .find(|(s, _)| *s == w[1])
+                        .expect("path traverses a non-existent edge")
+                        .1
+                });
+                assert_eq!(walked.sum::<usize>(), cost);
+            }
+        }
+    }
+
     #[test]
     fn fringe_ok() {
         for target in 0..9 {
@@ -346,6 +376,21 @@ mod ex2 {
         assert!(path.iter().all(|&(nx, ny)| OPEN[ny][nx]));
     }
 
+    #[test]
+    fn dijkstra_bidirectional_path_ok() {
+        const GOAL: (usize, usize) = (6, 3);
+        // The maze is undirected, so successors and predecessors coincide.
+        let (path, cost) =
+            dijkstra_bidirectional(&(2, 3), &GOAL, successors, successors).expect("path not found");
+        assert_eq!(cost, 8);
+        assert_eq!(path.first(), Some(&(2, 3)));
+        assert_eq!(path.last(), Some(&GOAL));
+        assert!(path.iter().all(|&(nx, ny)| OPEN[ny][nx]));
+        // The result must match the unidirectional Dijkstra shortest path cost.
+        let (_, reference) = dijkstra(&(2, 3), successors, |n| n == &GOAL).unwrap();
+        assert_eq!(cost, reference);
+    }
+
     #[test]
     fn dfs_path_ok() {
         const GOAL: (usize, usize) = (6, 3);
@@ -430,6 +475,15 @@ mod ex2 {
         );
     }
 
+    #[test]
+    fn dijkstra_bidirectional_no_path() {
+        const GOAL: (usize, usize) = (1, 1);
+        assert_eq!(
+            dijkstra_bidirectional(&(2, 3), &GOAL, successors, successors),
+            None
+        );
+    }
+
     #[test]
     fn dfs_no_path() {
         const GOAL: (usize, usize) = (1, 1);