AIOJ

Before reading this article, first read the definitions section of the Introduction to Network Flow wiki.

Minimum-cost flow

Given a network $G=(V,E)$ , each edge has not only a capacity constraint $c(u,v)$ but also a cost $w(u,v)$ per unit of flow.

When the flow on $(u,v)$ is $f(u,v)$ , the cost incurred is $f(u,v)\times w(u,v)$ .

$w$ also satisfies skew-symmetry, namely $w(u,v)=-w(v,u)$ .

The maximum flow with the minimum total cost in this network is called minimum-cost maximum flow; that is, minimize $\sum_{(u,v)\in E}f(u,v)\times w(u,v)$ subject to maximizing $\sum_{(s,v)\in E}f(s,v)$ .

SSP algorithm

The SSP (Successive Shortest Path) algorithm is a greedy algorithm. Its idea is to find an augmenting path with minimum unit cost and augment along it each time, until no augmenting path remains in the graph.

If the graph contains a cycle with negative unit cost, the SSP algorithm cannot correctly find the minimum-cost maximum flow of the network. In this case, a cycle-canceling algorithm must first be used to eliminate the negative cycle in the graph.

Proof

We prove the correctness of the SSP algorithm using mathematical induction and proof by contradiction.

Let $f_i$ be the minimum cost when the flow is $i$ . We assume that the initial network has no negative cycles; in this case, $f_0=0$ .

Suppose that the $f_i$ found by the SSP algorithm is a minimum cost. Starting from $f_i$ , we find a shortest augmenting path and thereby obtain $f_{i+1}$ . Then $f_{i+1}-f_i$ is the length of this shortest augmenting path.

Suppose that a smaller $f_{i+1}$ exists, and denote it by $f'_{i+1}$ . Since $f_{i+1}-f_i$ is already the shortest augmenting path, $f'_{i+1}-f_i$ must correspond to an augmenting path that passes through at least one negative cycle.

Here is the contradiction: if there is an augmenting path passing through at least one negative cycle, then $f_i$ is not a minimum cost. By adding flow to this negative cycle, we can reduce the cost corresponding to $f_i$ without increasing the flow leaving $s$ .

Therefore, the SSP algorithm can correctly find the minimum-cost maximum flow of a network without negative cycles.

Time complexity

If the Bellman–Ford algorithm is used to find shortest paths, the time complexity of finding an augmenting path once is $O(nm)$ . Let the maximum flow of the network be $f$ ; then the worst-case time complexity is $O(nmf)$ . In fact, the SSP algorithm runs in pseudo-polynomial time.

Why is the SSP algorithm pseudo-polynomial time?

The SSP algorithm has an upper bound of $O(nmf)$ on its time complexity. This is a polynomial in the value range, so it runs in pseudo-polynomial time.

One can construct a network with $m=n^2,f=2^{n/2}$ ¹ such that the SSP algorithm takes $O(n^3 2^{n/2})$ time, so the SSP algorithm is not polynomial-time.

Implementation

Simply replace the process of finding an augmenting path in the EK or Dinic algorithm with a shortest-path algorithm that finds an augmenting path with minimum unit cost.

Implementation based on the EK algorithm

struct qxx {
  int nex, t, v, c;
};

qxx e[M];
int h[N], cnt = 1;

void add_path(int f, int t, int v, int c) {
  e[++cnt] = qxx{h[f], t, v, c}, h[f] = cnt;
}

void add_flow(int f, int t, int v, int c) {
  add_path(f, t, v, c);
  add_path(t, f, 0, -c);
}

int dis[N], pre[N], incf[N];
bool vis[N];

bool spfa() {
  memset(dis, 0x3f, sizeof(dis));
  queue<int> q;
  q.push(s), dis[s] = 0, incf[s] = INF, incf[t] = 0;
  while (q.size()) {
    int u = q.front();
    q.pop();
    vis[u] = false;
    for (int i = h[u]; i; i = e[i].nex) {
      const int &v = e[i].t, &w = e[i].v, &c = e[i].c;
      if (!w || dis[v] <= dis[u] + c) continue;
      dis[v] = dis[u] + c, incf[v] = min(w, incf[u]), pre[v] = i;
      if (!vis[v]) q.push(v), vis[v] = true;
    }
  }
  return incf[t];
}

