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

#include "minFlow.h"

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

        // Create separate flow graph for use in first phase
        // Just like fg1, but with extra sink-to-source edge
        Mflograph fg2(fg1.n(),fg1.m()+1, fg1.src(), fg1.snk());
        fg2.copyFrom(fg1);
        int tcap = 0;
        for (edge e = fg2.first(); e != 0; e = fg2.next(e))
                tcap += fg2.cap(fg2.tail(e),e);
        edge snkSrcEdge = fg2.join(fg2.snk(),fg2.src());
        if (snkSrcEdge == 0)
                fatal("minFlow: internal error, can't create sink/source edge");
        fg2.setCapacity(snkSrcEdge,tcap); fg2.setMinFlo(snkSrcEdge,0);
        fg = &fg2;

        // 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);
        }
        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; }
                add2cycle(e);
        }
        floVal = fg2.f(fg2.snk(),snkSrcEdge);

        // Transfer flow from fg2 to fg1
        fg2.remove(snkSrcEdge);
        fg1.copyFrom(fg2);
        fg = &fg1;

        // 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);
                        }
                }
        }
        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 = u; edge px = pEdge[x];
        while (x != v) {
                vertex y = fg->mate(x,px);
                f = min(f,fg->res(y,px));
                x = y; px = pEdge[x];
        }
        // add flow to saturate cycle
        fg->addFlow(u,e,f);
        x = u; px = pEdge[x];
        while (x != v) {
                vertex y = fg->mate(x,px);
                fg->addFlow(y,px,f);
                x = y; px = pEdge[x];
        }
        return f;
}