Flujo máximo - algoritmo MPM

El algoritmo MPM (Malhotra, Pramodh-Kumar y Maheshwari) resuelve el problema de flujo máximo en $O(V^3)$ . Este algoritmo es similar al algoritmo de Dinic.

Algoritmo

Como el algoritmo de Dinic, MPM corre en fases; durante cada fase encontramos el flujo bloqueante en la red por capas de la red residual de $G$ . La diferencia principal con Dinic es cómo encontramos el flujo bloqueante. Consideremos la red por capas $L$ . Para cada nodo definimos su potencial interno (inner potential) y su potencial externo (outer potential) como:

$\begin{aligned} p_{i n} (v) & = \sum_{(u, v) \in L} (c (u, v) - f (u, v)) \\ p_{o u t} (v) & = \sum_{(v, u) \in L} (c (v, u) - f (v, u)) \end{aligned}$

También fijamos $p_{in}(s) = p_{out}(t) = \infty$ . Dados $p_{in}$ y $p_{out}$ definimos el potencial como $p(v) = min(p_{in}(v), p_{out}(v))$ . Llamamos a un nodo $r$ un nodo de referencia (reference node) si $p(r) = min{p(v)}$ . Consideremos un nodo de referencia $r$ . Afirmamos que el flujo se puede aumentar en $p(r)$ de modo que $p(r)$ pase a ser $0$ . Es cierto porque $L$ es acíclica, así que podemos empujar el flujo hacia afuera de $r$ por las aristas salientes y llegará a $t$ porque cada nodo tiene suficiente potencial externo para empujar el flujo hacia afuera cuando lo alcanza. De forma similar, podemos tirar del flujo desde $s$ . La construcción del flujo bloqueante se basa en este hecho. En cada iteración encontramos un nodo de referencia y empujamos el flujo de $s$ a $t$ a través de $r$ . Este proceso se puede simular con BFS. Todos los arcos completamente saturados se pueden eliminar de $L$ , ya que de todos modos no se usarán más adelante en esta fase. Del mismo modo, se pueden eliminar todos los nodos distintos de $s$ y $t$ sin arcos salientes o entrantes.

Cada fase funciona en $O(V^2)$ porque hay a lo sumo $V$ iteraciones (porque al menos se elimina el nodo de referencia elegido), y en cada iteración eliminamos todas las aristas por las que pasamos excepto a lo sumo $V$ . Sumando, obtenemos $O(V^2 + E) = O(V^2)$ . Como hay menos de $V$ fases (véase la demostración aquí), MPM funciona en $O(V^3)$ en total.

Implementación


struct MPM{
    struct FlowEdge{
        int v, u;
        long long cap, flow;
        FlowEdge(){}
        FlowEdge(int _v, int _u, long long _cap, long long _flow)
            : v(_v), u(_u), cap(_cap), flow(_flow){}
        FlowEdge(int _v, int _u, long long _cap)
            : v(_v), u(_u), cap(_cap), flow(0ll){}
    };
    const long long flow_inf = 1e18;
    vector<FlowEdge> edges;
    vector<char> alive;
    vector<long long> pin, pout;
    vector<list<int> > in, out;
    vector<vector<int> > adj;
    vector<long long> ex;
    int n, m = 0;
    int s, t;
    vector<int> level;
    vector<int> q;
    int qh, qt;
    void resize(int _n){
        n = _n;
        ex.resize(n);
        q.resize(n);
        pin.resize(n);
        pout.resize(n);
        adj.resize(n);
        level.resize(n);
        in.resize(n);
        out.resize(n);
    }
    MPM(){}
    MPM(int _n, int _s, int _t){resize(_n); s = _s; t = _t;}
    void add_edge(int v, int u, long long cap){
        edges.push_back(FlowEdge(v, u, cap));
        edges.push_back(FlowEdge(u, v, 0));
        adj[v].push_back(m);
        adj[u].push_back(m + 1);
        m += 2;
    }
    bool bfs(){
        while(qh < qt){
            int v = q[qh++];
            for(int id : adj[v]){
                if(edges[id].cap - edges[id].flow < 1)continue;
                if(level[edges[id].u] != -1)continue;
                level[edges[id].u] = level[v] + 1;
                q[qt++] = edges[id].u;
            }
        }
        return level[t] != -1;
    }
    long long pot(int v){
        return min(pin[v], pout[v]);
    }
    void remove_node(int v){
        for(int i : in[v]){
            int u = edges[i].v;
            auto it = find(out[u].begin(), out[u].end(), i);
            out[u].erase(it);
            pout[u] -= edges[i].cap - edges[i].flow;
        }
        for(int i : out[v]){
            int u = edges[i].u;
            auto it = find(in[u].begin(), in[u].end(), i);
            in[u].erase(it);
            pin[u] -= edges[i].cap - edges[i].flow;
        }
    }
    void push(int from, int to, long long f, bool forw){
        qh = qt = 0;
        ex.assign(n, 0);
        ex[from] = f;
        q[qt++] = from;
        while(qh < qt){
            int v = q[qh++];
            if(v == to)
                break;
            long long must = ex[v];
            auto it = forw ? out[v].begin() : in[v].begin();
            while(true){
                int u = forw ? edges[*it].u : edges[*it].v;
                long long pushed = min(must, edges[*it].cap - edges[*it].flow);
                if(pushed == 0)break;
                if(forw){
                    pout[v] -= pushed;
                    pin[u] -= pushed;
                }
                else{
                    pin[v] -= pushed;
                    pout[u] -= pushed;
                }
                if(ex[u] == 0)
                    q[qt++] = u;
                ex[u] += pushed;
                edges[*it].flow += pushed;
                edges[(*it)^1].flow -= pushed;
                must -= pushed;
                if(edges[*it].cap - edges[*it].flow == 0){
                    auto jt = it;
                    ++jt;
                    if(forw){
                        in[u].erase(find(in[u].begin(), in[u].end(), *it));
                        out[v].erase(it);
                    }
                    else{
                        out[u].erase(find(out[u].begin(), out[u].end(), *it));
                        in[v].erase(it);
                    }
                    it = jt;
                }
                else break;
                if(!must)break;
            }
        }
    }
    long long flow(){
        long long ans = 0;
        while(true){
            pin.assign(n, 0);
            pout.assign(n, 0);
            level.assign(n, -1);
            alive.assign(n, true);
            level[s] = 0;
            qh = 0; qt = 1;
            q[0] = s;
            if(!bfs())
                break;
            for(int i = 0; i < n; i++){
                out[i].clear();
                in[i].clear();
            }
            for(int i = 0; i < m; i++){
                if(edges[i].cap - edges[i].flow == 0)
                    continue;
                int v = edges[i].v, u = edges[i].u;
                if(level[v] + 1 == level[u] && (level[u] < level[t] || u == t)){
                    in[u].push_back(i);
                    out[v].push_back(i);
                    pin[u] += edges[i].cap - edges[i].flow;
                    pout[v] += edges[i].cap - edges[i].flow;
                }
            }
            pin[s] = pout[t] = flow_inf;
            while(true){
                int v = -1;
                for(int i = 0; i < n; i++){
                    if(!alive[i])continue;
                    if(v == -1 || pot(i) < pot(v))
                        v = i;
                }
                if(v == -1)
                    break;
                if(pot(v) == 0){
                    alive[v] = false;
                    remove_node(v);
                    continue;
                }
                long long f = pot(v);
                ans += f;
                push(v, s, f, false);
                push(v, t, f, true);
                alive[v] = false;
                remove_node(v);
            }
        }
        return ans;
    }
};