competitive/graph/
minimum_cost_flow.rs1use super::{BidirectionalSparseGraph, Graph};
2use std::{cmp::Reverse, collections::BinaryHeap};
3
4#[derive(Debug, Clone)]
5pub struct PrimalDualBuilder {
6 vsize: usize,
7 edges: Vec<(usize, usize)>,
8 capacities: Vec<u64>,
9 costs: Vec<i64>,
10 has_negedge: bool,
11}
12impl PrimalDualBuilder {
13 pub fn new(vsize: usize, esize_expect: usize) -> Self {
14 Self {
15 vsize,
16 edges: Vec::with_capacity(esize_expect),
17 capacities: Vec::with_capacity(esize_expect * 2),
18 costs: Vec::with_capacity(esize_expect * 2),
19 has_negedge: false,
20 }
21 }
22 pub fn add_edge(&mut self, from: usize, to: usize, cap: u64, cost: i64) {
23 self.edges.push((from, to));
24 self.capacities.push(cap);
25 self.capacities.push(0);
26 self.has_negedge |= cost < 0;
27 self.costs.push(cost);
28 self.costs.push(-cost);
29 }
30 pub fn gen_graph(&mut self) -> BidirectionalSparseGraph {
31 let edges = std::mem::take(&mut self.edges);
32 BidirectionalSparseGraph::from_edges(self.vsize, edges)
33 }
34 pub fn build(self, graph: &BidirectionalSparseGraph) -> PrimalDual<'_> {
35 let PrimalDualBuilder {
36 vsize,
37 capacities,
38 costs,
39 has_negedge,
40 ..
41 } = self;
42 PrimalDual {
43 graph,
44 capacities,
45 costs,
46 potential: std::iter::repeat_n(0, vsize).collect(),
47 dist: Vec::with_capacity(vsize),
48 prev_vertex: std::iter::repeat_n(0, vsize).collect(),
49 prev_edge: std::iter::repeat_n(0, vsize).collect(),
50 has_negedge,
51 queue: BinaryHeap::new(),
52 }
53 }
54}
55impl Extend<(usize, usize, u64, i64)> for PrimalDualBuilder {
56 fn extend<T: IntoIterator<Item = (usize, usize, u64, i64)>>(&mut self, iter: T) {
57 for (from, to, cap, cost) in iter {
58 self.add_edge(from, to, cap, cost)
59 }
60 }
61}
62
63#[derive(Debug)]
64pub struct PrimalDual<'a> {
65 graph: &'a BidirectionalSparseGraph,
66 capacities: Vec<u64>,
67 costs: Vec<i64>,
68 potential: Vec<i64>,
69 dist: Vec<i64>,
70 prev_vertex: Vec<usize>,
71 prev_edge: Vec<usize>,
72 has_negedge: bool,
73 queue: BinaryHeap<(Reverse<i64>, usize)>,
74}
75impl PrimalDual<'_> {
76 pub fn builder(vsize: usize, esize_expect: usize) -> PrimalDualBuilder {
77 PrimalDualBuilder::new(vsize, esize_expect)
78 }
79 fn bellman_ford(&mut self, s: usize) {
80 self.potential.clear();
81 self.potential.resize(self.graph.vertices_size(), i64::MAX);
82 self.potential[s] = 0;
83 for _ in 1..self.graph.vertices_size() {
84 let mut end = true;
85 for u in self.graph.vertices() {
86 if self.potential[u] == i64::MAX {
87 continue;
88 }
89 for a in self.graph.neighbors(u) {
90 if self.capacities[a.label] == 0 {
91 continue;
92 }
93 let ncost = self.potential[u].saturating_add(self.costs[a.label]);
94 if self.potential[a.to] > ncost {
95 self.potential[a.to] = ncost;
96 end = false;
97 }
98 }
99 }
100 if end {
101 break;
102 }
103 }
104 }
105 fn dijkstra(&mut self, s: usize, t: usize) -> bool {
106 self.dist.clear();
107 self.dist.resize(self.graph.vertices_size(), i64::MAX);
108 self.dist[s] = 0;
109 self.queue.clear();
110 self.queue.push((Reverse(0), s));
111 while let Some((Reverse(d), u)) = self.queue.pop() {
112 if d >= self.dist[t] {
113 break;
114 }
115 if self.dist[u] < d {
116 continue;
117 }
118 for a in self.graph.neighbors(u) {
119 if self.capacities[a.label] == 0 {
120 continue;
121 }
122 let ncost = (d.saturating_add(self.costs[a.label]))
123 .saturating_add(self.potential[u].saturating_sub(self.potential[a.to]));
124 if self.dist[a.to] > ncost {
125 debug_assert!(ncost >= d);
126 self.dist[a.to] = ncost;
127 self.prev_vertex[a.to] = u;
128 self.prev_edge[a.to] = a.label;
129 self.queue.push((Reverse(ncost), a.to));
130 }
131 }
132 }
133 self.dist[t] != i64::MAX
134 }
135 pub fn minimum_cost_flow_limited(&mut self, s: usize, t: usize, limit: u64) -> (u64, i64) {
137 let mut flow = 0;
138 let mut cost = 0;
139 if self.has_negedge {
140 self.bellman_ford(s);
141 }
142 while flow < limit && self.dijkstra(s, t) {
143 let shortest = self.dist[t];
144 for (p, d) in self.potential.iter_mut().zip(self.dist.iter()) {
145 *p = p.saturating_add((*d).min(shortest));
146 }
147 let mut f = limit - flow;
148 let mut v = t;
149 while v != s {
150 f = f.min(self.capacities[self.prev_edge[v]]);
151 v = self.prev_vertex[v];
152 }
153 flow += f;
154 cost += f as i64 * (self.potential[t] - self.potential[s]);
155 let mut v = t;
156 while v != s {
157 self.capacities[self.prev_edge[v]] -= f;
158 self.capacities[self.prev_edge[v] ^ 1] += f;
159 v = self.prev_vertex[v];
160 }
161 }
162 (flow, cost)
163 }
164 pub fn minimum_cost_flow(&mut self, s: usize, t: usize) -> (u64, i64) {
166 self.minimum_cost_flow_limited(s, t, u64::MAX)
167 }
168 pub fn get_flow(&self, eid: usize) -> u64 {
169 self.capacities[eid * 2 + 1]
170 }
171}