oxicuda_graphalg/max_flow/
edmonds_karp.rs

1//! Edmonds-Karp max flow: BFS augmenting paths. O(VE^2).
2
3use std::collections::VecDeque;
4
5use crate::error::{GraphalgError, GraphalgResult};
6
7/// Residual graph as an n x n capacity matrix (row-major).
8#[derive(Debug, Clone)]
9pub struct FlowNetwork {
10    pub n: usize,
11    pub cap: Vec<f64>,
12}
13
14impl FlowNetwork {
15    pub fn new(n: usize) -> Self {
16        Self {
17            n,
18            cap: vec![0.0; n * n],
19        }
20    }
21
22    pub fn add_edge(&mut self, u: usize, v: usize, c: f64) -> GraphalgResult<()> {
23        if u >= self.n || v >= self.n {
24            return Err(GraphalgError::IndexOutOfBounds {
25                index: u.max(v),
26                len: self.n,
27            });
28        }
29        if c < 0.0 {
30            return Err(GraphalgError::InvalidEdgeWeight(format!(
31                "negative capacity {c} on ({u},{v})"
32            )));
33        }
34        self.cap[u * self.n + v] += c;
35        Ok(())
36    }
37
38    pub fn c(&self, u: usize, v: usize) -> f64 {
39        self.cap[u * self.n + v]
40    }
41
42    pub fn c_mut(&mut self, u: usize, v: usize) -> &mut f64 {
43        &mut self.cap[u * self.n + v]
44    }
45}
46
47pub fn edmonds_karp(net: &FlowNetwork, source: usize, sink: usize) -> GraphalgResult<f64> {
48    if source >= net.n || sink >= net.n {
49        return Err(GraphalgError::SourceOutOfRange {
50            node: source.max(sink),
51            n: net.n,
52        });
53    }
54    let mut residual = net.clone();
55    let mut max_flow = 0.0;
56    loop {
57        // BFS to find augmenting path
58        let mut parent = vec![usize::MAX; net.n];
59        parent[source] = source;
60        let mut q: VecDeque<usize> = VecDeque::new();
61        q.push_back(source);
62        while let Some(u) = q.pop_front() {
63            if u == sink {
64                break;
65            }
66            for v in 0..net.n {
67                if parent[v] == usize::MAX && residual.c(u, v) > 1e-12 {
68                    parent[v] = u;
69                    q.push_back(v);
70                }
71            }
72        }
73        if parent[sink] == usize::MAX {
74            break;
75        }
76        // Find bottleneck
77        let mut bottleneck = f64::INFINITY;
78        let mut v = sink;
79        while v != source {
80            let u = parent[v];
81            bottleneck = bottleneck.min(residual.c(u, v));
82            v = u;
83        }
84        // Augment
85        let mut v = sink;
86        while v != source {
87            let u = parent[v];
88            *residual.c_mut(u, v) -= bottleneck;
89            *residual.c_mut(v, u) += bottleneck;
90            v = u;
91        }
92        max_flow += bottleneck;
93    }
94    Ok(max_flow)
95}
96
97#[cfg(test)]
98mod tests {
99    use super::*;
100
101    #[test]
102    fn ek_simple() {
103        // s -> a (3), s -> b (2), a -> t (3), b -> t (2)
104        let mut net = FlowNetwork::new(4);
105        net.add_edge(0, 1, 3.0).expect("ok");
106        net.add_edge(0, 2, 2.0).expect("ok");
107        net.add_edge(1, 3, 3.0).expect("ok");
108        net.add_edge(2, 3, 2.0).expect("ok");
109        let f = edmonds_karp(&net, 0, 3).expect("ok");
110        assert!((f - 5.0).abs() < 1e-9);
111    }
112
113    #[test]
114    fn ek_bottleneck() {
115        let mut net = FlowNetwork::new(3);
116        net.add_edge(0, 1, 10.0).expect("ok");
117        net.add_edge(1, 2, 1.0).expect("ok");
118        let f = edmonds_karp(&net, 0, 2).expect("ok");
119        assert!((f - 1.0).abs() < 1e-9);
120    }
121}
oxicuda_graphalg/max_flow/edmonds_karp.rs

oxicuda_graphalg/max_flow/
edmonds_karp.rs
