drwnMaxFlow.h

/******************************************************************************
** DARWIN: A FRAMEWORK FOR MACHINE LEARNING RESEARCH AND DEVELOPMENT
** Distributed under the terms of the BSD license (see the LICENSE file)
** Copyright (c) 2007-2015, Stephen Gould
** All rights reserved.
**
******************************************************************************
** FILENAME:    drwnMaxFlow.h
** AUTHOR(S):   Stephen Gould <stephen.gould@anu.edu.au>
**
*****************************************************************************/

#pragma once

#include <cassert>
#include <vector>
#include <map>
#include <deque>

using namespace std;

class drwnMaxFlow;

// drwnMaxFlow --------------------------------------------------------------

class drwnMaxFlow {
 protected:
    static const unsigned char FREE = 0x00;
    static const unsigned char SOURCE = 0x01;
    static const unsigned char TARGET = 0x02;

    typedef map<int, pair<unsigned, unsigned> > _drwnCapacitatedEdge;

    vector<double> _sourceEdges; 
    vector<double> _targetEdges; 
    vector<double> _edgeWeights; 
    vector<_drwnCapacitatedEdge> _nodes; 
    double _flowValue;           

    vector<unsigned char> _cut;  

 public:
    drwnMaxFlow(unsigned maxNodes = 0) : _flowValue(0.0) {
        _sourceEdges.reserve(maxNodes);
        _targetEdges.reserve(maxNodes);
        _edgeWeights.reserve(2 * maxNodes);
        _nodes.reserve(maxNodes);
    }

    virtual ~drwnMaxFlow() {
        // do nothing
    }

    size_t numNodes() const { return _nodes.size(); }

    virtual void reset();
    virtual void clear();

    inline int addNodes(unsigned n = 1) {
        int nodeId = (int)_nodes.size();
        _nodes.resize(_nodes.size() + n);
        _sourceEdges.resize(_nodes.size(), 0.0);
        _targetEdges.resize(_nodes.size(), 0.0);
        return nodeId;
    }

    void addConstant(double c) {
        _flowValue += c;
    }

    inline void addSourceEdge(int u, double cap) {
        DRWN_ASSERT((u >= 0) && (u < (int)_nodes.size()));
        if (cap < 0.0) { _flowValue += cap; _targetEdges[u] -= cap; }
        else _sourceEdges[u] += cap;
    }

    inline void addTargetEdge(int u, double cap) {
        DRWN_ASSERT((u >= 0) && (u < (int)_nodes.size()));
        if (cap < 0.0) { _flowValue += cap; _sourceEdges[u] -= cap; }
        else _targetEdges[u] += cap;
    }

    inline void addEdge(int u, int v, double cap_uv, double cap_vu = 0.0) {
        DRWN_ASSERT((u >= 0) && (u < (int)_nodes.size()));
        DRWN_ASSERT((v >= 0) && (v < (int)_nodes.size()));
        DRWN_ASSERT(u != v);

        _drwnCapacitatedEdge::const_iterator it = _nodes[u].find(v);
        if (it == _nodes[u].end()) {
            DRWN_ASSERT(cap_uv + cap_vu >= 0.0);
            if (cap_uv < 0.0) {
                _nodes[u].insert(make_pair(v, make_pair(_edgeWeights.size(), _edgeWeights.size() + 1)));
                _nodes[v].insert(make_pair(u, make_pair(_edgeWeights.size() + 1, _edgeWeights.size())));
                _edgeWeights.push_back(0.0);
                _edgeWeights.push_back(cap_uv + cap_vu);
                _sourceEdges[u] -= cap_uv;
                _targetEdges[v] -= cap_uv;
                _flowValue += cap_uv;
            } else if (cap_vu < 0.0) {
                _nodes[u].insert(make_pair(v, make_pair(_edgeWeights.size(), _edgeWeights.size() + 1)));
                _nodes[v].insert(make_pair(u, make_pair(_edgeWeights.size() + 1, _edgeWeights.size())));
                _edgeWeights.push_back(cap_uv + cap_vu);
                _edgeWeights.push_back(0.0);
                _sourceEdges[v] -= cap_vu;
                _targetEdges[u] -= cap_vu;
                _flowValue += cap_vu;
            } else {
                _nodes[u].insert(make_pair(v, make_pair(_edgeWeights.size(), _edgeWeights.size() + 1)));
                _nodes[v].insert(make_pair(u, make_pair(_edgeWeights.size() + 1, _edgeWeights.size())));
                _edgeWeights.push_back(cap_uv);
                _edgeWeights.push_back(cap_vu);
            }
        } else {
            const double w_u = _edgeWeights[it->second.first] += cap_uv;
            const double w_v = _edgeWeights[it->second.second] += cap_vu;
            DRWN_ASSERT(w_u + w_v >= 0.0);
            if (w_u < 0.0) {
                _edgeWeights[it->second.first] = 0.0;
                _edgeWeights[it->second.second] += w_u;
                _sourceEdges[u] -= w_u;
                _targetEdges[v] -= w_u;
                _flowValue += w_u;
            } else if (w_v < 0.0) {
                _edgeWeights[it->second.first] += w_v;
                _edgeWeights[it->second.second] = 0.0;
                _sourceEdges[v] -= w_v;
                _targetEdges[u] -= w_v;
                _flowValue += w_v;
            }
        }
    }

    virtual double solve() = 0;

    bool inSetS(int u) const { return (_cut[u] == SOURCE); }
    bool inSetT(int u) const { return (_cut[u] == TARGET); }

    double operator()(int u, int v) const;

 protected:
    void preAugmentPaths();
};

// drwnEdmondsKarpMaxFlow ---------------------------------------------------

class drwnEdmondsKarpMaxFlow : public drwnMaxFlow {
 public:
    drwnEdmondsKarpMaxFlow(unsigned maxNodes = 0) : drwnMaxFlow(maxNodes) {
        // do nothing
    }
    virtual ~drwnEdmondsKarpMaxFlow() {
        // do nothing
    }

    virtual double solve();

 protected:
    void findCuts();
};

// drwnBKMaxFlow ------------------------------------------------------------

class drwnBKMaxFlow : public drwnMaxFlow {
 protected:
    static const int TERMINAL = -1;  

    vector<int> _parents;
    vector<pair<double *, double *> > _weightptrs;
    drwnIndexQueue _activeList;

 public:
    drwnBKMaxFlow(unsigned maxNodes = 0) : drwnMaxFlow(maxNodes) {
        // do nothing
    }
    virtual ~drwnBKMaxFlow() {
        // do nothing
    }

    virtual void reset();
    virtual void clear();
    virtual double solve();

 protected:
    void initializeTrees();
    pair<int, int> expandTrees();
    void augmentBKPath(const pair<int, int>& path, deque<int>& orphans);
    void adoptOrphans(deque<int>& orphans);
    bool isAncestor(int u, int v) const;
};