rs-graph 0.14.1

// Copyright (c) 2015, 2016, 2017, 2018 Frank Fischer <frank-fischer@shadow-soft.de>
//
// This program is free software: you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation, either version 3 of the
// License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see  <http://www.gnu.org/licenses/>
//

//! This module implements a push relabel algorithm for solving max
//! flow problems.
//!
//! This implementation uses the gap heuristic and the global
//! relabelling heuristic.
//!
//! # Example
//!
//! ```
//! use rs_graph::{Builder, Graph, Digraph, Net, EdgeVec, NumberedNode};
//! use rs_graph::maxflow::pushrelabel;
//!
//! let mut g = Net::new();
//! let mut upper = vec![];
//! let s = g.add_node();
//! let t = g.add_node();
//! let v1 = g.add_node();
//! let v2 = g.add_node();
//! let v3 = g.add_node();
//! let v4 = g.add_node();
//! g.add_edge(s, v1); upper.push(5);
//! g.add_edge(s, v3); upper.push(5);
//! g.add_edge(v1, v2); upper.push(2);
//! g.add_edge(v1, v3); upper.push(1);
//! g.add_edge(v1, v4); upper.push(1);
//! g.add_edge(v2, t); upper.push(4);
//! g.add_edge(v3, v4); upper.push(2);
//! g.add_edge(v4, v2); upper.push(2);
//! g.add_edge(v4, t); upper.push(5);
//! let upper: EdgeVec<_> = upper.into();
//!
//! let (value, flow, mut mincut) = pushrelabel(&g, s, t, &upper);
//!
//! assert_eq!(value, 5);
//! assert!(g.edges().all(|e| flow[e] >= 0 && flow[e] <= upper[e]));
//! assert!(g.nodes().filter(|&u| u != s && u != t).all(|u| {
//!     g.outedges(u).map(|(e,_)| flow[e])
//!                  .sum::<isize>() ==
//!     g.inedges(u).map(|(e,_)| flow[e]).sum::<isize>()
//! }));
//!
//! mincut.sort_by_key(|u| u.index());
//! assert_eq!(mincut, vec![s, v1, v3]);
//! ```

use maxflow::MaxFlow;

use graph::{IndexBiEdgeVec, IndexNetwork, IndexNodeVec};
use EdgeMap;

use std::cmp::min;
use std::collections::VecDeque;

use num::traits::NumAssign;

pub struct PushRelabel<'a, G, F>
where
    G: 'a + IndexNetwork<'a>,
{
    g: &'a G,

    nodes: IndexNodeVec<'a, G, NodeInfo<G::Node, F>>,
    flow: IndexBiEdgeVec<'a, G, F>,
    buckets: Vec<Bucket<G::Node>>,
    queue: VecDeque<G::Node>,
    largest_act: usize,
    value: F,
    pub cnt_relabel: usize,
    pub use_global_relabelling: bool,
}

#[derive(Clone)]
struct NodeInfo<Node, Flow> {
    height: usize,
    excess: Flow,
    next_act: Option<Node>,
}

#[derive(Clone)]
struct Bucket<Node> {
    // first active node
    first_act: Option<Node>,
    // number of inactive nodes
    num_inact: usize,
}

impl<Node> Bucket<Node> {
    fn is_empty(&self) -> bool {
        self.first_act.is_none() && self.num_inact == 0
    }
}

