rs-graph 0.20.1

// Copyright (c) 2015-2021 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::traits::*;
//! use rs_graph::maxflow::pushrelabel;
//! use rs_graph::{EdgeVec, Net};
//! use rs_graph::string::{Data, from_ascii};
//!
//! let Data { graph: g, weights: upper, nodes } = from_ascii::<Net>(r"
//!      a---2-->b
//!     @|\      ^\
//!    / | \     | 4
//!   5  |  \    |  \
//!  /   |   |   |   @
//! s    1   1   2    t
//!  \   |   |   |   @
//!   5  |    \  |  /
//!    \ |     \ | 5
//!     @v      @|/
//!      c---2-->d
//!     ").unwrap();
//!
//! let s = nodes[&'s'];
//! let t = nodes[&'t'];
//! let v1 = nodes[&'a'];
//! let v2 = nodes[&'b'];
//! let v3 = nodes[&'c'];
//! let v4 = nodes[&'d'];
//!
//! let (value, flow, mut mincut) = pushrelabel(&g, s, t, |e| upper[e.index()]);
//!
//! assert_eq!(value, 5);
//! assert!(g.edges().all(|e| flow[e] >= 0 && flow[e] <= upper[e.index()]));
//! assert!(g.nodes().filter(|&u| u != s && u != t).all(|u| {
//!     g.outedges(u).map(|(e,_)| flow[e]).sum::<usize>() ==
//!     g.inedges(u).map(|(e,_)| flow[e]).sum::<usize>()
//! }));
//!
//! mincut.sort_by_key(|u| u.index());
//! assert_eq!(mincut, vec![v1, s, v3]);
//! ```
//!
//! ```
//! use rs_graph::traits::*;
//! use rs_graph::maxflow::pushrelabel;
//! use rs_graph::{EdgeVec, Net};
//! use rs_graph::string::{Data, from_ascii};
//!
//! let Data { graph: g, weights: upper, nodes } = from_ascii::<Net>(r"
//!                ---8-->a---10---
//!               /       |        \
//!              /        1  --3--  |
//!             /         | /     \ |
//!            /          v@       \v
//!      ---->b-----9---->c----8--->d----
//!     /      \         @         @^    \
//!   18        ---6--  /         / |     33
//!   /               \/         /  |      \
//!  /                /\    -----   |       @
//! s           --5--- |   /        |        t
//!  \         /       |  /         |       @
//!   27      |  ----2-|--         /       /
//!    \      | /      |  /----8---       7
//!     \     |/       @  |              /
//!      ---->e----9-->f------6---->g----
//!            \          |        @
//!             \         |       /
//!              --5----->h---4---
//!     ").unwrap();
//!
//! let s = nodes[&'s'];
//! let t = nodes[&'t'];
//!
//! assert_eq!(g.num_edges(), 18);
//!
//! let (value, flow, mut mincut) = pushrelabel(&g, s, t, |e| upper[e.index()]);
//! assert_eq!(value, 29);
//!
//! let mincutval = g
//!     .edges()
//!     .filter(|&e| mincut.contains(&g.src(e)) && !mincut.contains(&g.snk(e)))
//!     .map(|e| upper[e.index()])
//!     .sum::<usize>();
//! assert_eq!(value, mincutval);
//!
//! mincut.sort_by_key(|u| u.index());
//! assert_eq!(mincut, "bcsef".chars().map(|v| nodes[&v]).collect::<Vec<_>>());
//!
//! assert!(g.edges().all(|e| flow[e] >= 0 && flow[e] <= upper[e.index()]));
//! assert!(g.nodes().filter(|&u| u != s && u != t).all(|u| {
//!     g.outedges(u).map(|(e,_)| flow[e]).sum::<usize>() ==
//!     g.inedges(u).map(|(e,_)| flow[e]).sum::<usize>()
//! }));
//! ```

use crate::maxflow::MaxFlow;

use crate::adapters::Network;
use crate::traits::{Directed, GraphIterator, GraphSize, GraphType, IndexDigraph, IndexGraph};
use crate::vec::{EdgeVec, NodeVec};

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

use crate::num::traits::NumAssign;

