|
| 1 | +// Licensed under the Apache License, Version 2.0 (the "License"); you may |
| 2 | +// not use this file except in compliance with the License. You may obtain |
| 3 | +// a copy of the License at |
| 4 | +// |
| 5 | +// http://www.apache.org/licenses/LICENSE-2.0 |
| 6 | +// |
| 7 | +// Unless required by applicable law or agreed to in writing, software |
| 8 | +// distributed under the License is distributed on an "AS IS" BASIS, WITHOUT |
| 9 | +// WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the |
| 10 | +// License for the specific language governing permissions and limitations |
| 11 | +// under the License. |
| 12 | +// |
| 13 | +// This Library is for finding out the shortest paths in increasing cost order. |
| 14 | + |
| 15 | +use crate::petgraph::algo::Measure; |
| 16 | +use crate::petgraph::graph::{Graph, Node, NodeIndex}; |
| 17 | +use crate::petgraph::visit::{EdgeRef, IntoEdgeReferences, IntoEdges, VisitMap, Visitable}; |
| 18 | +use crate::petgraph::EdgeType; |
| 19 | + |
| 20 | +use std::collections::hash_map::Entry::{Occupied, Vacant}; |
| 21 | +use std::collections::{BinaryHeap, HashMap, HashSet}; |
| 22 | +use std::fmt::Debug; |
| 23 | +use std::hash::Hash; |
| 24 | + |
| 25 | +use crate::min_scored::MinScored; |
| 26 | +use petgraph::data::DataMap; |
| 27 | +use petgraph::graph::{Edge, EdgeIndex, EdgeReference, GraphIndex}; |
| 28 | +use petgraph::visit::{GraphBase, GraphRef}; |
| 29 | +use rand::Rng; |
| 30 | + |
| 31 | +#[derive(Debug)] |
| 32 | +struct NodeData<N, K> { |
| 33 | + pub node_id: N, // node identifier |
| 34 | + pub parent_id: Option<N>, // parent identifier |
| 35 | + pub cost: K, // cost of minimum path from source |
| 36 | +} |
| 37 | + |
| 38 | +pub fn dijkstra_shortest_path_with_excluded_prefix<G, F, K>( |
| 39 | + graph: G, |
| 40 | + start: G::NodeId, |
| 41 | + target: G::NodeId, |
| 42 | + mut edge_cost: F, |
| 43 | + excluded_prefix: &HashSet<G::NodeId>, |
| 44 | +) -> (Vec<G::NodeId>, K) |
| 45 | +where |
| 46 | + G: Visitable + IntoEdges, |
| 47 | + G::NodeId: Eq + Hash, |
| 48 | + F: FnMut(G::EdgeRef) -> K, |
| 49 | + K: Measure + Copy, |
| 50 | + G::NodeId: Debug, |
| 51 | + <G as IntoEdgeReferences>::EdgeRef: PartialEq, |
| 52 | +{ |
| 53 | + let mut visit_map = graph.visit_map(); |
| 54 | + let mut scores: HashMap<G::NodeId, K> = HashMap::new(); |
| 55 | + let mut pathnodes: HashMap<G::NodeId, NodeData<G::NodeId, K>> = HashMap::new(); |
| 56 | + let mut visit_next: BinaryHeap<MinScored<K, G::NodeId>> = BinaryHeap::new(); |
| 57 | + visit_next.push(MinScored(K::default(), start)); |
| 58 | + scores.insert(start, K::default()); |
| 59 | + pathnodes.insert( |
| 60 | + start, |
| 61 | + NodeData { |
| 62 | + node_id: start, |
| 63 | + parent_id: None, |
| 64 | + cost: K::default(), |
| 65 | + }, |
| 66 | + ); |
| 67 | + |
| 68 | + // In the loop below, all the nodes which have been assigned the shortest path from source |
| 69 | + // are marked as visited. The visisted nodes are not present in the heap, and hence we do not |
| 70 | + // consider edges to visited nodes. |
| 71 | + while let Some(MinScored(node_score, node)) = visit_next.pop() { |
| 72 | + visit_map.visit(node); |
| 73 | + |
| 74 | + // traverse the unvisited neighbors of node. |
| 75 | + for edge in graph.edges(node) { |
| 76 | + let v = edge.target(); |
| 77 | + |
| 78 | + // don't traverse the nodes marked for exclusion, or which have been visited. |
| 79 | + if visit_map.is_visited(&v) || excluded_prefix.contains(&v) { |
| 80 | + continue; |
| 81 | + } |
| 82 | + |
| 83 | + let edge_weight = edge_cost(edge); |
| 84 | + match scores.entry(v) { |
| 85 | + Occupied(mut ent) => { |
| 86 | + if node_score + edge_weight < *ent.get() { |
| 87 | + // the node leads to shorter path to v. We update the score in the priority heap |
| 88 | + // Since this parent defines the new shortest path, we update the entry in pathnodes |
| 89 | + // with this parent, discarding any earlier entry with other parent(s). |
| 90 | + scores.insert(v, node_score + edge_weight); |
| 91 | + visit_next.push(MinScored(node_score + edge_weight, v)); |
| 92 | + if let Some(v_pnode) = pathnodes.get(&v) { |
| 93 | + let new_node: NodeData<G::NodeId, K> = NodeData { |
| 94 | + node_id: v_pnode.node_id, |
| 95 | + parent_id: Some(node), |
| 96 | + cost: node_score + edge_weight, |
| 97 | + }; |
| 98 | + pathnodes.insert(v_pnode.node_id, new_node); |
| 99 | + } else { |
| 100 | + assert_eq!( |
| 101 | + true, false, |
| 102 | + "Invariant not satisfied, Node {:?} not present in pathnodes", |
| 103 | + v |
| 104 | + ); |
| 105 | + } |
| 106 | + } |
| 107 | + } |
| 108 | + Vacant(_) => { |
| 109 | + // the node v has no entry in the priority queue so far. |
| 110 | + // We must be visiting it the first time. |
| 111 | + // Create the entry in the priority queue and an entry in the pathnodes map. |
| 112 | + scores.insert(v, node_score + edge_weight); |
| 113 | + visit_next.push(MinScored(node_score + edge_weight, v)); |
| 114 | + pathnodes.insert( |
| 115 | + v, |
| 116 | + NodeData { |
| 117 | + node_id: v, |
| 118 | + parent_id: Some(node), |
| 119 | + cost: node_score + edge_weight, |
| 120 | + }, |
| 121 | + ); |
| 122 | + } |
| 123 | + } |
| 124 | + } |
| 125 | + } |
| 126 | + |
| 127 | + // Now return the path |
| 128 | + let mut shortest_path: Vec<G::NodeId> = Vec::new(); |
| 129 | + if let Some(target_node) = pathnodes.get(&target) { |
| 130 | + let min_cost = pathnodes.get(&target).unwrap().cost; |
| 131 | + |
| 132 | + // Let us just return one shortest path |
| 133 | + let mut current_node = target_node; |
| 134 | + shortest_path.push(current_node.node_id); |
| 135 | + while current_node.node_id != start { |
| 136 | + current_node = pathnodes.get(¤t_node.parent_id.unwrap()).unwrap(); |
| 137 | + shortest_path.push(current_node.node_id); |
| 138 | + } |
| 139 | + shortest_path.reverse(); |
| 140 | + (shortest_path, min_cost) |
| 141 | + } else { |
| 142 | + (vec![], K::default()) |
| 143 | + } |
| 144 | +} |
| 145 | + |
| 146 | +/// Implementation of Yen's Algorithm to find k shortest paths. |
| 147 | +/// More on Yen's algorithm - https://people.csail.mit.edu/minilek/yen_kth_shortest.pdf |
| 148 | +
|
| 149 | +pub fn get_smallest_k_paths_yen<N, K, T>( |
| 150 | + graph: &mut Graph<N, K, T>, |
| 151 | + start: NodeIndex, |
| 152 | + target: NodeIndex, |
| 153 | + max_paths: usize, |
| 154 | +) -> (Vec<Vec<NodeIndex>>) |
| 155 | +where |
| 156 | + K: Measure + Copy, |
| 157 | + T: EdgeType, |
| 158 | +{ |
| 159 | + let mut listA: Vec<Vec<NodeIndex>> = Vec::new(); // list to contain shortest paths |
| 160 | + let (shortest_path, min_cost) = dijkstra_shortest_path_with_excluded_prefix( |
| 161 | + &*graph, |
| 162 | + start, |
| 163 | + target, |
| 164 | + |e| *e.weight(), |
| 165 | + &HashSet::new(), |
| 166 | + ); |
| 167 | + |
| 168 | + println!("Inserting path of cost {:?} in listA", min_cost); |
| 169 | + listA.push(shortest_path); |
| 170 | + // A binary heap that contains the candidate paths. In each iteration the candidate paths with |
| 171 | + // their costs are pushed on to the heap. The best path from the heap at the end of the iteration |
| 172 | + // is added to listA. |
| 173 | + let mut listB: BinaryHeap<MinScored<K, Vec<NodeIndex>>> = BinaryHeap::new(); |
| 174 | + let mut paths_in_listB: HashSet<Vec<NodeIndex>> = HashSet::new(); |
| 175 | + |
| 176 | + for i in 1usize..max_paths { |
| 177 | + // listA contains the i shortest paths, while listB contains candidate paths. |
| 178 | + // To determine the (i+1)^th shortest path we proceed as follows (according to Yen's algorithm) |
| 179 | + // We set the last added path in listA, i.e, the i^{th} shortest path as the current path. |
| 180 | + // Let x[0],...,x[ell-1] be the vertices in the current path. |
| 181 | + // We search for (i+1)^th path by considering paths which "diverge" from the current path at one |
| 182 | + // of the nodes x[0],...,x[ell-2]. However, we restrict these paths not to take diverging edge |
| 183 | + // which has appeared in any of the previous paths in listA, which are also identical to current |
| 184 | + // path till the node j. The new path is chosen as the minimum cost path from the following collection. |
| 185 | + // 1. Paths already in list B, |
| 186 | + // 2. The collection of shortest diverging paths from each of the nodes x[0],...x[ell-2]. |
| 187 | + // A diverging path at x[j] is formed as union of current path till node x[j], and the shortest path |
| 188 | + // from x[j] to target after removing prohibited edges (see above). |
| 189 | + |
| 190 | + let mut current_path = listA[i - 1].to_owned(); // current path is the last added path in listA |
| 191 | + let mut excluded_edges_map: HashMap<(NodeIndex, NodeIndex), K> = HashMap::new(); // map to keep track of removed edges, which are added back later |
| 192 | + let mut root_cost = K::default(); // keep track of the cost of current path till diversion point. |
| 193 | + |
| 194 | + for j in 0usize..current_path.len() - 1 { |
| 195 | + // we are looking for diversion at current_path[j] |
| 196 | + let root_node = current_path[j]; |
| 197 | + let next_edge = graph.find_edge(root_node, current_path[j + 1]).unwrap(); |
| 198 | + let next_edge_cost = *graph.edge_weight(next_edge).unwrap(); |
| 199 | + |
| 200 | + for path in listA.iter() { |
| 201 | + // for each path that agrees with current path till node j, |
| 202 | + // remove the neighbor of j^{th} node on that path. |
| 203 | + if path.len() < j + 1 { |
| 204 | + continue; |
| 205 | + } |
| 206 | + |
| 207 | + if path[0..=j] == current_path[0..=j] { |
| 208 | + if let Some(edge) = graph.find_edge(root_node, path[j + 1]) { |
| 209 | + excluded_edges_map |
| 210 | + .insert((root_node, path[j + 1]), *graph.edge_weight(edge).unwrap()); |
| 211 | + graph.remove_edge(edge); |
| 212 | + } |
| 213 | + } |
| 214 | + } |
| 215 | + // find the shortest path form root_node to target in the graph after removing prohibited edges. |
| 216 | + let (shortest_root_target_path, path_cost) = |
| 217 | + dijkstra_shortest_path_with_excluded_prefix( |
| 218 | + &*graph, |
| 219 | + root_node, |
| 220 | + target, |
| 221 | + |e| *e.weight(), |
| 222 | + &HashSet::from_iter(current_path[0..=j].to_vec().into_iter()), |
| 223 | + ); |
| 224 | + |
| 225 | + if shortest_root_target_path.len() > 0 { |
| 226 | + // create new_path by appending current_path till divergence point and the shortest path after that. |
| 227 | + let mut new_path = current_path[0..j].to_owned(); |
| 228 | + new_path.extend_from_slice(&shortest_root_target_path); |
| 229 | + // Add the path to listB if it has not been considered already for inclusion in listB |
| 230 | + if !paths_in_listB.contains(&new_path) { |
| 231 | + listB.push(MinScored(root_cost + path_cost, new_path.clone())); |
| 232 | + paths_in_listB.insert(new_path); |
| 233 | + } |
| 234 | + } |
| 235 | + |
| 236 | + // @todo We should prune listB to size at most max_paths to be more efficient. |
| 237 | + |
| 238 | + // add current edge cost to the root cost |
| 239 | + root_cost = root_cost + next_edge_cost; |
| 240 | + } |
| 241 | + |
| 242 | + // finally restore the edges we removed to reset the graph for next iteration |
| 243 | + for (edge, wt) in &excluded_edges_map { |
| 244 | + graph.add_edge(edge.0, edge.1, *wt); |
| 245 | + } |
| 246 | + |
| 247 | + // remove the path of least cost from listB, and add to listA |
| 248 | + if let Some(MinScored(path_cost, min_path)) = listB.pop() { |
| 249 | + println!("Adding path of cost {:?} to listA", path_cost); |
| 250 | + listA.push(min_path); |
| 251 | + } else { |
| 252 | + // we have run out of candidates paths now. |
| 253 | + return listA; |
| 254 | + } |
| 255 | + } |
| 256 | + |
| 257 | + listA |
| 258 | +} |
| 259 | + |
| 260 | + |
| 261 | +#[cfg(test)] |
| 262 | +mod tests { |
| 263 | + use crate::shortest_path::simple_shortest_paths::get_smallest_k_paths_yen; |
| 264 | + use petgraph::graph::NodeIndex; |
| 265 | + use petgraph::{Graph, Undirected}; |
| 266 | + use rand::Rng; |
| 267 | + |
| 268 | + // The function below generates a graph consisting of n hexagons between two terminal vertices. |
| 269 | + // All edges have weights 1. The illustration below shows the graph for n=2. |
| 270 | + // In general there are 2^n shortest paths from 0 to 6n+1. |
| 271 | + // 2 ----- 3 8 ----- 9 |
| 272 | + // / \ / \ |
| 273 | + // 0 ---- 1 4 ---- 7 10 ----- 13 |
| 274 | + // \ / \ / |
| 275 | + // 5 ----- 6 11 ----- 12 |
| 276 | + // |
| 277 | + fn generate_n_cycle_example(n: usize) -> (Graph<usize, u32, Undirected>, Vec<NodeIndex>) { |
| 278 | + let mut g: Graph<usize, u32, Undirected> = Graph::new_undirected(); |
| 279 | + let num_nodes = 6 * n + 2; |
| 280 | + let mut node_names: Vec<usize> = Vec::new(); |
| 281 | + for i in 0..num_nodes { |
| 282 | + node_names.push(i); |
| 283 | + } |
| 284 | + let mut nodes = (0..num_nodes) |
| 285 | + .into_iter() |
| 286 | + .map(|i| g.add_node(node_names[i])) |
| 287 | + .collect::<Vec<_>>(); |
| 288 | + // build cycles |
| 289 | + for i in 0..n { |
| 290 | + let base = 6 * i + 1; |
| 291 | + for j in 0..6 { |
| 292 | + g.add_edge( |
| 293 | + nodes[base + (j % 6)], |
| 294 | + nodes[base + ((j + 1) % 6)], |
| 295 | + (j + 1) as u32, |
| 296 | + ); |
| 297 | + } |
| 298 | + g.add_edge(nodes[base + 3], nodes[base + 6], 1); |
| 299 | + } |
| 300 | + |
| 301 | + g.add_edge(nodes[0], nodes[1], 1); |
| 302 | + |
| 303 | + (g, nodes) |
| 304 | + } |
| 305 | + |
| 306 | + // Generates an undirected cycle of n edges, each with weight 1. |
| 307 | + fn generate_n_gon_example(n: usize) -> (Graph<usize, u32, Undirected>, Vec<NodeIndex>) { |
| 308 | + let mut g: Graph<usize, u32, Undirected> = Graph::new_undirected(); |
| 309 | + let num_nodes = n; |
| 310 | + let mut node_names: Vec<usize> = Vec::new(); |
| 311 | + for i in 0..num_nodes { |
| 312 | + node_names.push(i); |
| 313 | + } |
| 314 | + let mut nodes = (0..num_nodes) |
| 315 | + .into_iter() |
| 316 | + .map(|i| g.add_node(node_names[i])) |
| 317 | + .collect::<Vec<_>>(); |
| 318 | + // build cycles |
| 319 | + for i in 0..n { |
| 320 | + g.add_edge(nodes[i], nodes[(i + 1) % n], 1); |
| 321 | + } |
| 322 | + (g, nodes) |
| 323 | + } |
| 324 | + |
| 325 | + // This function generates an undirected n-gon as in previous example, with additional chords. |
| 326 | + // The argument step specifies the clockwise distance between vertices connected by the chords. |
| 327 | + fn generate_n_gon_with_chords_example( |
| 328 | + n: usize, |
| 329 | + step: usize, |
| 330 | + ) -> (Graph<usize, u32, Undirected>, Vec<NodeIndex>) { |
| 331 | + let mut g: Graph<usize, u32, Undirected> = Graph::new_undirected(); |
| 332 | + let num_nodes = n; |
| 333 | + let mut node_names: Vec<usize> = Vec::new(); |
| 334 | + for i in 0..num_nodes { |
| 335 | + node_names.push(i); |
| 336 | + } |
| 337 | + let mut nodes = (0..num_nodes) |
| 338 | + .into_iter() |
| 339 | + .map(|i| g.add_node(node_names[i])) |
| 340 | + .collect::<Vec<_>>(); |
| 341 | + // build cycles |
| 342 | + for i in 0..n { |
| 343 | + g.add_edge(nodes[i], nodes[(i + 1) % n], 1); |
| 344 | + g.add_edge(nodes[i], nodes[(i + step) % n], 1); |
| 345 | + } |
| 346 | + (g, nodes) |
| 347 | + } |
| 348 | + |
| 349 | + #[test] |
| 350 | + fn test_k_shortest_paths() { |
| 351 | + let (mut graph, nodes) = generate_n_gon_with_chords_example(6, 2); |
| 352 | + let paths = get_smallest_k_paths_yen(&mut graph, nodes[0], nodes[1], 10); |
| 353 | + for path in paths { |
| 354 | + println!("{:#?}", path); |
| 355 | + } |
| 356 | + } |
| 357 | + |
| 358 | + #[test] |
| 359 | + fn test_k_shortest_random_graph() { |
| 360 | + let mut g = Graph::new_undirected(); |
| 361 | + let nodes: Vec<NodeIndex> = (0..5000).map(|_| g.add_node(())).collect(); |
| 362 | + let mut rng: rand::prelude::ThreadRng = rand::thread_rng(); |
| 363 | + for _ in 0..62291 { |
| 364 | + // Adjust the number of edges as desired |
| 365 | + let a = rng.gen_range(0..nodes.len()); |
| 366 | + let b = rng.gen_range(0..nodes.len()); |
| 367 | + let weight = rng.gen_range(1..100); // Random weight between 1 and 100 |
| 368 | + if a != b { |
| 369 | + // Prevent self-loops |
| 370 | + g.add_edge(nodes[a], nodes[b], weight as f32); |
| 371 | + } |
| 372 | + } |
| 373 | + println!("Graph created"); |
| 374 | + let source = nodes[1]; |
| 375 | + let target = nodes[4000]; |
| 376 | + let shortest_path_get: usize = 5; |
| 377 | + for p in get_smallest_k_paths_yen(&mut g, source, target, 10) { |
| 378 | + println!("Path: {:#?} ", p); |
| 379 | + } |
| 380 | + } |
| 381 | +} |
0 commit comments