impl<'a, G, F> MaxFlow<'a> for PushRelabel<'a, G, F>
where
    G: IndexNetwork<'a>,
    F: NumAssign + Ord + Copy,
{
    type Graph = G;
    type Flow = F;

    fn new(g: &'a G) -> Self {
        let n = g.num_nodes();
        PushRelabel {
            g: g,

            nodes: idxnodevec![g; NodeInfo { height: 0, excess: F::zero(), next_act: None}],
            flow: idxbiedgevec![g; F::zero()],
            buckets: vec![
                Bucket {
                    first_act: None,
                    num_inact: 0,
                };
                n * 2
            ],
            queue: VecDeque::with_capacity(n),
            largest_act: 0,
            value: F::zero(),

            cnt_relabel: 0,
            use_global_relabelling: true,
        }
    }

    fn as_graph(&self) -> &'a Self::Graph {
        self.g
    }

    fn value(&self) -> F {
        self.value
    }

    fn flow(&self, e: G::Edge) -> F {
        self.flow[self.g.forward(e)]
    }

    fn solve<Us>(&mut self, src: G::Node, snk: G::Node, upper: Us)
    where
        Us: EdgeMap<'a, G, F>,
    {
        assert_ne!(
            self.g.node_id(src),
            self.g.node_id(snk),
            "Source and sink node must not be equal"
        );

        self.cnt_relabel = 0;
        self.init_preflow(src, upper);

        self.update_heights(snk);

        let mut lvl_relabel = if self.use_global_relabelling {
            self.g.num_nodes()
        } else {
            usize::max_value()
        };

        let mut next_edges = idxnodevec![self.g, self.g.nodes().map(|u| self.g.neighs(u)).collect()];

        while self.largest_act != 0 && self.largest_act != self.g.num_nodes() {
            let l = self.largest_act;
            let u = self.buckets[l].first_act.unwrap();
            self.buckets[l].first_act = self.nodes[u].next_act;
            self.discharge(u, &mut next_edges);
            if self.cnt_relabel >= lvl_relabel && self.largest_act < self.g.num_nodes() {
                self.update_heights(snk);
                lvl_relabel += self.g.num_nodes();
            }

            if self.largest_act == 0 {
                // end of phase I
                // we put now all nodes in bucket n+1
                let mut first_act = None;
                let mut num_inact = 0;
                for u in self.g.nodes() {
                    if u != src && u != snk {
                        self.nodes[u].height = self.g.num_nodes() + 1;
                        if self.nodes[u].excess > F::zero() {
                            self.nodes[u].next_act = first_act;
                            first_act = Some(u);
                        } else {
                            num_inact += 1;
                        }
                    }
                }
                self.buckets[self.g.num_nodes()] = Bucket {
                    first_act: None,
                    num_inact: 1,
                };
                self.buckets[self.g.num_nodes() + 1] = Bucket { first_act, num_inact };
                if first_act.is_some() {
                    self.largest_act = self.g.num_nodes() + 1;
                }
            }
        }

        self.value = self.nodes[snk].excess;
    }

    fn mincut(&self) -> Vec<G::Node> {
        // find min-cut with bfs from source
        let n = self.g.num_nodes();
        if let Some(src) = self.g.nodes().find(|&u| self.nodes[u].height == n) {
            let mut queue = VecDeque::with_capacity(n);
            let mut seen = idxnodevec![self.g; false];
            seen[src] = true;
            queue.push_back(src);
            while let Some(u) = queue.pop_front() {
                for (e, v) in self.g.outedges(u) {
                    if !seen[v] && !self.flow[self.g.reverse(e)].is_zero() {
                        seen[v] = true;
                        queue.push_back(v);
                    }
                }
            }
            self.g.nodes().filter(|&u| seen[u]).collect()
        } else {
            vec![]
        }
    }
}