/// The push-relabel algorithm.
///
/// This struct contains all algorithmic working data.
pub struct PushRelabel<'a, G, F>
where
    G: 'a + IndexDigraph<'a>,
{
    /// The underlying graph (extended to a network) the flow problem is solved
    /// on.
    g: Network<'a, G>,

    /// Data associated with each node.
    nodes: NodeVec<'a, &'a G, NodeInfo<'a, G, F>>,
    /// Current flow on each edge.
    flow: EdgeVec<'a, Network<'a, G>, F>,
    /// The buckets containing the nodes of a specific height.
    buckets: Vec<Bucket<G::Node>>,
    /// The queue of nodes for a BFS.
    queue: VecDeque<G::Node>,
    /// The largest height of an active node.
    largest_act: usize,
    /// The flow value.
    value: F,
    /// The number of relabel operations performed during the algorithm.
    pub cnt_relabel: usize,
    /// Whether to use the global relabelling heuristic.
    pub use_global_relabelling: bool,
}

/// Data associated with a node.
struct NodeInfo<'a, G, Flow>
where
    G: IndexDigraph<'a>,
{
    /// The current height of the node.
    height: usize,
    /// The excess of flow of the node.
    excess: Flow,
    /// The next active node in the linked list for the current height.
    next_act: Option<G::Node>,
    /// The current edge to be considered in [`discharge`].
    current_edge: Option<(
        <Network<'a, G> as GraphType<'a>>::Edge,
        <Network<'a, G> as Directed<'a>>::OutIt,
    )>,
}

/// A bucket containing nodes of some height.
#[derive(Clone)]
struct Bucket<Node> {
    /// The first active node of this height.
    ///
    /// The active nodes are kept in a singly linked list, this is the first
    /// node in that list.
    first_act: Option<Node>,
    /// Number of inactive nodes of this height.
    num_inact: usize,
}

