rs_graph/maxflow/edmondskarp.rs
1/*
2 * Copyright (c) 2017-2022 Frank Fischer <frank-fischer@shadow-soft.de>
3 *
4 * This program is free software: you can redistribute it and/or
5 * modify it under the terms of the GNU General Public License as
6 * published by the Free Software Foundation, either version 3 of the
7 * License, or (at your option) any later version.
8 *
9 * This program is distributed in the hope that it will be useful, but
10 * WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 * General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with this program. If not, see <http://www.gnu.org/licenses/>
16 */
17
18//! This module implements the max flow algorithm of Edmonds-Karp.
19//!
20//! # Example
21//!
22//! ```
23//! use rs_graph::traits::*;
24//! use rs_graph::maxflow::edmondskarp;
25//! use rs_graph::Net;
26//! use rs_graph::string::{Data, from_ascii};
27//!
28//! let Data { graph: g, weights: upper, nodes } = from_ascii::<Net>(r"
29//! a---2-->b
30//! @|\ ^\
31//! / | \ | 4
32//! 5 | \ | \
33//! / | | | @
34//! s 1 1 2 t
35//! \ | | | @
36//! 5 | \ | /
37//! \ | \ | 5
38//! @v @|/
39//! c---2-->d
40//! ").unwrap();
41//!
42//! let s = g.id2node(nodes[&'s']);
43//! let t = g.id2node(nodes[&'t']);
44//! let v1 = g.id2node(nodes[&'a']);
45//! let v2 = g.id2node(nodes[&'b']);
46//! let v3 = g.id2node(nodes[&'c']);
47//! let v4 = g.id2node(nodes[&'d']);
48//!
49//! let (value, flow, mut mincut) = edmondskarp(&g, s, t, |e| upper[e.index()]);
50//!
51//! assert_eq!(value, 5);
52//! assert!(flow.iter().all(|&(e, f)| f >= 0 && f <= upper[e.index()]));
53//! assert!(g.nodes().filter(|&u| u != s && u != t).all(|u| {
54//! g.outedges(u).map(|(e,_)| flow[g.edge_id(e)].1).sum::<usize>() ==
55//! g.inedges(u).map(|(e,_)| flow[g.edge_id(e)].1).sum::<usize>()
56//! }));
57//!
58//! mincut.sort_by_key(|u| u.index());
59//! assert_eq!(mincut, vec![v1, s, v3]);
60//! ```
61//!
62//! ```
63//! use rs_graph::traits::*;
64//! use rs_graph::maxflow::edmondskarp;
65//! use rs_graph::Net;
66//! use rs_graph::string::{Data, from_ascii};
67//!
68//! let Data { graph: g, weights: upper, nodes } = from_ascii::<Net>(r"
69//! ---8-->a---10---
70//! / | \
71//! / 1 --3-- |
72//! / | / \ |
73//! / v@ \v
74//! ---->b-----9---->c----8--->d----
75//! / \ @ @^ \
76//! 18 ---6-- / / | 33
77//! / \/ / | \
78//! / /\ ----- | @
79//! s --5--- | / | t
80//! \ / | / | @
81//! 27 | ----2-|-- / /
82//! \ | / | /----8--- 6
83//! \ |/ @ | /
84//! ---->e----9-->f------6---->g----
85//! \ | @
86//! \ | /
87//! --5----->h---4---
88//! ").unwrap();
89//!
90//! let s = g.id2node(nodes[&'s']);
91//! let t = g.id2node(nodes[&'t']);
92//!
93//! assert_eq!(g.num_edges(), 18);
94//!
95//! let (value, flow, mut mincut) = edmondskarp(&g, s, t, |e| upper[e.index()]);
96//! assert_eq!(value, 29);
97//!
98//! mincut.sort_by_key(|u| u.index());
99//! assert_eq!(mincut, "bcsef".chars().map(|v| g.id2node(nodes[&v])).collect::<Vec<_>>());
100//! ```
101
102use crate::traits::IndexDigraph;
103
104use std::cmp::min;
105use std::collections::VecDeque;
106
107use crate::num::traits::NumAssign;
108
109/// Max-flow algorithm of Edmonds and Karp.
110pub struct EdmondsKarp<'a, G, F>
111where
112 G: IndexDigraph,
113{
114 g: &'a G,
115 neighs: Vec<Vec<(usize, usize)>>,
116 pred: Vec<(usize, usize)>,
117 flow: Vec<F>,
118 queue: VecDeque<usize>,
119 value: F,
120}
121
122impl<'a, G, F> EdmondsKarp<'a, G, F>
123where
124 G: IndexDigraph,
125 F: NumAssign + Ord + Copy,
126{
127 /// Create a new Dinic algorithm instance for a graph.
128 pub fn new(g: &'a G) -> Self {
129 EdmondsKarp {
130 g,
131 neighs: g
132 .nodes()
133 .map(|u| {
134 g.outedges(u)
135 .map(|(e, v)| (g.edge_id(e) << 1, g.node_id(v)))
136 .chain(g.inedges(u).map(|(e, v)| ((g.edge_id(e) << 1) | 1, g.node_id(v))))
137 .collect()
138 })
139 .collect(),
140 pred: vec![(usize::max_value(), usize::max_value()); g.num_nodes()],
141 flow: vec![F::zero(); g.num_edges() * 2],
142 queue: VecDeque::with_capacity(g.num_nodes()),
143 value: F::zero(),
144 }
145 }
146
147 /// Return the underlying graph.
148 pub fn as_graph(&self) -> &'a G {
149 self.g
150 }
151
152 /// Return the value of the latest computed maximum flow.
153 pub fn value(&self) -> F {
154 self.value
155 }
156
157 /// Return the flow value on edge `e`
158 pub fn flow(&self, e: G::Edge<'_>) -> F {
159 self.flow[self.g.edge_id(e) << 1]
160 }
161
162 /// Return an iterator over all (edge, flow) pairs.
163 pub fn flow_iter<'b>(&'b self) -> impl Iterator<Item = (G::Edge<'a>, F)> + 'b {
164 self.g.edges().enumerate().map(move |(i, e)| (e, self.flow[i << 1]))
165 }
166
167 pub fn solve<Us>(&mut self, src: G::Node<'_>, snk: G::Node<'_>, upper: Us)
168 where
169 Us: Fn(G::Edge<'a>) -> F,
170 {
171 let src = self.g.node_id(src);
172 let snk = self.g.node_id(snk);
173 assert_ne!(src, snk, "Source and sink node must not be equal");
174
175 // initialize network flow
176 for (e, flw) in self.flow.iter_mut().enumerate() {
177 *flw = if (e & 1) == 0 {
178 F::zero()
179 } else {
180 upper(self.g.id2edge(e >> 1))
181 };
182 }
183 self.value = F::zero();
184
185 // nothing to do if there is no edge
186 if self.g.num_edges() == 0 {
187 return;
188 }
189
190 loop {
191 // do bfs from source to sink
192 self.pred.fill((usize::max_value(), usize::max_value()));
193
194 // just some dummy edge
195 self.pred[src] = (0, 0);
196 self.queue.clear();
197 self.queue.push_back(src);
198 'bfs: while let Some(u) = self.queue.pop_front() {
199 for &(e, v) in &self.neighs[u] {
200 if self.pred[v].0 == usize::max_value() && !self.flow[e ^ 1].is_zero() {
201 self.pred[v] = (e, u);
202 self.queue.push_back(v);
203 if v == snk {
204 break 'bfs;
205 }
206 }
207 }
208 }
209
210 // sink cannot be reached -> stop
211 if self.pred[snk].0 == usize::max_value() {
212 break;
213 }
214
215 // compute augmentation value
216 let mut v = snk;
217 let mut df = self.flow[self.pred[v].0 ^ 1];
218 while v != src {
219 let (e, u) = self.pred[v];
220 df = min(df, self.flow[e ^ 1]);
221 v = u;
222 }
223
224 debug_assert!(!df.is_zero());
225
226 // now augment the flow
227 let mut v = snk;
228 while v != src {
229 let (e, u) = self.pred[v];
230 self.flow[e] += df;
231 self.flow[e ^ 1] -= df;
232 v = u;
233 }
234
235 self.value += df;
236 }
237 }
238
239 /// Return the minimal cut associated with the last maximum flow.
240 pub fn mincut(&self) -> Vec<G::Node<'a>> {
241 self.g
242 .nodes()
243 .filter(|&u| self.pred[self.g.node_id(u)].0 != usize::max_value())
244 .collect()
245 }
246}
247
248/// Solve the maxflow problem using the algorithm of Edmonds-Karp.
249///
250/// The function solves the max flow problem from the source nodes
251/// `src` to the sink node `snk` with the given `upper` bounds on
252/// the edges.
253///
254/// The function returns the flow value, the flow on each edge and the
255/// nodes in a minimal cut.
256pub fn edmondskarp<'a, G, F, Us>(
257 g: &'a G,
258 src: G::Node<'a>,
259 snk: G::Node<'a>,
260 upper: Us,
261) -> (F, Vec<(G::Edge<'a>, F)>, Vec<G::Node<'a>>)
262where
263 G: IndexDigraph,
264 F: 'a + NumAssign + Ord + Copy,
265 Us: Fn(G::Edge<'_>) -> F,
266{
267 let mut maxflow = EdmondsKarp::new(g);
268 maxflow.solve(src, snk, upper);
269 (maxflow.value(), maxflow.flow_iter().collect(), maxflow.mincut())
270}