logicaffeine_proof/
min_cost_flow.rs

1//! Minimum-cost flow, exact and certified: the dual of every scheduling and retiming problem whose constraints are
2//! differences of two variables, which is how a register cut of least area at a target clock is found.
3//!
4//! [`min_cost_flow`] runs the primal network simplex on a strongly feasible spanning tree (Cunningham, *A network
5//! simplex method*, 1976): an artificial root joined to every node, block-search pricing, and the tree kept as a
6//! thread list so a pivot touches only the subtree it moves (the layout of LEMON's `NetworkSimplex`, Kovács,
7//! *Minimum-cost flow algorithms: an experimental evaluation*, 2015). Integer data give an integer optimum.
8//!
9//! Every answer carries what proves it, and [`check`] replays it trusting nothing the search did:
10//! - an optimum carries its flow and a potential `π` under which every arc below capacity has reduced cost
11//!   `cost + π(from) − π(to) ≥ 0` and every arc with flow has it `≤ 0`, which is optimality by complementary
12//!   slackness;
13//! - an infeasible problem carries a node set whose supply exceeds the capacity leaving it;
14//! - one without an optimum carries a cycle of uncapacitated arcs of negative cost.
15
16/// The bound on an arc's cost.
17pub const COST_LIMIT: i64 = 1 << 32;
18
19/// The bound on an arc's capacity and on the supplies' total size.
20pub const FLOW_LIMIT: i64 = 1 << 62;
21
22/// An arc: `capacity` units at most (`None` for no bound), at `cost` a unit.
23#[derive(Clone, Copy, Debug, PartialEq, Eq)]
24pub struct Arc {
25    pub from: usize,
26    pub to: usize,
27    pub capacity: Option<i64>,
28    pub cost: i64,
29}
30
31/// A problem's answer, with its certificate.
32#[derive(Clone, Debug, PartialEq, Eq)]
33pub enum Solution {
34    Optimal { flow: Vec<i64>, potential: Vec<i128>, cost: i128 },
35    /// Nodes whose supply is more than the arcs leaving them can carry.
36    Infeasible { cut: Vec<usize> },
37    /// No optimum: uncapacitated arcs, in order, around a cycle of negative cost, around which any flow meeting the
38    /// supplies gets cheaper without end.
39    Unbounded { cycle: Vec<usize> },
40}
41
42/// Why a problem was not solved.
43#[derive(Clone, Copy, Debug, PartialEq, Eq)]
44pub enum FlowError {
45    /// An arc names a node past the count.
46    OutOfRange { arc: usize },
47    /// An arc's capacity is negative or past [`FLOW_LIMIT`], or its cost past [`COST_LIMIT`].
48    OutOfBounds { arc: usize },
49    /// The supplies do not sum to zero, or their sizes sum past [`FLOW_LIMIT`].
50    Unbalanced,
51    /// More than 2³¹ nodes or arcs.
52    TooLarge,
53}
54
55fn validate(nodes: usize, supply: &[i64], arcs: &[Arc]) -> Result<(), FlowError> {
56    if nodes >= 1 << 31 || arcs.len() >= 1 << 31 || supply.len() != nodes {
57        return Err(if supply.len() != nodes { FlowError::Unbalanced } else { FlowError::TooLarge });
58    }
59    for (k, a) in arcs.iter().enumerate() {
60        if a.from >= nodes || a.to >= nodes {
61            return Err(FlowError::OutOfRange { arc: k });
62        }
63        if a.capacity.is_some_and(|c| !(0..FLOW_LIMIT).contains(&c)) || a.cost.abs() >= COST_LIMIT {
64            return Err(FlowError::OutOfBounds { arc: k });
65        }
66    }
67    let size: i128 = supply.iter().map(|&s| (s as i128).abs()).sum();
68    if supply.iter().map(|&s| s as i128).sum::<i128>() != 0 || size >= FLOW_LIMIT as i128 {
69        return Err(FlowError::Unbalanced);
70    }
71    Ok(())
72}
73
74const INF: i64 = i64::MAX;
75const UP: i8 = 1;
76const DOWN: i8 = -1;
77const TREE: i8 = 0;
78const LOWER: i8 = 1;
79const UPPER: i8 = -1;
80
81/// The simplex's state: arcs `0..m` the problem's, `m..m + n` the artificial ones, node `n` the root.
82struct Simplex {
83    m: usize,
84    source: Vec<usize>,
85    target: Vec<usize>,
86    cap: Vec<i64>,
87    cost: Vec<i128>,
88    flow: Vec<i64>,
89    state: Vec<i8>,
90    pi: Vec<i128>,
91    parent: Vec<usize>,
92    pred: Vec<usize>,
93    pred_dir: Vec<i8>,
94    thread: Vec<usize>,
95    rev_thread: Vec<usize>,
96    succ_num: Vec<usize>,
97    last_succ: Vec<usize>,
98    dirty: Vec<usize>,
99    next_arc: usize,
100    block: usize,
101    in_arc: usize,
102    join: usize,
103    u_in: usize,
104    v_in: usize,
105    u_out: usize,
106    delta: i64,
107}
108
109const NONE: usize = usize::MAX;
110
111impl Simplex {
112    fn new(nodes: usize, supply: &[i64], arcs: &[Arc]) -> Simplex {
113        let (n, m) = (nodes, arcs.len());
114        let all = m + n;
115        let max_cost = arcs.iter().map(|a| (a.cost as i128).abs()).max().unwrap_or(0);
116        let art = (max_cost + 1) * (n as i128 + 1);
117        let mut s = Simplex {
118            m,
119            source: Vec::with_capacity(all),
120            target: Vec::with_capacity(all),
121            cap: Vec::with_capacity(all),
122            cost: Vec::with_capacity(all),
123            flow: vec![0; all],
124            state: vec![LOWER; all],
125            pi: vec![0; n + 1],
126            parent: vec![NONE; n + 1],
127            pred: vec![NONE; n + 1],
128            pred_dir: vec![0; n + 1],
129            thread: vec![0; n + 1],
130            rev_thread: vec![0; n + 1],
131            succ_num: vec![1; n + 1],
132            last_succ: vec![0; n + 1],
133            dirty: Vec::new(),
134            next_arc: 0,
135            block: ((m as f64).sqrt() as usize).max(10),
136            in_arc: 0,
137            join: 0,
138            u_in: 0,
139            v_in: 0,
140            u_out: 0,
141            delta: 0,
142        };
143        for a in arcs {
144            s.source.push(a.from);
145            s.target.push(a.to);
146            s.cap.push(a.capacity.unwrap_or(INF));
147            s.cost.push(a.cost as i128);
148        }
149        let root = n;
150        s.thread[root] = 0;
151        s.rev_thread[0] = root;
152        s.succ_num[root] = n + 1;
153        s.last_succ[root] = if n == 0 { root } else { root - 1 };
154        if n == 0 {
155            s.thread[root] = root;
156            s.rev_thread[root] = root;
157        }
158        for u in 0..n {
159            let e = m + u;
160            s.parent[u] = root;
161            s.pred[u] = e;
162            s.thread[u] = u + 1;
163            s.rev_thread[u + 1] = u;
164            s.last_succ[u] = u;
165            s.cap.push(INF);
166            s.state[e] = TREE;
167            if supply[u] >= 0 {
168                s.pred_dir[u] = UP;
169                s.source.push(u);
170                s.target.push(root);
171                s.flow[e] = supply[u];
172                s.cost.push(0);
173                s.pi[u] = 0;
174            } else {
175                s.pred_dir[u] = DOWN;
176                s.source.push(root);
177                s.target.push(u);
178                s.flow[e] = -supply[u];
179                s.cost.push(art);
180                s.pi[u] = art;
181            }
182        }
183        s
184    }
185
186    fn reduced(&self, e: usize) -> i128 {
187        self.cost[e] + self.pi[self.source[e]] - self.pi[self.target[e]]
188    }
189
190    /// The arc of the most negative violation in the next block that has one.
191    fn entering(&mut self) -> bool {
192        let m = self.m;
193        if m == 0 {
194            return false;
195        }
196        let mut best: i128 = 0;
197        let mut count = self.block;
198        let mut e = self.next_arc;
199        for _ in 0..m {
200            let c = self.state[e] as i128 * self.reduced(e);
201            if c < best {
202                best = c;
203                self.in_arc = e;
204            }
205            e = if e + 1 == m { 0 } else { e + 1 };
206            count -= 1;
207            if count == 0 {
208                if best < 0 {
209                    break;
210                }
211                count = self.block;
212            }
213        }
214        self.next_arc = e;
215        best < 0
216    }
217
218    fn join_node(&mut self) {
219        let (mut u, mut v) = (self.source[self.in_arc], self.target[self.in_arc]);
220        while u != v {
221            if self.succ_num[u] < self.succ_num[v] {
222                u = self.parent[u];
223            } else {
224                v = self.parent[v];
225            }
226        }
227        self.join = u;
228    }
229
230    /// The arc leaving the tree: the first bottleneck on the cycle's way back to the join from the side the
231    /// entering arc starts, else the last on the other side, so the tree stays strongly feasible. Whether one does.
232    fn leaving(&mut self) -> bool {
233        let e = self.in_arc;
234        let (first, second) = if self.state[e] == LOWER { (self.source[e], self.target[e]) } else { (self.target[e], self.source[e]) };
235        self.delta = self.cap[e];
236        let mut result = 0;
237        let mut u = first;
238        while u != self.join {
239            let a = self.pred[u];
240            let d = if self.pred_dir[u] == DOWN { if self.cap[a] == INF { INF } else { self.cap[a] - self.flow[a] } } else { self.flow[a] };
241            if d < self.delta {
242                self.delta = d;
243                self.u_out = u;
244                result = 1;
245            }
246            u = self.parent[u];
247        }
248        let mut u = second;
249        while u != self.join {
250            let a = self.pred[u];
251            let d = if self.pred_dir[u] == UP { if self.cap[a] == INF { INF } else { self.cap[a] - self.flow[a] } } else { self.flow[a] };
252            if d <= self.delta {
253                self.delta = d;
254                self.u_out = u;
255                result = 2;
256            }
257            u = self.parent[u];
258        }
259        if result == 1 {
260            self.u_in = first;
261            self.v_in = second;
262        } else {
263            self.u_in = second;
264            self.v_in = first;
265        }
266        result != 0
267    }
268
269    fn change_flow(&mut self, change: bool) {
270        let e = self.in_arc;
271        if self.delta > 0 {
272            let val = self.state[e] as i64 * self.delta;
273            self.flow[e] += val;
274            let mut u = self.source[e];
275            while u != self.join {
276                let a = self.pred[u];
277                self.flow[a] -= self.pred_dir[u] as i64 * val;
278                u = self.parent[u];
279            }
280            let mut u = self.target[e];
281            while u != self.join {
282                let a = self.pred[u];
283                self.flow[a] += self.pred_dir[u] as i64 * val;
284                u = self.parent[u];
285            }
286        }
287        if change {
288            self.state[e] = TREE;
289            let out = self.pred[self.u_out];
290            self.state[out] = if self.flow[out] == 0 { LOWER } else { UPPER };
291        } else {
292            self.state[e] = -self.state[e];
293        }
294    }
295
296    /// Hang the subtree under `u_out` from `v_in` by the entering arc, re-rooted at `u_in`.
297    fn update_tree(&mut self) {
298        let (u_in, v_in, u_out, join, in_arc) = (self.u_in, self.v_in, self.u_out, self.join, self.in_arc);
299        let old_rev_thread = self.rev_thread[u_out];
300        let old_succ_num = self.succ_num[u_out];
301        let old_last_succ = self.last_succ[u_out];
302        let v_out = self.parent[u_out];
303        if u_in == u_out {
304            self.parent[u_in] = v_in;
305            self.pred[u_in] = in_arc;
306            self.pred_dir[u_in] = if u_in == self.source[in_arc] { UP } else { DOWN };
307            if self.thread[v_in] != u_out {
308                let mut after = self.thread[old_last_succ];
309                self.thread[old_rev_thread] = after;
310                self.rev_thread[after] = old_rev_thread;
311                after = self.thread[v_in];
312                self.thread[v_in] = u_out;
313                self.rev_thread[u_out] = v_in;
314                self.thread[old_last_succ] = after;
315                self.rev_thread[after] = old_last_succ;
316            }
317        } else {
318            let thread_continue = if old_rev_thread == v_in { self.thread[old_last_succ] } else { self.thread[v_in] };
319            let mut stem = u_in;
320            let mut par_stem = v_in;
321            let mut last = self.last_succ[u_in];
322            let mut after = self.thread[last];
323            self.thread[v_in] = u_in;
324            self.dirty.clear();
325            self.dirty.push(v_in);
326            while stem != u_out {
327                let next_stem = self.parent[stem];
328                self.thread[last] = next_stem;
329                self.dirty.push(last);
330                let before = self.rev_thread[stem];
331                self.thread[before] = after;
332                self.rev_thread[after] = before;
333                self.parent[stem] = par_stem;
334                par_stem = stem;
335                stem = next_stem;
336                last = if self.last_succ[stem] == self.last_succ[par_stem] { self.rev_thread[par_stem] } else { self.last_succ[stem] };
337                after = self.thread[last];
338            }
339            self.parent[u_out] = par_stem;
340            self.thread[last] = thread_continue;
341            self.rev_thread[thread_continue] = last;
342            self.last_succ[u_out] = last;
343            if old_rev_thread != v_in {
344                self.thread[old_rev_thread] = after;
345                self.rev_thread[after] = old_rev_thread;
346            }
347            for i in 0..self.dirty.len() {
348                let u = self.dirty[i];
349                let t = self.thread[u];
350                self.rev_thread[t] = u;
351            }
352            let mut tmp_sc = 0usize;
353            let tmp_ls = self.last_succ[u_out];
354            let mut u = u_out;
355            while u != u_in {
356                let p = self.parent[u];
357                self.pred[u] = self.pred[p];
358                self.pred_dir[u] = -self.pred_dir[p];
359                tmp_sc = tmp_sc + self.succ_num[u] - self.succ_num[p];
360                self.succ_num[u] = tmp_sc;
361                self.last_succ[p] = tmp_ls;
362                u = p;
363            }
364            self.pred[u_in] = in_arc;
365            self.pred_dir[u_in] = if u_in == self.source[in_arc] { UP } else { DOWN };
366            self.succ_num[u_in] = old_succ_num;
367        }
368        let up_limit_out = if self.last_succ[join] == v_in { join } else { NONE };
369        let last_succ_out = self.last_succ[u_out];
370        let mut u = v_in;
371        while u != NONE && self.last_succ[u] == v_in {
372            self.last_succ[u] = last_succ_out;
373            u = self.parent[u];
374        }
375        if join != old_rev_thread && v_in != old_rev_thread {
376            let mut u = v_out;
377            while u != up_limit_out && u != NONE && self.last_succ[u] == old_last_succ {
378                self.last_succ[u] = old_rev_thread;
379                u = self.parent[u];
380            }
381        } else if last_succ_out != old_last_succ {
382            let mut u = v_out;
383            while u != up_limit_out && u != NONE && self.last_succ[u] == old_last_succ {
384                self.last_succ[u] = last_succ_out;
385                u = self.parent[u];
386            }
387        }
388        let mut u = v_in;
389        while u != join {
390            self.succ_num[u] += old_succ_num;
391            u = self.parent[u];
392        }
393        let mut u = v_out;
394        while u != join {
395            self.succ_num[u] -= old_succ_num;
396            u = self.parent[u];
397        }
398    }
399
400    fn update_potential(&mut self) {
401        let u_in = self.u_in;
402        let sigma = self.pi[self.v_in] - self.pi[u_in] - self.pred_dir[u_in] as i128 * self.cost[self.in_arc];
403        let end = self.thread[self.last_succ[u_in]];
404        let mut u = u_in;
405        while u != end {
406            self.pi[u] += sigma;
407            u = self.thread[u];
408        }
409    }
410
411    /// The cycle the entering arc closes through the tree, its arcs in the direction flow would go around it.
412    fn cycle(&self) -> Vec<usize> {
413        let e = self.in_arc;
414        let (first, second) = if self.state[e] == LOWER { (self.source[e], self.target[e]) } else { (self.target[e], self.source[e]) };
415        let mut back = Vec::new();
416        let mut u = second;
417        while u != self.join {
418            back.push(self.pred[u]);
419            u = self.parent[u];
420        }
421        let mut down = Vec::new();
422        let mut u = first;
423        while u != self.join {
424            down.push(self.pred[u]);
425            u = self.parent[u];
426        }
427        down.reverse();
428        let mut cycle = vec![e];
429        cycle.extend(back);
430        cycle.extend(down);
431        cycle
432    }
433}
434
435/// The least-cost flow meeting `supply` (positive where flow starts, negative where it ends) over `arcs`.
436pub fn min_cost_flow(nodes: usize, supply: &[i64], arcs: &[Arc]) -> Result<Solution, FlowError> {
437    validate(nodes, supply, arcs)?;
438    let mut s = Simplex::new(nodes, supply, arcs);
439    while s.entering() {
440        s.join_node();
441        let change = s.leaving();
442        if s.delta == INF {
443            return Ok(Solution::Unbounded { cycle: s.cycle() });
444        }
445        s.change_flow(change);
446        if change {
447            s.update_tree();
448            s.update_potential();
449        }
450    }
451    let m = arcs.len();
452    if (m..m + nodes).any(|e| s.flow[e] != 0) {
453        return Ok(Solution::Infeasible { cut: infeasible_cut(nodes, arcs, &s) });
454    }
455    let flow: Vec<i64> = s.flow[..m].to_vec();
456    let cost = arcs.iter().zip(&flow).map(|(a, &f)| a.cost as i128 * f as i128).sum();
457    Ok(Solution::Optimal { flow, potential: s.pi[..nodes].to_vec(), cost })
458}
459
460/// The nodes the residual network reaches from those whose supply did not leave: no residual path reaches a node
461/// whose demand went unmet, or the optimum would have used it, so their supply exceeds what can leave them.
462fn infeasible_cut(nodes: usize, arcs: &[Arc], s: &Simplex) -> Vec<usize> {
463    let m = arcs.len();
464    let mut reached = vec![false; nodes];
465    let mut stack: Vec<usize> = (0..nodes).filter(|&u| s.source[m + u] == u && s.flow[m + u] > 0).collect();
466    for &u in &stack {
467        reached[u] = true;
468    }
469    let mut out: Vec<Vec<usize>> = vec![Vec::new(); nodes];
470    let mut inn: Vec<Vec<usize>> = vec![Vec::new(); nodes];
471    for (k, a) in arcs.iter().enumerate() {
472        out[a.from].push(k);
473        inn[a.to].push(k);
474    }
475    while let Some(u) = stack.pop() {
476        for &k in &out[u] {
477            let v = arcs[k].to;
478            if !reached[v] && arcs[k].capacity.map_or(true, |c| s.flow[k] < c) {
479                reached[v] = true;
480                stack.push(v);
481            }
482        }
483        for &k in &inn[u] {
484            let v = arcs[k].from;
485            if !reached[v] && s.flow[k] > 0 {
486                reached[v] = true;
487                stack.push(v);
488            }
489        }
490    }
491    (0..nodes).filter(|&u| reached[u]).collect()
492}
493
494/// Whether `solution` is the problem's answer, by its certificate alone.
495pub fn check(nodes: usize, supply: &[i64], arcs: &[Arc], solution: &Solution) -> bool {
496    if validate(nodes, supply, arcs).is_err() {
497        return false;
498    }
499    match solution {
500        Solution::Optimal { flow, potential, cost } => {
501            if flow.len() != arcs.len() || potential.len() != nodes {
502                return false;
503            }
504            let mut net: Vec<i128> = supply.iter().map(|&s| s as i128).collect();
505            for (a, &f) in arcs.iter().zip(flow) {
506                if f < 0 || a.capacity.is_some_and(|c| f > c) {
507                    return false;
508                }
509                net[a.from] -= f as i128;
510                net[a.to] += f as i128;
511            }
512            if net.iter().any(|&x| x != 0) {
513                return false;
514            }
515            let slack = arcs.iter().zip(flow).all(|(a, &f)| {
516                let Some(rc) = (a.cost as i128).checked_add(potential[a.from]).and_then(|x| x.checked_sub(potential[a.to])) else { return false };
517                (a.capacity.is_some_and(|c| f == c) || rc >= 0) && (f == 0 || rc <= 0)
518            });
519            slack && *cost == arcs.iter().zip(flow).map(|(a, &f)| a.cost as i128 * f as i128).sum::<i128>()
520        }
521        Solution::Infeasible { cut } => {
522            let mut inside = vec![false; nodes];
523            for &u in cut {
524                if u >= nodes || std::mem::replace(&mut inside[u], true) {
525                    return false;
526                }
527            }
528            let mut leaving: i128 = 0;
529            for a in arcs {
530                if inside[a.from] && !inside[a.to] {
531                    match a.capacity {
532                        Some(c) => leaving += c as i128,
533                        None => return false,
534                    }
535                }
536            }
537            cut.iter().map(|&u| supply[u] as i128).sum::<i128>() > leaving
538        }
539        Solution::Unbounded { cycle } => {
540            !cycle.is_empty()
541                && cycle.iter().all(|&k| k < arcs.len() && arcs[k].capacity.is_none())
542                && cycle.iter().enumerate().all(|(i, &k)| arcs[k].to == arcs[cycle[(i + 1) % cycle.len()]].from)
543                && cycle.iter().map(|&k| arcs[k].cost as i128).sum::<i128>() < 0
544        }
545    }
546}
547
548#[cfg(test)]
549#[path = "../../logicaffeine_tests/unit/logicaffeine_proof/min_cost_flow/tests.rs"]
550mod tests;
logicaffeine_proof/min_cost_flow.rs

logicaffeine_proof/
min_cost_flow.rs