impl<Node> Bucket<Node> {
    /// Returns `true` if the bucket is empty.
    ///
    /// This means that there are neither active nor inactive nodes of the
    /// buckets height.
    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: IndexDigraph<'a>,
    F: NumAssign + Ord + Copy,
{
    type Graph = G;
    type Flow = F;

    /// Return a new push-relabel algorithm data structure for the digraph `g`.
    fn new(g: &'a G) -> Self {
        let n = g.num_nodes();
        PushRelabel {
            g: Network::new(g),
            nodes: NodeVec::new_with(g, |_| NodeInfo {
                height: 0,
                excess: F::zero(),
                next_act: None,
                current_edge: None,
            }),
            flow: EdgeVec::new(Network::new(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,
        }
    }

    /// Return a reference to the underlying graph.
    fn as_graph(&self) -> &'a Self::Graph {
        self.g.as_graph()
    }

    /// Return the flow value.
    ///
    /// The function returns 0 if the flow has not been computed, yet.
    fn value(&self) -> F {
        self.value
    }

    /// Return the flow value over some edge.
    ///
    /// The function returns 0 if the flow has not been computed, yet.
    fn flow(&self, e: G::Edge) -> F {
        self.flow[self.g.from_edge(e)]
    }

    /// Run the push-relabel algorithm from some source to some sink node.
    fn solve<Us>(&mut self, src: G::Node, snk: G::Node, upper: Us)
    where
        Us: Fn(G::Edge) -> Self::Flow,
    {
        assert_ne!(
            self.g.node_id(src),
            self.g.node_id(snk),
            "Source and sink node must not be equal"
        );

        let n = self.g.num_nodes();
        self.cnt_relabel = 0;
        self.init_preflow(src, upper);

        self.update_heights(src, snk, false);

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

        while self.largest_act != 0 && self.largest_act != n {
            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);
            if self.cnt_relabel >= lvl_relabel && self.largest_act < n {
                self.update_heights(src, snk, false);
                lvl_relabel += n;
            }

            if self.largest_act == 0 {
                // End of phase I. Compute exact labels one last time, this time
                // from the source.
                self.update_heights(src, snk, true);
            }
        }

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

    /// Return the minimal cut.
    fn mincut(&self) -> Vec<G::Node> {
        // The minimal cut consists of the nodes with height >= n
        let n = self.g.num_nodes();
        self.g.nodes().filter(|u| self.nodes[*u].height >= n).collect()
    }
}

impl<'a, G, F> PushRelabel<'a, G, F>
where
    G: IndexDigraph<'a>,
    F: NumAssign + Ord + Copy,
{
    /// Initialize preflow algorithm.
    ///
    /// All edges leaving the source node are saturated, the source's
    /// height is set to `n`, all other heights are set to `0`.
    fn init_preflow<Us>(&mut self, src: G::Node, upper: Us)
    where
        Us: Fn(G::Edge) -> F,
    {
        for u in self.g.nodes() {
            self.nodes[u] = NodeInfo {
                height: 0,
                excess: F::zero(),
                next_act: None,
                current_edge: None,
            };
        }
        for e in self.g.edges() {
            if e.is_forward() {
                self.flow[e] = F::zero();
            } else {
                self.flow[e] = upper(e.edge());
            }
        }

        // send maximal flow out of source
        for (e, v) in self.g.outedges(src) {
            let f = e.reverse();
            let ub = self.flow[f];
            self.flow[e] = ub;
            self.flow[f] = F::zero();
            self.nodes[v].excess += ub;
        }

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

    /// Compute exact labels.
    ///
    /// This function does a bfs from the sink (phase I) or source (phase II) to
    /// compute exact labels.
    fn update_heights(&mut self, src: G::Node, snk: G::Node, from_src: bool) {
        let n = self.g.num_nodes();

        // Put all nodes to height 2*n
        for u in self.g.nodes() {
            self.nodes[u].height = n + n;
            // we need to reset the iterators for correctness
            self.nodes[u].current_edge = None;
        }

        // Except for the source and and the sink
        self.nodes[snk].height = 0;
        self.nodes[src].height = n;

        for b in if from_src {
            &mut self.buckets[n + 1..]
        } else {
            &mut self.buckets[..n]
        } {
            b.first_act = None;
            b.num_inact = 0;
        }

        // source and sink count as inactive, just in case
        self.buckets[0].num_inact = 1;
        self.buckets[n].num_inact = 1;

        // find correct labels by BFS from sink
        self.queue.clear();
        self.queue.push_back(snk);

        self.largest_act = 0;
        loop {
            while let Some(v) = self.queue.pop_front() {
                let h = self.nodes[v].height + 1;
                for (e, u) in self.g.outedges(v) {
                    let udata = &mut self.nodes[u];
                    if udata.height > h && self.flow[e] > F::zero() {
                        udata.height = h;
                        self.queue.push_back(u);
                        // add node to appropriate bucket
                        if udata.excess > F::zero() {
                            // u is active
                            udata.next_act = self.buckets[h].first_act;
                            self.buckets[h].first_act = Some(u);
                            self.largest_act = h;
                            debug_assert!(from_src || self.largest_act < n);
                        } else {
                            // u is inactive
                            self.buckets[h].num_inact += 1;
                        }
                    }
                }
            }
            // possibly start again from the source to find the perfect
            // labels for phase II
            if self.largest_act >= n || !from_src {
                break;
            }
            self.largest_act = n;
            self.queue.push_back(src);
        }
    }

    /// Discharges node `u`.
    ///
    /// This function does a sequence of push and relabel operations for an
    /// active node `u` until its excess reaches 0.
    ///
    /// In phase I, the function may also with `u` having nonzero excess if the
    /// height of u is at least `n`. In this case `u` gets disconnected from the
    /// sink and will not be considered again until phase II.
    #[allow(clippy::many_single_char_names)]
    fn discharge(&mut self, u: G::Node) {
        let n = self.g.num_nodes();

        loop {
            // The minimal height of the reachable neighbors.
            // If `u` gets relabelled, it gets height `n_neighbor + 1`.
            let mut h_neighbor = 2 * n;
            let h_u = self.nodes[u].height;

            // Start at current edge (if any) or restart at beginning.
            let (first_edge, mut rest) = if let Some((e, r)) = self.nodes[u].current_edge.take() {
                // just return the current edge
                ((e, self.g.snk(e)), r)
            } else {
                // start at the beginning
                let mut it = self.g.out_iter(u);
                if let Some(e) = it.next(&self.g) {
                    (e, it)
                } else {
                    // there is not beginning, just stop
                    return;
                }
            };

            let mut cur = first_edge;

            loop {
                let (e, v) = cur;
                // skip non-admissible edges
                let f = e.reverse();
                if !self.flow[f].is_zero() {
                    if h_u == self.nodes[v].height + 1 {
                        // Push along edge e
                        let df = min(self.nodes[u].excess, self.flow[f]);

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

                        if self.nodes[v].excess.is_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);
                            debug_assert!(self.buckets[h].num_inact > 0, "height:{} n:{}", h, n);
                            self.buckets[h].num_inact -= 1;
                        }

                        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
                            self.buckets[h_u].num_inact += 1;
                            self.largest_act = (0..h_u + 1)
                                .rev()
                                .find(|&h| self.buckets[h].first_act.is_some())
                                .unwrap_or(0);
                            // save current edge
                            self.nodes[u].current_edge = Some((cur.0, rest));
                            return;
                        }
                    } else {
                        h_neighbor = h_neighbor.min(self.nodes[v].height);
                    }
                }

                // edge is full, go to next
                if let Some(e) = rest.next(&self.g) {
                    cur = e
                } else {
                    break;
                }
            }

            // Update new height for neighbors before `first_edge`.
            let end = first_edge.0;
            for (e, v) in self.g.outedges(u) {
                if e == end {
                    break;
                }
                if !self.flow[e.reverse()].is_zero() {
                    h_neighbor = h_neighbor.min(self.nodes[v].height);
                }
            }

            // we ran out of admissible edges but node still has positive excess, relabel node
            if !self.relabel(u, h_neighbor + 1) {
                // node has not got active again (i.e. its label is too high now), stop
                break;
            }
        }
    }

