toolbox_rs/
ford_fulkerson.rs

1use crate::bfs::BFS;
2use crate::{
3    edge::InputEdge,
4    graph::{Graph, NodeID},
5    max_flow::{MaxFlow, ResidualEdgeData},
6    static_graph::StaticGraph,
7};
8use bitvec::vec::BitVec;
9use itertools::Itertools;
10use log::{debug, warn};
11use std::sync::atomic::{AtomicI32, Ordering};
12use std::sync::Arc;
13use std::time::Instant;
14
15pub struct FordFulkerson {
16    residual_graph: StaticGraph<ResidualEdgeData>,
17    max_flow: i32,
18    finished: bool,
19    source: NodeID,
20    target: NodeID,
21    bound: Option<Arc<AtomicI32>>,
22}
23
24impl MaxFlow for FordFulkerson {
25    fn from_edge_list(
26        mut edge_list: Vec<InputEdge<ResidualEdgeData>>,
27        source: usize,
28        target: usize,
29    ) -> Self {
30        let number_of_edges = edge_list.len();
31
32        debug!("extending {} edges", edge_list.len());
33        // blindly generate reverse edges for all edges with zero capacity
34        edge_list.extend_from_within(..);
35        edge_list.iter_mut().skip(number_of_edges).for_each(|edge| {
36            edge.reverse();
37            edge.data.capacity = 0;
38        });
39        debug!("into {} edges", edge_list.len());
40
41        // dedup-merge edge set, by using the following trick: not the dedup(.) call
42        // below takes the second argument as mut. When deduping equivalent values
43        // a and b, then a is accumulated onto b.
44        edge_list.sort_unstable();
45        edge_list.dedup_by(|a, b| {
46            // edges a and b are assumed to be equivalent in the residual graph if
47            // (and only if) they are parallel. In other words, this removes parallel
48            // edges in the residual graph and accumulates capacities on the remaining
49            // egde.
50            let result = a.source == b.source && a.target == b.target;
51            if result {
52                b.data.capacity += a.data.capacity;
53            }
54            result
55        });
56        debug!("dedup-merged {} edges", edge_list.len());
57
58        // at this point the edge set of the residual graph doesn't have any
59        // duplicates anymore. note that this is fine, as we are looking to
60        // compute a node partition.
61        Self {
62            residual_graph: StaticGraph::new(edge_list),
63            max_flow: 0,
64            finished: false,
65            source,
66            target,
67            bound: None,
68        }
69    }
70
71    fn run_with_upper_bound(&mut self, bound: Arc<AtomicI32>) {
72        warn!("Upper bound {} is discarded", bound.load(Ordering::Relaxed));
73        self.bound = Some(bound);
74        self.run()
75    }
76
77    fn run(&mut self) {
78        let mut bfs = BFS::new(
79            &[self.source],
80            &[self.target],
81            self.residual_graph.number_of_nodes(),
82        );
83        let filter = |graph: &StaticGraph<ResidualEdgeData>, edge| graph.data(edge).capacity <= 0;
84        // let mut iteration = 0;
85        while bfs.run_with_filter(&self.residual_graph, filter) {
86            let start = Instant::now();
87            // retrieve node path. The path is unambiguous, as we removed all duplicate edges
88            // find min capacity on edges of the path
89            let bootleneck_head_tail = bfs
90                .path_iter()
91                .tuple_windows()
92                .min_by_key(|(a, b)| {
93                    let edge = self.residual_graph.find_edge_unchecked(*b, *a);
94                    self.residual_graph.data(edge).capacity
95                })
96                .expect("graph is broken, couldn't find min edge");
97            let duration = start.elapsed();
98            debug!(" flow assignment1 took: {:?} (done)", duration);
99
100            let bottleneck_edge = self
101                .residual_graph
102                .find_edge_unchecked(bootleneck_head_tail.1, bootleneck_head_tail.0);
103            debug!("  bottleneck edge: {}", bottleneck_edge);
104            let path_flow = self.residual_graph.data(bottleneck_edge).capacity;
105            debug_assert!(path_flow > 0);
106            debug!("min edge: {}, capacity: {}", bottleneck_edge, path_flow);
107            // sum up flow
108            self.max_flow += path_flow;
109            let duration = start.elapsed();
110            debug!(" flow assignment2 took: {:?} (done)", duration);
111
112            // assign flow to residual graph
113            for (a, b) in bfs.path_iter().tuple_windows() {
114                let rev_edge = self.residual_graph.find_edge_unchecked(a, b);
115                let fwd_edge = self.residual_graph.find_edge_unchecked(b, a);
116
117                self.residual_graph.data_mut(fwd_edge).capacity -= path_flow;
118                self.residual_graph.data_mut(rev_edge).capacity += path_flow;
119            }
120            let duration = start.elapsed();
121            debug!(" flow assignment3 took: {:?} (done)", duration);
122        }
123
124        self.finished = true;
125    }
126
127    fn max_flow(&self) -> Result<i32, String> {
128        if !self.finished {
129            return Err("Assigment was not computed.".to_string());
130        }
131        Ok(self.max_flow)
132    }
133
134    fn assignment(&self, source: NodeID) -> Result<BitVec, String> {
135        if !self.finished {
136            return Err("Assigment was not computed.".to_string());
137        }
138
139        // run a reachability analysis
140        let mut reachable = BitVec::with_capacity(self.residual_graph.number_of_nodes());
141        reachable.resize(self.residual_graph.number_of_nodes(), false);
142        let mut stack = vec![source];
143        stack.reserve(self.residual_graph.number_of_nodes());
144        reachable.set(source, true);
145        while let Some(node) = stack.pop() {
146            for edge in self.residual_graph.edge_range(node) {
147                let target = self.residual_graph.target(edge);
148                let reached = reachable.get(target).unwrap();
149                if !reached && self.residual_graph.data(edge).capacity > 0 {
150                    stack.push(target);
151                    reachable.set(target, true);
152                }
153            }
154        }
155        Ok(reachable)
156    }
157}
158
159#[cfg(test)]
160mod tests {
161
162    use std::sync::atomic::AtomicI32;
163    use std::sync::Arc;
164
165    use crate::edge::InputEdge;
166    use crate::ford_fulkerson::FordFulkerson;
167    use crate::max_flow::MaxFlow;
168    use crate::max_flow::ResidualEdgeData;
169    use bitvec::bits;
170    use bitvec::prelude::Lsb0;
171
172    #[test]
173    fn max_flow_clr() {
174        let edges = vec![
175            InputEdge::new(0, 1, ResidualEdgeData::new(16)),
176            InputEdge::new(0, 2, ResidualEdgeData::new(13)),
177            InputEdge::new(1, 2, ResidualEdgeData::new(10)),
178            InputEdge::new(1, 3, ResidualEdgeData::new(12)),
179            InputEdge::new(2, 1, ResidualEdgeData::new(4)),
180            InputEdge::new(2, 4, ResidualEdgeData::new(14)),
181            InputEdge::new(3, 2, ResidualEdgeData::new(9)),
182            InputEdge::new(3, 5, ResidualEdgeData::new(20)),
183            InputEdge::new(4, 3, ResidualEdgeData::new(7)),
184            InputEdge::new(4, 5, ResidualEdgeData::new(4)),
185        ];
186
187        let source = 0;
188        let target = 5;
189        let mut max_flow_solver = FordFulkerson::from_edge_list(edges, source, target);
190        max_flow_solver.run();
191
192        // it's OK to expect the solver to have run
193        let max_flow = max_flow_solver
194            .max_flow()
195            .expect("max flow computation did not run");
196        assert_eq!(23, max_flow);
197
198        // it's OK to expect the solver to have run
199        let assignment = max_flow_solver
200            .assignment(source)
201            .expect("assignment computation did not run");
202
203        assert_eq!(assignment, bits![1, 1, 1, 0, 1, 0]);
204    }
205
206    #[test]
207    fn run_with_upper_bound_no_effect() {
208        let edges = vec![
209            InputEdge::new(0, 1, ResidualEdgeData::new(16)),
210            InputEdge::new(0, 2, ResidualEdgeData::new(13)),
211            InputEdge::new(1, 2, ResidualEdgeData::new(10)),
212            InputEdge::new(1, 3, ResidualEdgeData::new(12)),
213            InputEdge::new(2, 1, ResidualEdgeData::new(4)),
214            InputEdge::new(2, 4, ResidualEdgeData::new(14)),
215            InputEdge::new(3, 2, ResidualEdgeData::new(9)),
216            InputEdge::new(3, 5, ResidualEdgeData::new(20)),
217            InputEdge::new(4, 3, ResidualEdgeData::new(7)),
218            InputEdge::new(4, 5, ResidualEdgeData::new(4)),
219        ];
220
221        let source = 0;
222        let target = 5;
223        let mut max_flow_solver = FordFulkerson::from_edge_list(edges, source, target);
224        let upper_bound = Arc::new(AtomicI32::new(1));
225        max_flow_solver.run_with_upper_bound(upper_bound);
226
227        // it's OK to expect the solver to have run
228        let max_flow = max_flow_solver
229            .max_flow()
230            .expect("max flow computation did not run");
231        assert_eq!(23, max_flow);
232    }
233
234    #[test]
235    fn max_flow_ita() {
236        let edges = vec![
237            InputEdge::new(0, 1, ResidualEdgeData::new(5)),
238            InputEdge::new(0, 4, ResidualEdgeData::new(7)),
239            InputEdge::new(0, 5, ResidualEdgeData::new(6)),
240            InputEdge::new(1, 2, ResidualEdgeData::new(4)),
241            InputEdge::new(1, 7, ResidualEdgeData::new(3)),
242            InputEdge::new(4, 7, ResidualEdgeData::new(4)),
243            InputEdge::new(4, 6, ResidualEdgeData::new(1)),
244            InputEdge::new(5, 6, ResidualEdgeData::new(5)),
245            InputEdge::new(2, 3, ResidualEdgeData::new(3)),
246            InputEdge::new(7, 3, ResidualEdgeData::new(7)),
247            InputEdge::new(6, 7, ResidualEdgeData::new(1)),
248            InputEdge::new(6, 3, ResidualEdgeData::new(6)),
249        ];
250
251        let source = 0;
252        let target = 3;
253        let mut max_flow_solver = FordFulkerson::from_edge_list(edges, source, target);
254        max_flow_solver.run();
255
256        // it's OK to expect the solver to have run
257        let max_flow = max_flow_solver
258            .max_flow()
259            .expect("max flow computation did not run");
260        assert_eq!(15, max_flow);
261
262        // it's OK to expect the solver to have run
263        let assignment = max_flow_solver
264            .assignment(source)
265            .expect("assignment computation did not run");
266        assert_eq!(assignment, bits![1, 0, 0, 0, 1, 1, 0, 0]);
267    }
268
269    #[test]
270    fn max_flow_yt() {
271        let edges = vec![
272            InputEdge::new(9, 0, ResidualEdgeData::new(5)),
273            InputEdge::new(9, 1, ResidualEdgeData::new(10)),
274            InputEdge::new(9, 2, ResidualEdgeData::new(15)),
275            InputEdge::new(0, 3, ResidualEdgeData::new(10)),
276            InputEdge::new(1, 0, ResidualEdgeData::new(15)),
277            InputEdge::new(1, 4, ResidualEdgeData::new(20)),
278            InputEdge::new(2, 5, ResidualEdgeData::new(25)),
279            InputEdge::new(3, 4, ResidualEdgeData::new(25)),
280            InputEdge::new(3, 6, ResidualEdgeData::new(10)),
281            InputEdge::new(4, 2, ResidualEdgeData::new(5)),
282            InputEdge::new(4, 7, ResidualEdgeData::new(30)),
283            InputEdge::new(5, 7, ResidualEdgeData::new(20)),
284            InputEdge::new(5, 8, ResidualEdgeData::new(10)),
285            InputEdge::new(7, 8, ResidualEdgeData::new(15)),
286            InputEdge::new(6, 10, ResidualEdgeData::new(5)),
287            InputEdge::new(7, 10, ResidualEdgeData::new(15)),
288            InputEdge::new(8, 10, ResidualEdgeData::new(10)),
289        ];
290
291        let source = 9;
292        let target = 10;
293        let mut max_flow_solver = FordFulkerson::from_edge_list(edges, source, target);
294        max_flow_solver.run();
295
296        // it's OK to expect the solver to have run
297        let max_flow = max_flow_solver
298            .max_flow()
299            .expect("max flow computation did not run");
300        assert_eq!(30, max_flow);
301
302        // it's OK to expect the solver to have run
303        let assignment = max_flow_solver
304            .assignment(source)
305            .expect("assignment computation did not run");
306        assert_eq!(assignment, bits![0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0]);
307    }
308
309    #[test]
310    fn max_flow_ff() {
311        let edges = vec![
312            InputEdge::new(0, 1, ResidualEdgeData::new(7)),
313            InputEdge::new(0, 2, ResidualEdgeData::new(3)),
314            InputEdge::new(1, 2, ResidualEdgeData::new(1)),
315            InputEdge::new(1, 3, ResidualEdgeData::new(6)),
316            InputEdge::new(2, 4, ResidualEdgeData::new(8)),
317            InputEdge::new(3, 5, ResidualEdgeData::new(2)),
318            InputEdge::new(3, 2, ResidualEdgeData::new(3)),
319            InputEdge::new(4, 3, ResidualEdgeData::new(2)),
320            InputEdge::new(4, 5, ResidualEdgeData::new(8)),
321        ];
322
323        let source = 0;
324        let target = 5;
325        let mut max_flow_solver = FordFulkerson::from_edge_list(edges, source, target);
326        max_flow_solver.run();
327
328        // it's OK to expect the solver to have run
329        let max_flow = max_flow_solver
330            .max_flow()
331            .expect("max flow computation did not run");
332        assert_eq!(9, max_flow);
333
334        // it's OK to expect the solver to have run
335        let assignment = max_flow_solver
336            .assignment(source)
337            .expect("assignment computation did not run");
338        assert_eq!(assignment, bits![1, 1, 0, 1, 0, 0]);
339    }
340
341    #[test]
342    #[should_panic]
343    fn flow_not_computed() {
344        let edges = vec![
345            InputEdge::new(0, 1, ResidualEdgeData::new(7)),
346            InputEdge::new(0, 2, ResidualEdgeData::new(3)),
347            InputEdge::new(1, 2, ResidualEdgeData::new(1)),
348            InputEdge::new(1, 3, ResidualEdgeData::new(6)),
349            InputEdge::new(2, 4, ResidualEdgeData::new(8)),
350            InputEdge::new(3, 5, ResidualEdgeData::new(2)),
351            InputEdge::new(3, 2, ResidualEdgeData::new(3)),
352            InputEdge::new(4, 3, ResidualEdgeData::new(2)),
353            InputEdge::new(4, 5, ResidualEdgeData::new(8)),
354        ];
355
356        // the expect(.) call is being tested
357        FordFulkerson::from_edge_list(edges, 0, 1)
358            .max_flow()
359            .expect("max flow computation did not run");
360    }
361
362    #[test]
363    #[should_panic]
364    fn assignment_not_computed() {
365        let edges = vec![
366            InputEdge::new(0, 1, ResidualEdgeData::new(7)),
367            InputEdge::new(0, 2, ResidualEdgeData::new(3)),
368            InputEdge::new(1, 2, ResidualEdgeData::new(1)),
369            InputEdge::new(1, 3, ResidualEdgeData::new(6)),
370            InputEdge::new(2, 4, ResidualEdgeData::new(8)),
371            InputEdge::new(3, 5, ResidualEdgeData::new(2)),
372            InputEdge::new(3, 2, ResidualEdgeData::new(3)),
373            InputEdge::new(4, 3, ResidualEdgeData::new(2)),
374            InputEdge::new(4, 5, ResidualEdgeData::new(8)),
375        ];
376
377        // the expect(.) call is being tested
378        FordFulkerson::from_edge_list(edges, 0, 1)
379            .assignment(0)
380            .expect("assignment computation did not run");
381    }
382}
toolbox_rs/ford_fulkerson.rs

toolbox_rs/
ford_fulkerson.rs
