/Users/jst/src/grafalgo/cpp/graphAlgorithms/maxFlo/minFlowx.cpp

#include "minFlow.h"

minFlow::minFlow(Mflograph& fg1, int& floVal) : fg(&fg1) {
        pEdge = new edge[fg->n()+1];

        // attempt to find a flow that satisfies all min capacities
        // first, add edges with unsatisfied min capacities to a todo list
        UiList todo(fg->m());
        for (edge e = fg->first(); e != 0; e = fg->next(e)) {
                vertex u = fg->tail(e);
                if (fg->f(u,e) < fg->minFlo(e))
                        todo.addLast(e);
        }
        floVal = 0;
        while (!todo.empty()) {
                edge e = todo.first();
                vertex u = fg->tail(e); vertex v = fg->head(e);
                if (fg->f(u,e) >= fg->minFlo(e)) {
                        todo.removeFirst(); continue;
                }
                if (!findCycle(e)) { floVal = -1; return; }
                floVal += add2cycle(e);
        }

        // now, push additional flow from source to sink
        while(findPath()) {
                floVal += augment(); 
        }
}

minFlow::~minFlow() { delete [] pEdge; }

bool minFlow::findPath() {
        vertex u,v; edge e;
        UiList queue(fg->n());

        for (u = 1; u <= fg->n(); u++) pEdge[u] = 0;
        queue.addLast(fg->src());
        while (!queue.empty()) {
                u = queue.first(); queue.removeFirst();
                for (e = fg->firstAt(u); e != 0; e = fg->nextAt(u,e)) {
                        v = fg->mate(u,e);
                        if (fg->res(u,e) > 0 && pEdge[v] == 0 && 
                            v != fg->src()) {
                                pEdge[v] = e;
                                if (v == fg->snk()) return true;
                                queue.addLast(v);
                        }
                }
        }
        return false;
}

int minFlow::augment() {
        vertex u, v; edge e; flow f;

        // determine residual capacity of path
        f = BIGINT;
        u = fg->snk(); e = pEdge[u];
        while (u != fg->src()) {
                v = fg->mate(u,e);
                f = min(f,fg->res(v,e));
                u = v; e = pEdge[u];
        }
        // add flow to saturate path
        u = fg->snk(); e = pEdge[u];
        while (u != fg->src()) {
                v = fg->mate(u,e);
                fg->addFlow(v,e,f);
                u = v; e = pEdge[u];
        }
        return f;
}

bool minFlow::findCycle(edge e) {
        vertex u = fg->tail(e); vertex v = fg->head(e);
        for (vertex x = 1; x <= fg->n(); x++) pEdge[x] = 0;

        UiList queue(fg->n()); queue.addLast(v);
        while (!queue.empty()) {
                vertex x = queue.first(); queue.removeFirst();
                for (edge ex = fg->firstAt(x); ex != 0; ex = fg->nextAt(x,ex)) {
                        vertex y = fg->mate(x,ex);
                        if (fg->res(x,ex) > 0 && pEdge[y] == 0 && y != v) {
                                pEdge[y] = ex;
                                if (y == u) return true;
                                queue.addLast(y);
                        }
                }
                // if search reaches sink, continue searching from source
                if (x != fg->snk() || pEdge[fg->src()] != 0) continue;
                pEdge[fg->src()] = -1; // special value to indicate src/snk link
                if (u == fg->src()) return true;
                x = fg->src();
                for (edge ex = fg->firstAt(x); ex != 0; ex = fg->nextAt(x,ex)) {
                        vertex y = fg->mate(x,ex);
                        if (fg->res(x,ex) > 0 && pEdge[y] == 0 && y != v) {
                                pEdge[y] = ex;
                                if (y == u) return true;
                                queue.addLast(y);
                        }
                }
        }
        return false;
}

int minFlow::add2cycle(edge e) {
        vertex u = fg->tail(e); vertex v = fg->head(e);

        // determine residual capacity of cycle
        flow f = fg->res(u,e);
        vertex x = (pEdge[u] != -1 ? u : fg->snk());
        edge px = pEdge[x];
        while (x != v) {
                vertex y = fg->mate(x,px);
                f = min(f,fg->res(y,px));
                x = (pEdge[y] != -1 ? y : fg->snk());
                px = pEdge[x];
        }
        // add flow to saturate cycle
        fg->addFlow(u,e,f);
        x = (pEdge[u] != -1 ? u : fg->snk());
        px = pEdge[x];
        while (x != v) {
                vertex y = fg->mate(x,px);
                fg->addFlow(y,px,f);
                x = (pEdge[y] != -1 ? y : fg->snk());
                px = pEdge[x];
        }
        return f;
}