impl<'a, G, F> PushRelabel<'a, G, F>
where
    G: IndexNetwork<'a>,
    F: NumAssign + Ord + Copy,
{
    /// Initialize preflow algorithm.
    ///
    /// All edges leaving the source node are saturated, the source's
    /// label is set to `n`, all other labels are set to `0`.
    ///
    /// The function returns the initialized data structures.
    fn init_preflow<Us>(&mut self, src: G::Node, upper: Us)
    where
        Us: EdgeMap<'a, G, F>,
    {
        for u in self.g.nodes() {
            self.nodes[u] = NodeInfo {
                height: 0,
                excess: F::zero(),
                next_act: None,
            };
        }
        for e in self.g.edges() {
            self.flow[self.g.forward(e)] = F::zero();
            self.flow[self.g.backward(e)] = upper[e];
        }

        for (e, v) in self.g.outedges(src) {
            self.flow[e] = upper[e];
            self.flow[self.g.reverse(e)] = F::zero();
            self.nodes[v].excess += upper[e];
        }

        self.nodes[src].height = self.g.num_nodes();
    }

    /// Compute exact labels.
    ///
    /// This function does a bfs from the sink to compute exact labels.
    fn update_heights(&mut self, snk: G::Node) {
        for u in self.g.nodes() {
            self.nodes[u].height = self.g.num_nodes();
        }
        for b in &mut self.buckets {
            b.first_act = None;
            b.num_inact = 0;
        }

        self.queue.clear();
        self.queue.push_back(snk);
        self.nodes[snk].height = 0;
        self.largest_act = 0;
        while let Some(v) = self.queue.pop_front() {
            let h = self.nodes[v].height;
            for (f, u) in self.g.neighs(v) {
                if self.flow[f] > F::zero() && self.nodes[u].height > h + 1 {
                    self.nodes[u].height = h + 1;
                    self.queue.push_back(u);
                }
            }
            // add node to appropriate bucket
            if self.nodes[v].excess > F::zero() {
                self.nodes[v].next_act = self.buckets[h].first_act;
                self.buckets[h].first_act = Some(v);
                self.largest_act = h;
            } else {
                self.buckets[h].num_inact += 1;
            }
        }
    }

    #[cfg_attr(feature = "cargo-clippy", allow(many_single_char_names))]
    fn discharge(&mut self, u: G::Node, next_edges: &mut IndexNodeVec<'a, G, G::NeighIter>) {
        while self.nodes[u].height < self.g.num_nodes()  // node is active in phase I
            || self.largest_act > self.g.num_nodes()
        // node is active in phase II
        {
            while let Some((e, v)) = next_edges[u].next() {
                // skip non-admissible edges
                let f = self.g.reverse(e);
                if self.flow[f].is_zero() || self.nodes[u].height != self.nodes[v].height + 1 {
                    continue;
                }

                let df = min(self.nodes[u].excess, self.flow[f]);

                debug_assert_eq!(self.nodes[u].height, self.nodes[v].height + 1);
                debug_assert!(df > F::zero());

                if self.nodes[v].excess == F::zero() {
                    // sink node becomes active
                    let h = self.nodes[v].height;
                    self.nodes[v].next_act = self.buckets[h].first_act;
                    self.buckets[h].first_act = Some(v);
                    self.buckets[h].num_inact -= 1;
                    next_edges[v] = self.g.neighs(v);
                }

                self.flow[e] += df;
                self.flow[f] -= df;
                self.nodes[u].excess -= df;
                self.nodes[v].excess += df;

                // check if node is fully discharged
                if self.nodes[u].excess.is_zero() {
                    // node is inactive now
                    let mut h = self.nodes[u].height;
                    self.buckets[h].num_inact += 1;
                    while h > 0 && self.buckets[h].first_act.is_none() {
                        h -= 1;
                    }
                    self.largest_act = h;
                    return;
                }
            }

            // we ran out of admissible edges but node remains active, relabel node
            self.relabel(u);
            // reset iterator
            next_edges[u] = self.g.neighs(u);
        }
    }

    fn relabel(&mut self, u: G::Node) {
        debug_assert!(self.nodes[u].excess > F::zero());

        let h_old = self.nodes[u].height;
        let h_new = self
            .g
            .neighs(u)
            .filter_map(|(e, v)| {
                if self.flow[self.g.reverse(e)] > F::zero() {
                    Some(self.nodes[v].height)
                } else {
                    None
                }
            })
            .min()
            .unwrap() + 1;
        debug_assert!(h_new > h_old);

        self.nodes[u].height = h_new;

        // Test whether the node's old bucket is empty now.
        if self.buckets[h_old].is_empty() && h_old < self.g.num_nodes() {
            // It is empty now, so u is disconnected from the sink.
            // We remove all nodes with higher label from all buckets.
            for h in h_old + 1..self.g.num_nodes() {
                if self.buckets[h].is_empty() {
                    // all higher buckets are empty, too
                    break;
                }
                self.buckets[h] = Bucket {
                    first_act: None,
                    num_inact: 0,
                };
            }

            // find largest remaining bucket with an active node
            for h in (0..h_old).rev() {
                if self.buckets[h].first_act.is_some() {
                    self.largest_act = h;
                    return;
                }
            }

            // no active node found, so, end of phase I
            self.largest_act = 0;
            return;
        }

        // this node has now a too large label for phase I -> stop
        if self.largest_act < self.g.num_nodes() && h_new >= self.g.num_nodes() {
            debug_assert_eq!(self.largest_act, h_old);
            let mut h = self.largest_act;
            while h > 0 && self.buckets[h].first_act.is_none() {
                h -= 1;
            }
            self.largest_act = h;
            return;
        }

        // node remains active, is must have the largest height now
        self.largest_act = h_new;

        // count number of relabellings
        self.cnt_relabel += 1;
    }
}

#[cfg(test)]
mod tests {
    use maxflow::pushrelabel;
    use {Builder, Digraph, EdgeVec, Graph, Net};

    #[test]
    fn test_pushrelabel() {
        let mut g = Net::new();
        let mut upper = vec![];
        let s = g.add_node();
        let t = g.add_node();
        let v1 = g.add_node();
        let v2 = g.add_node();
        let v3 = g.add_node();
        let v4 = g.add_node();
        g.add_edge(s, v1);
        upper.push(15);
        g.add_edge(s, v3);
        upper.push(10);
        g.add_edge(v1, v2);
        upper.push(6);
        g.add_edge(v1, v3);
        upper.push(7);
        g.add_edge(v2, t);
        upper.push(5);
        g.add_edge(v2, v4);
        upper.push(2);
        g.add_edge(v3, v2);
        upper.push(11);
        g.add_edge(v3, v4);
        upper.push(4);
        g.add_edge(v4, v2);
        upper.push(4);
        g.add_edge(v4, t);
        upper.push(20);
        let upper: EdgeVec<_> = upper.into();
        let (value, flow, _) = pushrelabel(&g, s, t, &upper);

        assert_eq!(value, 11);
        assert!(g.edges().all(|e| flow[e] >= 0 && flow[e] <= upper[e]));
        assert!(
            g.nodes()
                .filter(|&u| u != s && u != t)
                .all(|u| g.outedges(u).map(|(e, _)| flow[e]).sum::<isize>()
                    == g.inedges(u).map(|(e, _)| flow[e]).sum::<isize>())
        );
    }
}