MaxFlowGoldbergTarjan.h

Go to the documentation of this file.

 
#pragma once
 
#include <ogdf/graphalg/MaxFlowModule.h>
 
//#define OGDF_GT_USE_GAP_RELABEL_HEURISTIC
#define OGDF_GT_USE_MAX_ACTIVE_LABEL
#ifdef OGDF_GT_USE_GAP_RELABEL_HEURISTIC
//#define OGDF_GT_GRH_STEPS 1   // gap relabel frequency: call gapRelabel() after OGDF_GT_GRH_STEPS relabel() operations (1 == off)
#endif
#define OGDF_GT_USE_PUSH_RELABEL_SECOND_STAGE
 
// world666 is much better without OGDF_GT_USE_PUSH_RELABEL_SECOND_STAGE
 
namespace ogdf {
 
 
template<typename TCap>
class MaxFlowGoldbergTarjan : public MaxFlowModule<TCap> {
    NodeArray<int> m_label;
    NodeArray<TCap> m_ex; // ex_f(v) values will be saved here to save runtime
#ifdef OGDF_GT_USE_MAX_ACTIVE_LABEL
    NodeArray<ListIterator<node>> m_activeLabelListPosition; // holds the iterator of every active node in the corresp. list of m_labeList
    Array<List<node>> m_activeLabelList; // array indexed by label, contains list of active nodes with that label
    int m_maxLabel = 0; // the maximum label among all active nodes
#endif
#ifdef OGDF_GT_USE_GAP_RELABEL_HEURISTIC
    NodeArray<ListIterator<node>> m_labelListPosition; // holds the iterator of every node in the corresp. list of m_labeList
    Array<List<node>> m_labelList; // array indexed by label, contains list of nodes with that label
#endif
 
    mutable List<node> m_cutNodes;
    mutable List<edge> m_cutEdges;
 
    inline TCap getCap(const edge e) const {
        return e->target() == *this->m_s ? 0 : (*this->m_cap)[e];
    }
 
    inline bool isResidualEdge(const adjEntry adj) const {
        const edge e = adj->theEdge();
        if (adj->theNode() == e->source()) {
            return this->m_et->less((*this->m_flow)[e], getCap(e));
        }
        return this->m_et->greater((*this->m_flow)[e], (TCap)0);
    }
 
    inline bool isAdmissible(const adjEntry adj) const {
        OGDF_ASSERT(adj);
        return isResidualEdge(adj) && m_label[adj->theNode()] == m_label[adj->twinNode()] + 1;
    }
 
    inline bool isActive(const node v) const {
        OGDF_ASSERT((v != *this->m_s && v != *this->m_t)
                || (m_label[*this->m_s] == this->m_G->numberOfNodes() && m_label[*this->m_t] == 0));
        return this->m_et->greater(m_ex[v], (TCap)0) && this->m_G->numberOfNodes() > m_label[v]
                && m_label[v] > 0;
    }
 
#ifdef OGDF_GT_USE_MAX_ACTIVE_LABEL
    inline void setActive(const node v) {
        const int label = m_label[v];
        OGDF_ASSERT(0 < label);
        OGDF_ASSERT(label < this->m_G->numberOfNodes());
        OGDF_ASSERT(!m_activeLabelListPosition[v].valid());
        m_activeLabelListPosition[v] = m_activeLabelList[label].pushBack(v);
        if (label > m_maxLabel) {
            m_maxLabel = label;
        }
    }
 
    inline void findNewMaxLabel() {
        while (m_maxLabel > 0 && m_activeLabelList[m_maxLabel].empty()) {
            --m_maxLabel;
        }
    }
 
    inline void setInactive(const node v) {
        OGDF_ASSERT(m_activeLabelListPosition[v].valid());
        m_activeLabelList[m_label[v]].del(m_activeLabelListPosition[v]);
        m_activeLabelListPosition[v] = nullptr;
        findNewMaxLabel();
    }
#endif
 
