toolbox_rs/
dinic.rs

1//! A Max-Flow computation implementing Cherkassky's variant of Dinitz' seminal
2//! algorithm. The implementation at hand is distinguished by three factors:
3//! 1) Computing the layer graph in a single BFS starting in t.
4//! 2) Omitting maintenance of the layer graph.
5//! 3) Running the augmentation phase as a single DFS.
6//!
7//! The DFS restarts after it found an augmenting path on the tail of the
8//! saturated edge that is closest to the source.
9use crate::{
10    edge::InputEdge,
11    graph::{Graph, NodeID},
12    max_flow::{MaxFlow, ResidualEdgeData},
13    static_graph::StaticGraph,
14};
15use bitvec::vec::BitVec;
16use core::cmp::min;
17use log::debug;
18use std::{
19    collections::VecDeque,
20    sync::{
21        atomic::{AtomicI32, Ordering},
22        Arc,
23    },
24    time::Instant,
25};
26
27pub struct Dinic {
28    residual_graph: StaticGraph<ResidualEdgeData>,
29    max_flow: i32,
30    finished: bool,
31    level: Vec<usize>,
32    parents: Vec<NodeID>,
33    stack: Vec<(NodeID, i32)>,
34    dfs_count: usize,
35    bfs_count: usize,
36    queue: VecDeque<NodeID>,
37    source: NodeID,
38    target: NodeID,
39    bound: Option<Arc<AtomicI32>>,
40}
41
42impl Dinic {
43    fn bfs(&mut self) -> bool {
44        let start = Instant::now();
45        self.bfs_count += 1;
46        // init
47        self.level.fill(usize::MAX);
48        self.level[self.target] = 0;
49
50        self.queue.clear();
51        self.queue.push_back(self.target);
52
53        let duration = start.elapsed();
54        debug!("BFS init: {:?}", duration);
55
56        // label residual graph nodes in BFS order, but in reverse starting from the target
57        while let Some(u) = self.queue.pop_front() {
58            for edge in self.residual_graph.edge_range(u) {
59                let v = self.residual_graph.target(edge);
60                if v != self.source && self.level[v] != usize::MAX {
61                    // node v is not the source, and is already visited. Note the source can be reached multiple times
62                    continue;
63                }
64
65                // check capacity of reverse edge
66                let rev_edge = self.residual_graph.find_edge_unchecked(v, u);
67                let edge_capacity = self.residual_graph.data(rev_edge).capacity;
68                if edge_capacity < 1 {
69                    // no capacity to use on this edge
70                    continue;
71                }
72                self.level[v] = self.level[u] + 1;
73                if v != self.source {
74                    self.queue.push_back(v);
75                }
76            }
77        }
78        let duration = start.elapsed();
79        debug!(
80            "BFS took: {:?}, upper bound on path length: {}",
81            duration, self.level[self.source]
82        );
83        self.level[self.source] != usize::MAX
84    }
85
86    fn dfs(&mut self) -> i32 {
87        let start = Instant::now();
88        self.dfs_count += 1;
89        self.stack.clear();
90        self.stack.push((self.source, i32::MAX));
91
92        self.parents.fill(NodeID::MAX);
93        self.parents[self.source] = self.source;
94
95        let duration = start.elapsed();
96        debug!(" DFS init: {:?}", duration);
97        let mut blocking_flow = 0;
98        while let Some((u, flow)) = self.stack.pop() {
99            for edge in self.residual_graph.edge_range(u) {
100                let v = self.residual_graph.target(edge);
101                if self.parents[v] != NodeID::MAX {
102                    // v already in queue
103                    continue;
104                }
105                if self.level[u] < self.level[v] {
106                    // edge is not leading to target on a path in the BFS tree
107                    continue;
108                }
109                let available_capacity = self.residual_graph.data(edge).capacity;
110                if available_capacity == 0 {
111                    // no capacity to use on this edge
112                    continue;
113                }
114                self.parents[v] = u;
115                let flow = min(flow, available_capacity);
116                if v == self.target {
117                    let duration = start.elapsed();
118                    debug!(" reached target {}: {:?}", v, duration);
119                    // reached a target. Unpack path in reverse order, assign flow
120                    let mut v = v; // mutable shadow
121                    let mut closest_tail = u;
122                    loop {
123                        let u = self.parents[v];
124                        if u == v {
125                            break;
126                        }
127                        let fwd_edge = self.residual_graph.find_edge_unchecked(u, v);
128                        self.residual_graph.data_mut(fwd_edge).capacity -= flow;
129                        if 0 == self.residual_graph.data_mut(fwd_edge).capacity {
130                            closest_tail = u;
131                        }
132                        let rev_edge = self.residual_graph.find_edge_unchecked(v, u);
133                        self.residual_graph.data_mut(rev_edge).capacity += flow;
134                        v = u;
135                    }
136                    let duration = start.elapsed();
137                    debug!(" augmentation took: {:?}", duration);
138
139                    // unwind stack till tail node, then continue the search
140                    let before = self.stack.len();
141                    while let Some((node, _)) = self.stack.pop() {
142                        if self.parents[node] == closest_tail {
143                            break; // while let
144                        }
145                    }
146                    blocking_flow += flow;
147                    debug!(" stack len before: {before}, after: {}", self.stack.len());
148
149                    // make target reachable again
150                    self.parents[self.target] = NodeID::MAX;
151                    self.dfs_count += 1;
152
153                    break; // for edge
154                } else {
155                    self.stack.push((v, flow));
156                }
157            }
158        }
159
160        let duration = start.elapsed();
161        debug!("DFS took: {:?} (unsuccessful)", duration);
162        blocking_flow
163    }
164}
165
166impl MaxFlow for Dinic {
167    fn from_edge_list(
168        mut edge_list: Vec<InputEdge<ResidualEdgeData>>,
169        source: usize,
170        target: usize,
171    ) -> Self {
172        debug_assert!(!edge_list.is_empty());
173        let number_of_edges = edge_list.len();
174
175        debug!("extending {} edges", edge_list.len());
176        // blindly generate reverse edges for all edges with zero capacity
177        edge_list.extend_from_within(..);
178        debug!("into {} edges", edge_list.len());
179
180        edge_list.iter_mut().skip(number_of_edges).for_each(|edge| {
181            edge.reverse();
182            edge.data.capacity = 0;
183        });
184        debug!("sorting after reversing");
185
186        // dedup-merge edge set, by using the following trick: not the dedup(.) call
187        // below takes the second argument as mut. When deduping equivalent values
188        // a and b, then a is accumulated onto b.
189        edge_list.sort_unstable_by(|a, b| {
190            if a.source == b.source {
191                return a.target.cmp(&b.target);
192            }
193            a.source.cmp(&b.source)
194        });
195        debug!("start dedup");
196        edge_list.dedup_by(|a, b| {
197            // edges a and b are assumed to be equivalent in the residual graph if
198            // (and only if) they are parallel. In other words, this removes parallel
199            // edges in the residual graph and accumulates capacities on the remaining
200            // egde.
201            let edges_are_parallel = a.is_parallel_to(b);
202            if edges_are_parallel {
203                b.data.capacity += a.data.capacity;
204            }
205            edges_are_parallel
206        });
207        edge_list.shrink_to_fit();
208        debug!("dedup done");
209
210        // at this point the edge set of the residual graph doesn't have any
211        // duplicates anymore. Note that this is fine, as we are looking to
212        // compute a node partition.
213
214        let residual_graph = StaticGraph::new_from_sorted_list(edge_list);
215        let number_of_nodes = residual_graph.number_of_nodes();
216
217        Self {
218            residual_graph,
219            max_flow: 0,
220            finished: false,
221            level: Vec::with_capacity(number_of_nodes),
222            parents: Vec::with_capacity(number_of_nodes),
223            stack: Vec::with_capacity(number_of_nodes),
224            dfs_count: 0,
225            bfs_count: 0,
226            queue: VecDeque::with_capacity(number_of_nodes),
227            source,
228            target,
229            bound: None,
230        }
231    }
232
233    fn run_with_upper_bound(&mut self, bound: Arc<AtomicI32>) {
234        debug!("upper bound: {}", bound.load(Ordering::Relaxed));
235
236        self.bound = Some(bound);
237        self.run()
238    }
239
240    fn run(&mut self) {
241        debug!(
242            "residual graph size: V {}, E {}",
243            self.residual_graph.number_of_nodes(),
244            self.residual_graph.number_of_edges()
245        );
246
247        let number_of_nodes = self.residual_graph.number_of_nodes();
248        self.parents.resize(number_of_nodes, 0);
249        self.level.resize(number_of_nodes, usize::MAX);
250        self.queue.reserve(number_of_nodes);
251
252        let mut flow = 0;
253        while self.bfs() {
254            flow += self.dfs();
255            if let Some(bound) = &self.bound {
256                // break early if an upper bound is known to the computation
257                if flow > bound.load(Ordering::Relaxed) {
258                    debug!("aborting max flow computation at {flow}");
259                    self.max_flow = flow;
260                    return;
261                }
262            }
263        }
264        if let Some(bound) = &self.bound {
265            bound.fetch_min(flow, Ordering::Relaxed);
266        }
267        self.max_flow = flow;
268        self.finished = true;
269    }
270
271    fn max_flow(&self) -> Result<i32, String> {
272        if !self.finished {
273            return Err("Assigment was not computed.".to_string());
274        }
275        debug!(
276            "finished in {} DFS, and {} BFS runs",
277            self.dfs_count, self.bfs_count
278        );
279        Ok(self.max_flow)
280    }
281
282    fn assignment(&self, source: NodeID) -> Result<BitVec, String> {
283        if !self.finished {
284            return Err("Assigment was not computed.".to_string());
285        }
286
287        // run a reachability analysis
288        let mut reachable = BitVec::new();
289        reachable.resize(self.residual_graph.number_of_nodes(), false);
290        let mut stack = vec![source];
291        stack.reserve(self.residual_graph.number_of_nodes());
292        reachable.set(source, true);
293        while let Some(node) = stack.pop() {
294            for edge in self.residual_graph.edge_range(node) {
295                let target = self.residual_graph.target(edge);
296                let reached = reachable.get(target).unwrap();
297                if !reached && self.residual_graph.data(edge).capacity > 0 {
298                    stack.push(target);
299                    reachable.set(target, true);
300                }
301            }
302        }
303        Ok(reachable)
304    }
305}
306
307#[cfg(test)]
308mod tests {
309
310    use crate::dinic::Dinic;
311    use crate::edge::EdgeData;
312    use crate::edge::InputEdge;
313    use crate::max_flow::MaxFlow;
314    use crate::max_flow::ResidualEdgeData;
315    use bitvec::bits;
316    use bitvec::prelude::Lsb0;
317
318    #[test]
319    fn max_flow_clr() {
320        let edges = vec![
321            InputEdge::new(0, 1, ResidualEdgeData::new(16)),
322            InputEdge::new(0, 2, ResidualEdgeData::new(13)),
323            InputEdge::new(1, 2, ResidualEdgeData::new(10)),
324            InputEdge::new(1, 3, ResidualEdgeData::new(12)),
325            InputEdge::new(2, 1, ResidualEdgeData::new(4)),
326            InputEdge::new(2, 4, ResidualEdgeData::new(14)),
327            InputEdge::new(3, 2, ResidualEdgeData::new(9)),
328            InputEdge::new(3, 5, ResidualEdgeData::new(20)),
329            InputEdge::new(4, 3, ResidualEdgeData::new(7)),
330            InputEdge::new(4, 5, ResidualEdgeData::new(4)),
331        ];
332
333        let source = 0;
334        let target = 5;
335        let mut max_flow_solver = Dinic::from_edge_list(edges, source, target);
336        max_flow_solver.run();
337
338        // it's OK to expect the solver to have run
339        let max_flow = max_flow_solver
340            .max_flow()
341            .expect("max flow computation did not run");
342        assert_eq!(23, max_flow);
343
344        // it's OK to expect the solver to have run
345        let assignment = max_flow_solver
346            .assignment(source)
347            .expect("assignment computation did not run");
348
349        assert_eq!(assignment, bits![1, 1, 1, 0, 1, 0]);
350    }
351
352    #[test]
353    fn max_flow_ita_from_generic_edge_list() {
354        let edges = vec![
355            InputEdge::new(0, 1, 5),
356            InputEdge::new(0, 4, 7),
357            InputEdge::new(0, 5, 6),
358            InputEdge::new(1, 2, 4),
359            InputEdge::new(1, 7, 3),
360            InputEdge::new(4, 7, 4),
361            InputEdge::new(4, 6, 1),
362            InputEdge::new(5, 6, 5),
363            InputEdge::new(2, 3, 3),
364            InputEdge::new(7, 3, 7),
365            InputEdge::new(6, 7, 1),
366            InputEdge::new(6, 3, 6),
367        ];
368
369        let source = 0;
370        let target = 3;
371        let mut max_flow_solver = Dinic::from_generic_edge_list(&edges, source, target, |edge| {
372            ResidualEdgeData::new(*edge.data())
373        });
374        max_flow_solver.run();
375
376        // it's OK to expect the solver to have run
377        let max_flow = max_flow_solver
378            .max_flow()
379            .expect("max flow computation did not run");
380        assert_eq!(15, max_flow);
381
382        // it's OK to expect the solver to have run
383        let assignment = max_flow_solver
384            .assignment(source)
385            .expect("assignment computation did not run");
386        assert_eq!(assignment, bits![1, 0, 0, 0, 1, 1, 0, 0]);
387    }
388
389    #[test]
390    fn max_flow_ita() {
391        let edges = vec![
392            InputEdge::new(0, 1, ResidualEdgeData::new(5)),
393            InputEdge::new(0, 4, ResidualEdgeData::new(7)),
394            InputEdge::new(0, 5, ResidualEdgeData::new(6)),
395            InputEdge::new(1, 2, ResidualEdgeData::new(4)),
396            InputEdge::new(1, 7, ResidualEdgeData::new(3)),
397            InputEdge::new(4, 7, ResidualEdgeData::new(4)),
398            InputEdge::new(4, 6, ResidualEdgeData::new(1)),
399            InputEdge::new(5, 6, ResidualEdgeData::new(5)),
400            InputEdge::new(2, 3, ResidualEdgeData::new(3)),
401            InputEdge::new(7, 3, ResidualEdgeData::new(7)),
402            InputEdge::new(6, 7, ResidualEdgeData::new(1)),
403            InputEdge::new(6, 3, ResidualEdgeData::new(6)),
404        ];
405
406        let source = 0;
407        let target = 3;
408        let mut max_flow_solver = Dinic::from_edge_list(edges, source, target);
409        max_flow_solver.run();
410
411        // it's OK to expect the solver to have run
412        let max_flow = max_flow_solver
413            .max_flow()
414            .expect("max flow computation did not run");
415        assert_eq!(15, max_flow);
416
417        // it's OK to expect the solver to have run
418        let assignment = max_flow_solver
419            .assignment(source)
420            .expect("assignment computation did not run");
421        assert_eq!(assignment, bits![1, 0, 0, 0, 1, 1, 0, 0]);
422    }
423
424    #[test]
425    fn max_flow_yt() {
426        let edges = vec![
427            InputEdge::new(9, 0, ResidualEdgeData::new(5)),
428            InputEdge::new(9, 1, ResidualEdgeData::new(10)),
429            InputEdge::new(9, 2, ResidualEdgeData::new(15)),
430            InputEdge::new(0, 3, ResidualEdgeData::new(10)),
431            InputEdge::new(1, 0, ResidualEdgeData::new(15)),
432            InputEdge::new(1, 4, ResidualEdgeData::new(20)),
433            InputEdge::new(2, 5, ResidualEdgeData::new(25)),
434            InputEdge::new(3, 4, ResidualEdgeData::new(25)),
435            InputEdge::new(3, 6, ResidualEdgeData::new(10)),
436            InputEdge::new(4, 2, ResidualEdgeData::new(5)),
437            InputEdge::new(4, 7, ResidualEdgeData::new(30)),
438            InputEdge::new(5, 7, ResidualEdgeData::new(20)),
439            InputEdge::new(5, 8, ResidualEdgeData::new(10)),
440            InputEdge::new(7, 8, ResidualEdgeData::new(15)),
441            InputEdge::new(6, 10, ResidualEdgeData::new(5)),
442            InputEdge::new(7, 10, ResidualEdgeData::new(15)),
443            InputEdge::new(8, 10, ResidualEdgeData::new(10)),
444        ];
445
446        let source = 9;
447        let target = 10;
448        let mut max_flow_solver = Dinic::from_edge_list(edges, source, target);
449        max_flow_solver.run();
450
451        // it's OK to expect the solver to have run
452        let max_flow = max_flow_solver
453            .max_flow()
454            .expect("max flow computation did not run");
455        assert_eq!(30, max_flow);
456
457        // it's OK to expect the solver to have run
458        let assignment = max_flow_solver
459            .assignment(source)
460            .expect("assignment computation did not run");
461        assert_eq!(assignment, bits![0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0]);
462    }
463
464    #[test]
465    fn max_flow_ff() {
466        let edges = vec![
467            InputEdge::new(0, 1, ResidualEdgeData::new(7)),
468            InputEdge::new(0, 2, ResidualEdgeData::new(3)),
469            InputEdge::new(1, 2, ResidualEdgeData::new(1)),
470            InputEdge::new(1, 3, ResidualEdgeData::new(6)),
471            InputEdge::new(2, 4, ResidualEdgeData::new(8)),
472            InputEdge::new(3, 5, ResidualEdgeData::new(2)),
473            InputEdge::new(3, 2, ResidualEdgeData::new(3)),
474            InputEdge::new(4, 3, ResidualEdgeData::new(2)),
475            InputEdge::new(4, 5, ResidualEdgeData::new(8)),
476        ];
477
478        let source = 0;
479        let target = 5;
480        let mut max_flow_solver = Dinic::from_edge_list(edges, source, target);
481        max_flow_solver.run();
482
483        // it's OK to expect the solver to have run
484        let max_flow = max_flow_solver
485            .max_flow()
486            .expect("max flow computation did not run");
487        assert_eq!(9, max_flow);
488
489        // it's OK to expect the solver to have run
490        let assignment = max_flow_solver
491            .assignment(source)
492            .expect("assignment computation did not run");
493        assert_eq!(assignment, bits![1, 1, 0, 1, 0, 0]);
494    }
495
496    #[test]
497    #[should_panic]
498    fn flow_not_computed() {
499        let edges = vec![
500            InputEdge::new(0, 1, ResidualEdgeData::new(7)),
501            InputEdge::new(0, 2, ResidualEdgeData::new(3)),
502            InputEdge::new(1, 2, ResidualEdgeData::new(1)),
503            InputEdge::new(1, 3, ResidualEdgeData::new(6)),
504            InputEdge::new(2, 4, ResidualEdgeData::new(8)),
505            InputEdge::new(3, 5, ResidualEdgeData::new(2)),
506            InputEdge::new(3, 2, ResidualEdgeData::new(3)),
507            InputEdge::new(4, 3, ResidualEdgeData::new(2)),
508            InputEdge::new(4, 5, ResidualEdgeData::new(8)),
509        ];
510
511        // the expect(.) call is being tested
512        Dinic::from_edge_list(edges, 1, 2)
513            .max_flow()
514            .expect("max flow computation did not run");
515    }
516
517    #[test]
518    #[should_panic]
519    fn assignment_not_computed() {
520        let edges = vec![
521            InputEdge::new(0, 1, ResidualEdgeData::new(7)),
522            InputEdge::new(0, 2, ResidualEdgeData::new(3)),
523            InputEdge::new(1, 2, ResidualEdgeData::new(1)),
524            InputEdge::new(1, 3, ResidualEdgeData::new(6)),
525            InputEdge::new(2, 4, ResidualEdgeData::new(8)),
526            InputEdge::new(3, 5, ResidualEdgeData::new(2)),
527            InputEdge::new(3, 2, ResidualEdgeData::new(3)),
528            InputEdge::new(4, 3, ResidualEdgeData::new(2)),
529            InputEdge::new(4, 5, ResidualEdgeData::new(8)),
530        ];
531
532        // the expect(.) call is being tested
533        Dinic::from_edge_list(edges, 1, 2)
534            .assignment(1)
535            .expect("assignment computation did not run");
536    }
537}
toolbox_rs/dinic.rs

toolbox_rs/
dinic.rs
