# Maximum flow - MPM algorithm

MPM (Malhotra, Pramodh-Kumar and Maheshwari) algorithm solves the maximum flow problem in $O(V^3)$. This algorithm is similar to Dinic's algorithm.

## Algorithm

Like Dinic's algorithm, MPM runs in phases, during each phase we find the blocking flow in the layered network of the residual network of $G$. The main difference from Dinic's is how we find the blocking flow. Consider the layered network $L$. For each node we define its' inner potential and outer potential as:

\begin{align} p_{in}(v) &= \sum\limits_{(u, v)\in L}(c(u, v) - f(u, v)) \\ p_{out}(v) &= \sum\limits_{(v, u)\in L}(c(v, u) - f(v, u)) \end{align}

Also we set $p_{in}(s) = p_{out}(t) = \infty$. Given $p_{in}$ and $p_{out}$ we define the potential as $p(v) = min(p_{in}(v), p_{out}(v))$. We call a node $r$ a reference node if $p(r) = min\{p(v)\}$. Consider a reference node $r$. We claim that the flow can be increased by $p(r)$ in such a way that $p(r)$ becomes $0$. It is true because $L$ is acyclic, so we can push the flow out of $r$ by outgoing edges and it will reach $t$ because each node has enough outer potential to push the flow out when it reaches it. Similarly, we can pull the flow from $s$. The construction of the blocked flow is based on this fact. On each iteration we find a reference node and push the flow from $s$ to $t$ through $r$. This process can be simulated by BFS. All completely saturated arcs can be deleted from $L$ as they won't be used later in this phase anyway. Likewise, all the nodes different from $s$ and $t$ without outgoing or incoming arcs can be deleted.

Each phase works in $O(V^2)$ because there are at most $V$ iterations (because at least the chosen reference node is deleted), and on each iteration we delete all the edges we passed through except at most $V$. Summing, we get $O(V^2 + E) = O(V^2)$. Since there are less than $V$ phases (see the proof here), MPM works in $O(V^3)$ total.

## Implementation

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<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);
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));
m += 2;
}
bool bfs(){
while(qh < qt){
int v = q[qh++];
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;
}
};