    // sets label of v, maintaining m_labelList (moves node v to the correct list in the array)
    inline void setLabel(const node v, int label) {
#ifdef OGDF_GT_USE_GAP_RELABEL_HEURISTIC
        if (m_labelListPosition[v].valid()) {
            m_labelList[m_label[v]].del(
                    m_labelListPosition[v]); // delete node from old list using noted position
        }
        m_labelListPosition[v] =
                m_labelList[label].pushBack(v); // push node to new list and update iterator
#endif
#ifdef OGDF_GT_USE_MAX_ACTIVE_LABEL
        if (m_activeLabelListPosition[v].valid()) {
            OGDF_ASSERT(0 < m_label[v]);
            OGDF_ASSERT(m_label[v] < this->m_G->numberOfNodes());
            setInactive(v);
        }
        m_label[v] = label; // update label
        if (v != *this->m_s && v != *this->m_t && isActive(v)) {
            setActive(v);
        }
#else
        m_label[v] = label; // update label
#endif
    }
 
#ifdef OGDF_GT_USE_GAP_RELABEL_HEURISTIC
    void gapRelabel() {
#   ifdef OGDF_GT_USE_MAX_ACTIVE_LABEL
        // XXX: this is a test but it seems to work and it seems to be fast!
        const int n = m_maxLabel + 1;
#   else
        const int n = this->m_G->numberOfNodes();
#   endif
        for (int i = 1; i < n - 1; ++i) {
            if (m_labelList[i].empty()) {
                for (int j = i + 1; j < n; ++j) {
                    while (!m_labelList[j].empty()) {
                        setLabel(m_labelList[j].front(), this->m_G->numberOfNodes());
                    }
                }
                break;
            }
        }
    }
#endif
 
    void push(const adjEntry adj) {
        const edge e = adj->theEdge();
        const node v = adj->theNode();
        if (v == e->source()) {
            const TCap value = min(m_ex[v], getCap(e) - (*this->m_flow)[e]);
            OGDF_ASSERT(this->m_et->geq(value, (TCap)0));
            (*this->m_flow)[e] += value;
            m_ex[v] -= value;
            m_ex[adj->twinNode()] += value;
        } else {
            const TCap value = min(m_ex[v], (*this->m_flow)[adj]);
            OGDF_ASSERT(this->m_et->geq(value, (TCap)0));
            (*this->m_flow)[adj] -= value;
            m_ex[v] -= value;
            m_ex[adj->twinNode()] += value;
        }
    }
 
    void globalRelabel() {
        // breadth-first search to relabel nodes with their respective distance to the sink in the residual graph
        const int n = this->m_G->numberOfNodes();
        NodeArray<int> dist(*this->m_G, n); // distance array
        List<node> queue; // reachable, not processed nodes
        dist[*this->m_t] = 0;
        queue.pushBack(*this->m_t);
        while (!queue.empty()) { // is there a node to check?
            node w = queue.popFrontRet();
            for (adjEntry adj : w->adjEntries) {
                node x = adj->twinNode();
                if (isResidualEdge(adj->twin()) && dist[x] == n) { // not already seen
                    dist[x] = dist[w] + 1; // set distance of node to sink
                    queue.pushBack(x);
                }
            }
        }
        // set distance of unreachable nodes to "number of nodes" thus making them inactive
        for (node w : this->m_G->nodes) {
            setLabel(w, dist[w]);
        }
    }
 
    void relabel(const node v) {
        int minLabel = this->m_G->numberOfNodes() - 1;
        for (adjEntry adj : v->adjEntries) {
            if (isResidualEdge(adj)) {
                const int label = m_label[adj->twinNode()];
                if (label < minLabel) {
                    minLabel = label;
                }
            }
        }
        if (minLabel + 1 != m_label[v]) { // == can happen after global relabel
            setLabel(v, minLabel + 1);
        }
    }
 