    /// The relabel operation.
    ///
    /// Relabel `u` to height `h_new`.
    ///
    /// The function returns `true` iff the node remains active (i.e. if it
    /// could be discharged again).
    ///
    /// In phase I the function returns `false` if the height of `u` gets at
    /// least `n`. In this case `u` is disconnected from the sink and will not
    /// be discharged again until phase II.
    fn relabel(&mut self, u: G::Node, h_new: usize) -> bool {
        debug_assert!(self.nodes[u].excess > F::zero());

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

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

        debug_assert!(h_new > h_old);

        // *** The GAP heuristic ***

        // Check if we are still in phase I.
        if h_old < n {
            let mut h_new = h_new;

            // Test whether the node's old bucket is empty now.
            // This would create a gap.
            if self.buckets[h_old].is_empty() {
                // 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..n {
                    if self.buckets[h].is_empty() {
                        // all higher buckets are empty, too
                        break;
                    }
                    self.buckets[h].first_act = None;
                    self.buckets[h].num_inact = 0;
                }

                // Relabel all nodes to n+1
                //
                // Actually, this is not very fast. We could keep a list of *all*
                // nodes of each height, then this would be possible in O(1)
                // amortized time (each node is relabelled to n+1 at most once).
                //
                // However, because each node can be the cause of the disconnection
                // at most once (when it is relabelled and creates a gap), the
                // running time is at most O(n^2) anyway.
                for u in self.g.nodes() {
                    if h_old < self.nodes[u].height {
                        self.nodes[u].height = n + 1;
                    }
                }

                // The node has now label n + 1.
                h_new = n + 1;
            }

            // This node has now a too large label for phase I.
            // There is no need to discharge the node again. Instead we look for
            // a new active node.
            if h_new >= n {
                debug_assert_eq!(self.largest_act, h_old);

                self.nodes[u].height = h_new;

                // Find largest remaining bucket with an active node.
                //
                // Phase I ends once this reaches 0.
                self.largest_act = (0..h_old + 1)
                    .rev()
                    .find(|&h| self.buckets[h].first_act.is_some())
                    .unwrap_or(0);

                // The current node has been relabelled to n+1, so no further discharge is necessary.
                return false;
            }
        }

        // The node remains active, is must have the largest height now.
        self.nodes[u].height = h_new;
        self.largest_act = h_new;
        true
    }
}

#[cfg(test)]
mod tests {
    use crate::maxflow::pushrelabel;
    use crate::traits::*;
    use crate::{Buildable, Builder, Net};

    #[test]
    fn test_pushrelabel() {
        let mut g = Net::new_builder();
        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 g = g.into_graph();
        let (value, flow, _) = pushrelabel(&g, s, t, |e| upper[e.index()]);

        assert_eq!(value, 11);
        assert!(g.edges().all(|e| flow[e] >= 0 && flow[e] <= upper[e.index()]));
        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>()));
    }
}