int maxflow, mincost;

void update() {
  maxflow += incf[t];
  for (int u = t; u != s; u = e[pre[u] ^ 1].t) {
    e[pre[u]].v -= incf[t], e[pre[u] ^ 1].v += incf[t];
    mincost += incf[t] * e[pre[u]].c;
  }
}

// Call: while(spfa())update();

Implementation based on the Dinic algorithm

#include <algorithm>
#include <cstdio>
#include <cstring>
#include <queue>

constexpr int N = 5e3 + 5, M = 1e5 + 5;
constexpr int INF = 0x3f3f3f3f;
int n, m, tot = 1, lnk[N], cur[N], ter[M], nxt[M], cap[M], cost[M], dis[N], ret;
bool vis[N];

void add(int u, int v, int w, int c) {
  ter[++tot] = v, nxt[tot] = lnk[u], lnk[u] = tot, cap[tot] = w, cost[tot] = c;
}

void addedge(int u, int v, int w, int c) { add(u, v, w, c), add(v, u, 0, -c); }

bool spfa(int s, int t) {
  memset(dis, 0x3f, sizeof(dis));
  memcpy(cur, lnk, sizeof(lnk));
  std::queue<int> q;
  q.push(s), dis[s] = 0, vis[s] = true;
  while (!q.empty()) {
    int u = q.front();
    q.pop(), vis[u] = false;
    for (int i = lnk[u]; i; i = nxt[i]) {
      int v = ter[i];
      if (cap[i] && dis[v] > dis[u] + cost[i]) {
        dis[v] = dis[u] + cost[i];
        if (!vis[v]) q.push(v), vis[v] = true;
      }
    }
  }
  return dis[t] != INF;
}

int dfs(int u, int t, int flow) {
  if (u == t) return flow;
  vis[u] = true;
  int ans = 0;
  for (int &i = cur[u]; i && ans < flow; i = nxt[i]) {
    int v = ter[i];
    if (!vis[v] && cap[i] && dis[v] == dis[u] + cost[i]) {
      int x = dfs(v, t, std::min(cap[i], flow - ans));
      if (x) ret += x * cost[i], cap[i] -= x, cap[i ^ 1] += x, ans += x;
    }
  }
  vis[u] = false;
  return ans;
}

int mcmf(int s, int t) {
  int ans = 0;
  while (spfa(s, t)) {
    int x;
    while ((x = dfs(s, t, INF))) ans += x;
  }
  return ans;
}

int main() {
  int s, t;
  scanf("%d%d%d%d", &n, &m, &s, &t);
  while (m--) {
    int u, v, w, c;
    scanf("%d%d%d%d", &u, &v, &w, &c);
    addedge(u, v, w, c);
  }
  int ans = mcmf(s, t);
  printf("%d %d\n", ans, ret);
  return 0;
}

Primal-Dual algorithm

Using Bellman–Ford to find shortest paths takes $O(nm)$ time, which is inferior to Dijkstra's algorithm² on both sparse and dense graphs. However, the network contains edges with negative unit cost, so Dijkstra's algorithm cannot be used directly.

The idea of the Primal-Dual algorithm is similar to that of the Johnson all-pairs shortest-path algorithm: by assigning a potential to each vertex, make the costs of all edges in the network (hereafter called edge weights) nonnegative, so that Dijkstra's algorithm can be applied to find an augmenting path with minimum unit cost.

First run a shortest-path algorithm once to find the shortest distance from the source to every vertex, which is also the initial potential $h_i$ of that vertex. Next, as in Johnson's algorithm, for an edge from $u$ to $v$ with unit cost $w$ , reset its edge weight to $w+h_u-h_v$ .

It can be seen that after setting the potentials this way, the shortest paths in the new network correspond exactly to those in the original network. The proof was given in the introduction of Johnson's algorithm and is omitted here.

Unlike an ordinary shortest-path problem, the shape of the graph changes after each augmentation, so the potentials of the vertices need to be updated.

How should they be updated? We first state the conclusion. Let the shortest distance from the source to vertex $i$ after augmentation be $d'_i$ (where the distance is computed after resetting every edge weight); it suffices to add $d'_i$ to $h_i$ . We now prove that after updating in this way, all edge weights in the graph are nonnegative.