    void relabelStage2(const node v) {
        int minLabel = this->m_G->numberOfNodes() - 1;
        for (adjEntry adj : v->adjEntries) {
            if (isResidualEdge(adj)) {
                const int label = m_label[adj->twinNode()];
                if (label < minLabel) {
                    minLabel = label;
                }
            }
        }
        OGDF_ASSERT(minLabel + 1 != m_label[v]);
        m_label[v] = minLabel + 1;
    }
 
public:
    // first stage: push excess towards sink
    TCap computeValue(const EdgeArray<TCap>& cap, const node& s, const node& t) {
        // TODO: init this stuff in the module?
        this->m_s = &s;
        this->m_t = &t;
        this->m_cap = &cap;
        this->m_flow->init(*this->m_G, (TCap)0);
        OGDF_ASSERT(this->isFeasibleInstance());
 
        m_label.init(*this->m_G);
        m_ex.init(*this->m_G, 0);
#ifdef OGDF_GT_USE_MAX_ACTIVE_LABEL
        m_activeLabelListPosition.init(*this->m_G, nullptr);
        m_activeLabelList.init(1, this->m_G->numberOfNodes() - 1);
        m_maxLabel = 0;
#endif
#ifdef OGDF_GT_USE_GAP_RELABEL_HEURISTIC
        m_labelListPosition.init(*this->m_G, nullptr);
        m_labelList.init(this->m_G->numberOfNodes() + 1);
#endif
        m_cutNodes.clear();
 
        // initialize residual graph for first preflow
        for (edge e : this->m_G->edges) {
            if (e->source() == *this->m_s && e->target() != *this->m_s) { // ignore loops
                (*this->m_flow)[e] = getCap(e);
                m_ex[e->target()] += getCap(e); // "+" needed for the case of multigraphs
            }
        }
 
        if (*this->m_t == *this->m_s) {
            return (TCap)0;
        }
 
        NodeArray<adjEntry> curr(*this->m_G);
        for (node v = this->m_G->firstNode(); v; v = v->succ()) {
            curr[v] = v->firstAdj();
        }
 
        globalRelabel(); // initialize distance labels
 
        int relCount = 0; // counts the relabel operations for the global relabeling heuristic
#ifdef OGDF_GT_USE_MAX_ACTIVE_LABEL
        while (m_maxLabel != 0) {
            OGDF_ASSERT(!m_activeLabelList[m_maxLabel].empty());
            const node v = m_activeLabelList[m_maxLabel].front();
            OGDF_ASSERT(m_maxLabel == m_label[v]);
            OGDF_ASSERT(m_activeLabelListPosition[v] == m_activeLabelList[m_maxLabel].begin());
#else
        List<node> active;
        for (adjEntry adj : (*this->m_s)->adjEntries) {
            node w = adj->theEdge()->target();
            if (w != *this->m_s) {
                active.pushBack(w);
            }
        }
        while (!active.empty()) {
            const node v = active.front();
#endif
            adjEntry& adj = curr[v];
            if (v == *this->m_s || v == *this->m_t || !isActive(v)) {
                // source, sink or not active: remove activity status
#ifdef OGDF_GT_USE_MAX_ACTIVE_LABEL
                setInactive(v);
#else
                active.popFront();
#endif
            } else {
                while (this->m_et->greater(m_ex[v], (TCap)0)) {
                    if (isAdmissible(adj)) {
                        // push and adjacent node becomes active
#ifdef OGDF_GT_USE_MAX_ACTIVE_LABEL
                        const node w = adj->twinNode();
                        if (w != *this->m_s && w != *this->m_t && !isActive(w)) {
                            // w will become active after push
                            setActive(w);
                        }
                        push(adj);
                        if (v != *this->m_s && !isActive(v)) {
                            setInactive(v);
                        }
#else
                        push(adj);
                        active.pushBack(adj->twinNode());
#endif
                    } else {
                        if (adj != v->lastAdj()) {
                            adj = adj->succ();
                        } else { // end of adjacency list
                            adj = v->firstAdj();
                            relabel(v);
                            ++relCount;
#ifdef OGDF_GT_USE_GAP_RELABEL_HEURISTIC
                            // only gapRelabel if we do not do a globalRelabel directly afterwards
                            if (relCount != this->m_G->numberOfNodes()
#   if (OGDF_GT_GRH_STEPS > 1)
                                    && relCount % OGDF_GT_GRH_STEPS
                                            == 0 // obey frequency of gap relabel heuristic
#   endif
                            ) {
                                gapRelabel();
                            }
#endif
                            break;
                        }
                    }
                }
                if (relCount == this->m_G->numberOfNodes()) {
                    relCount = 0;
                    globalRelabel();
                }
            }
        }
 
        TCap result = 0;
        for (adjEntry adj : (*this->m_t)->adjEntries) {
            edge e = adj->theEdge();
            if (e->target() == *this->m_t) {
                result += (*this->m_flow)[e];
            } else {
                result -= (*this->m_flow)[e];
            }
        }
        return result;
    }
 
    // second stage: push excess that has not reached the sink back towards source
    void computeFlowAfterValue() {
        List<node> active;
#ifdef OGDF_GT_USE_PUSH_RELABEL_SECOND_STAGE
        NodeArray<adjEntry> curr(*this->m_G);
        for (node v = this->m_G->firstNode(); v; v = v->succ()) {
            curr[v] = v->firstAdj();
            m_label[v] = 1;
            if (this->m_et->greater(m_ex[v], (TCap)0) && v != *this->m_s && v != *this->m_t) {
                active.pushBack(v);
            }
        }
        if (active.empty()) {
            return;
        }
 
        m_label[*this->m_s] = 0;
        while (!active.empty()) {
            node v = active.front();
            if (v == *this->m_s || v == *this->m_t || !isActive(v)) {
                active.popFront();
            } else {
                adjEntry& adj = curr[v];
                if (isAdmissible(adj)) {
                    push(adj);
                    active.pushBack(adj->twinNode());
                } else {
                    if (adj == v->lastAdj()) {
                        // no admissible outgoing edge found -> relabel node!
                        relabelStage2(v);
                        adj = v->firstAdj();
#   if 0
                        // node is still active but move it to the end of the queue
                        // (don't know if this is really necessary)
                        active.popFront();
                        active.pushBack(v);
#   endif
                    } else {
                        adj = adj->succ();
                    }
                }
            }
        }
#else // USE_PUSH_RELABEL_SECOND_STAGE
        m_ex[*this->m_s] = m_ex[*this->m_t] = 0;
        for (node v = this->m_G->firstNode(); v; v = v->succ()) {
            if (this->m_et->greater(m_ex[v], (TCap)0)) {
                active.pushBack(v);
            }
        }
        while (!active.empty()) {
            const node v = active.popFrontRet();
            if (this->m_et->greater(m_ex[v], (TCap)0) && v != *this->m_s && v != *this->m_t) {
                for (adjEntry adj = v->firstAdj(); adj; adj = adj->succ()) {
                    const edge e = adj->theEdge();
                    const node u = e->source();
                    if (u != v) { // e is incoming edge
                        if (this->m_et->greater(m_ex[v], (TCap)0) && isResidualEdge(adj)) {
                            push(adj);
                            if (u != *this->m_s) {
                                active.pushFront(u);
                            }
                        }
                    }
                }
            }
        }
#endif // USE_PUSH_RELABEL_SECOND_STAGE
    }
 
    using MaxFlowModule<TCap>::useEpsilonTest;
    using MaxFlowModule<TCap>::init;
    using MaxFlowModule<TCap>::computeFlow;
    using MaxFlowModule<TCap>::computeFlowAfterValue;
    using MaxFlowModule<TCap>::MaxFlowModule;
};
 
}