netoptim_rs/
lib.rs

1//! Bellman-Ford algorithms.
2
3pub mod neg_cycle;
4pub mod parametric;
5// pub mod min_cycle_ratio_ai;
6
7use petgraph::prelude::*;
8
9use petgraph::algo::{FloatMeasure, NegativeCycle};
10use petgraph::visit::{
11    IntoEdges, IntoNodeIdentifiers, NodeCount, NodeIndexable, VisitMap, Visitable,
12};
13
14#[derive(Debug, Clone)]
15pub struct Paths<NodeId, EdgeWeight> {
16    pub distances: Vec<EdgeWeight>,
17    pub predecessors: Vec<Option<NodeId>>,
18}
19
20/// \[Generic\] Compute shortest paths from node `source` to all other.
21///
22/// Using the [Bellman–Ford algorithm][bf]; negative edge costs are
23/// permitted, but the graph must not have a cycle of negative weights
24/// (in that case it will return an error).
25///
26/// On success, return one vec with path costs, and another one which points
27/// out the predecessor of a node along a shortest path. The vectors
28/// are indexed by the graph's node indices.
29///
30/// [bf]: https://en.wikipedia.org/wiki/Bellman%E2%80%93Ford_algorithm
31///
32/// # Example
33/// ```rust
34/// use petgraph::Graph;
35/// use petgraph::algo::bellman_ford;
36/// use petgraph::prelude::*;
37///
38/// let mut g = Graph::new();
39/// let a = g.add_node(()); // node with no weight
40/// let b = g.add_node(());
41/// let c = g.add_node(());
42/// let d = g.add_node(());
43/// let edge = g.add_node(());
44/// let f = g.add_node(());
45/// g.extend_with_edges(&[
46///     (0, 1, 2.0),
47///     (0, 3, 4.0),
48///     (1, 2, 1.0),
49///     (1, 5, 7.0),
50///     (2, 4, 5.0),
51///     (4, 5, 1.0),
52///     (3, 4, 1.0),
53/// ]);
54///
55/// // Graph represented with the weight of each edge
56/// //
57/// //     2       1
58/// // a ----- b ----- c
59/// // | 4     | 7     |
60/// // d       f       | 5
61/// // | 1     | 1     |
62/// // \------ edge ------/
63///
64/// let path = bellman_ford(&g, a);
65/// assert!(path.is_ok());
66/// let path = path.unwrap();
67/// assert_eq!(path.distances, vec![    0.0,     2.0,    3.0,    4.0,     5.0,     6.0]);
68/// assert_eq!(path.predecessors, vec![None, Some(a),Some(b),Some(a), Some(d), Some(edge)]);
69///
70/// // Node f (indice 5) can be reach from a with a path costing 6.
71/// // Predecessor of f is Some(edge) which predecessor is Some(d) which predecessor is Some(a).
72/// // Thus the path from a to f is a <-> d <-> edge <-> f
73///
74/// let graph_with_neg_cycle = Graph::<(), f32, Undirected>::from_edges(&[
75///         (0, 1, -2.0),
76///         (0, 3, -4.0),
77///         (1, 2, -1.0),
78///         (1, 5, -25.0),
79///         (2, 4, -5.0),
80///         (4, 5, -25.0),
81///         (3, 4, -1.0),
82/// ]);
83///
84/// assert!(bellman_ford(&graph_with_neg_cycle, NodeIndex::new(0)).is_err());
85/// ```
86pub fn bellman_ford<G>(
87    g: G,
88    source: G::NodeId,
89) -> Result<Paths<G::NodeId, G::EdgeWeight>, NegativeCycle>
90where
91    G: NodeCount + IntoNodeIdentifiers + IntoEdges + NodeIndexable,
92    G::EdgeWeight: FloatMeasure,
93{
94    let ix = |i| g.to_index(i);
95
96    // Step 1 and Step 2: initialize and relax
97    let (distances, predecessors) = bellman_ford_initialize_relax(g, source);
98
99    // Step 3: check for negative weight cycle
100    for i in g.node_identifiers() {
101        for edge in g.edges(i) {
102            let j = edge.target();
103            let w = *edge.weight();
104            if distances[ix(i)] + w < distances[ix(j)] {
105                return Err(NegativeCycle(()));
106            }
107        }
108    }
109
110    Ok(Paths {
111        distances,
112        predecessors,
113    })
114}
115
116/// \[Generic\] Find the path of a negative cycle reachable from node `source`.
117///
118/// Using the [find_negative_cycle][nc]; will search the Graph for negative cycles using
119/// [Bellman–Ford algorithm][bf]. If no negative cycle is found the function will return `None`.
120///
121/// If a negative cycle is found from source, return one vec with a path of `NodeId`s.
122///
123/// The time complexity of this algorithm should be the same as the Bellman-Ford (O(|V|·|E|)).
124///
125/// [nc]: https://blogs.asarkar.com/assets/docs/algorithms-curated/Negative-Weight%20Cycle%20Algorithms%20-%20Huang.pdf
126/// [bf]: https://en.wikipedia.org/wiki/Bellman%E2%80%93Ford_algorithm
127///
128/// # Example
129/// ```rust
130/// use petgraph::Graph;
131/// use petgraph::algo::find_negative_cycle;
132/// use petgraph::prelude::*;
133///
134/// let graph_with_neg_cycle = Graph::<(), f32, Directed>::from_edges(&[
135///         (0, 1, 1.),
136///         (0, 2, 1.),
137///         (0, 3, 1.),
138///         (1, 3, 1.),
139///         (2, 1, 1.),
140///         (3, 2, -3.),
141/// ]);
142///
143/// let path = find_negative_cycle(&graph_with_neg_cycle, NodeIndex::new(0));
144/// assert_eq!(
145///     path,
146///     Some([NodeIndex::new(1), NodeIndex::new(3), NodeIndex::new(2)].to_vec())
147/// );
148/// ```
149pub fn find_negative_cycle<G>(g: G, source: G::NodeId) -> Option<Vec<G::NodeId>>
150where
151    G: NodeCount + IntoNodeIdentifiers + IntoEdges + NodeIndexable + Visitable,
152    G::EdgeWeight: FloatMeasure,
153{
154    let ix = |i| g.to_index(i);
155    let mut path = Vec::<G::NodeId>::new();
156
157    // Step 1: initialize and relax
158    let (distance, predecessor) = bellman_ford_initialize_relax(g, source);
159
160    // Step 2: Check for negative weight cycle
161    'outer: for i in g.node_identifiers() {
162        for edge in g.edges(i) {
163            let j = edge.target();
164            let w = *edge.weight();
165            if distance[ix(i)] + w < distance[ix(j)] {
166                // Step 3: negative cycle found
167                let start = j;
168                let mut node = start;
169                let mut visited = g.visit_map();
170                // Go backward in the predecessor chain
171                loop {
172                    let ancestor = match predecessor[ix(node)] {
173                        Some(predecessor_node) => predecessor_node,
174                        None => node, // no predecessor, self cycle
175                    };
176                    // We have only 2 ways to find the cycle and break the loop:
177                    // 1. start is reached
178                    if ancestor == start {
179                        path.push(ancestor);
180                        break;
181                    }
182                    // 2. some node was reached twice
183                    else if visited.is_visited(&ancestor) {
184                        // Drop any node in path that is before the first ancestor
185                        let pos = path
186                            .iter()
187                            .position(|&p| p == ancestor)
188                            .expect("we should always have a position");
189                        path = path[pos..path.len()].to_vec();
190
191                        break;
192                    }
193
194                    // None of the above, some middle path node
195                    path.push(ancestor);
196                    visited.visit(ancestor);
197                    node = ancestor;
198                }
199                // We are done here
200                break 'outer;
201            }
202        }
203    }
204    if !path.is_empty() {
205        // Users will probably need to follow the path of the negative cycle
206        // so it should be in the reverse order than it was found by the algorithm.
207        path.reverse();
208        Some(path)
209    } else {
210        None
211    }
212}
213
214// Perform Step 1 and Step 2 of the Bellman-Ford algorithm.
215#[inline(always)]
216fn bellman_ford_initialize_relax<G>(
217    g: G,
218    source: G::NodeId,
219) -> (Vec<G::EdgeWeight>, Vec<Option<G::NodeId>>)
220where
221    G: NodeCount + IntoNodeIdentifiers + IntoEdges + NodeIndexable,
222    G::EdgeWeight: FloatMeasure,
223{
224    // Step 1: initialize graph
225    let mut predecessor = vec![None; g.node_bound()];
226    let mut distance = vec![<_>::infinite(); g.node_bound()];
227    let ix = |i| g.to_index(i);
228    distance[ix(source)] = <_>::zero();
229
230    // Step 2: relax edges repeatedly
231    for _ in 1..g.node_count() {
232        let mut did_update = false;
233        for i in g.node_identifiers() {
234            for edge in g.edges(i) {
235                let j = edge.target();
236                let w = *edge.weight();
237                if distance[ix(i)] + w < distance[ix(j)] {
238                    distance[ix(j)] = distance[ix(i)] + w;
239                    predecessor[ix(j)] = Some(i);
240                    did_update = true;
241                }
242            }
243        }
244        if !did_update {
245            break;
246        }
247    }
248    (distance, predecessor)
249}