It is easy to see that after one augmentation, because some $(i,j)$ edges lie on the augmenting path, corresponding $(j,i)$ edges are added to the residual network, and they necessarily satisfy $d'_i+(w(i,j)+h_i-h_j)=d'_j$ (otherwise the $(i,j)$ edge would not lie on the augmenting path). Rearranging gives $w(j,i)+(h_j+d'_j)-(h_i+d'_i)=0$ . Therefore, the weights of the newly added edges are nonnegative.

For an original edge, before augmentation, $d'_i+(w(i,j)+h_i-h_j) - d'_j \geq 0$ . Therefore, $w(i,j)+(d'_i+h_i)-(d'_j+h_j) \geq 0$ ; that is, using $h_i+d'_i$ as the new potential will not make the weight of $(i,j)$ negative.

Thus, after augmentation, all edge weights are nonnegative, and Dijkstra's algorithm can correctly find the shortest path in the graph.

Reference implementation

#include <algorithm>
#include <cstdio>
#include <cstring>
#include <queue>
constexpr int INF = 0x3f3f3f3f;
using namespace std;

struct edge {
  int v, f, c, next;
} e[100005];

struct node {
  int v, e;
} p[10005];

struct mypair {
  int dis, id;

  bool operator<(const mypair& a) const { return dis > a.dis; }

  mypair(int d, int x) { dis = d, id = x; }
};

int head[5005], dis[5005], vis[5005], h[5005];
int n, m, s, t, cnt = 1, maxf, minc;

void addedge(int u, int v, int f, int c) {
  e[++cnt].v = v;
  e[cnt].f = f;
  e[cnt].c = c;
  e[cnt].next = head[u];
  head[u] = cnt;
}

bool dijkstra() {
  priority_queue<mypair> q;
  for (int i = 1; i <= n; i++) dis[i] = INF;
  memset(vis, 0, sizeof(vis));
  dis[s] = 0;
  q.push(mypair(0, s));
  while (!q.empty()) {
    int u = q.top().id;
    q.pop();
    if (vis[u]) continue;
    vis[u] = 1;
    for (int i = head[u]; i; i = e[i].next) {
      int v = e[i].v, nc = e[i].c + h[u] - h[v];
      if (e[i].f && dis[v] > dis[u] + nc) {
        dis[v] = dis[u] + nc;
        p[v].v = u;
        p[v].e = i;
        if (!vis[v]) q.push(mypair(dis[v], v));
      }
    }
  }
  return dis[t] != INF;
}

void spfa() {
  queue<int> q;
  memset(h, 63, sizeof(h));
  h[s] = 0, vis[s] = 1;
  q.push(s);
  while (!q.empty()) {
    int u = q.front();
    q.pop();
    vis[u] = 0;
    for (int i = head[u]; i; i = e[i].next) {
      int v = e[i].v;
      if (e[i].f && h[v] > h[u] + e[i].c) {
        h[v] = h[u] + e[i].c;
        if (!vis[v]) {
          vis[v] = 1;
          q.push(v);
        }
      }
    }
  }
}

int main() {
  scanf("%d%d%d%d", &n, &m, &s, &t);
  for (int i = 1; i <= m; i++) {
    int u, v, f, c;
    scanf("%d%d%d%d", &u, &v, &f, &c);
    addedge(u, v, f, c);
    addedge(v, u, 0, -c);
  }
  spfa();  // Compute the initial potentials first
  while (dijkstra()) {
    int minf = INF;
    for (int i = 1; i <= n; i++) h[i] += dis[i];
    for (int i = t; i != s; i = p[i].v) minf = min(minf, e[p[i].e].f);
    for (int i = t; i != s; i = p[i].v) {
      e[p[i].e].f -= minf;
      e[p[i].e ^ 1].f += minf;
    }
    maxf += minf;
    minc += minf * h[t];
  }
  printf("%d %d\n", maxf, minc);
  return 0;
}

Exercises

References and notes

Footnotes

For the detailed construction method, see min_25's blog. ↩
On sparse graphs, heap optimization gives a time complexity of $O(m \log n)$ , while on dense graphs, omitting heap optimization gives a time complexity of $O(n^2)$ . ↩