最大流 (Goldberg-Tarjan)

説明

Goldberg Tarjan は最大流問題を解くアルゴリズムである．Edmond Karp が常に容量制限を満たすようにして流量を増加させていくアルゴリズムであるのに対し，Goldberg Tarjan は流量を最大に保った上で容量制限を満たすまで更新していく．匙で擦りきるようなイメージを持つと良い．

具体的には，アクティブな頂点に対して，PUSH と RELABEL のうち実行できるものを適当に実行していく，というアルゴリズムとなる．ただし PUSH はその頂点から出る枝に流せる限り流す関数，RELABEL は自分の高さを周りの高さよりも一つ高くする関数である．どの頂点を選ぶかによって実行速度が多少変化するが，理論的にはたいてい augmentation 系のアルゴリズムよりも高速になる．

アクティブな頂点の選び方として，高さが最高のものを常に取る (Highest-Level Selection) ことにすると，計算量は O(V^2 \sqrt{m}) となる．実用的にはここにいくつかのヒューリスティクスを突っ込むことで，相当に高速なアルゴリズムとなる．代表的なヒューリスティクスとして GLOBAL-RELABELING (定期的に最短路をまじめに計算する) と GAP-RELABELING (ギャップがあった場合上側を全消去する) がある．以下では GLOBAL-RELABELING のみを実装している．

ちなみに，辺数などにも依存するが，この実装ならば Dinic のアルゴリズムがこれと同程度以上の性能をたたき出す．好みで使い分けると良い．

計算量

O(V^2 \sqrt{m})．

ToDo

Runtime Error!!

GAP-RELABELING の導入

ソースコード

#define RESIDUE(s,t) (capacity[s][t]-flow[s][t])

#define GLOBAL_RELABELING() { \
  queue<int> Q; Q.push(t); \
  fill(ALL(d), INF); d[t] = 0; \
  while (!Q.empty()) { \
    int u = Q.front(); Q.pop(); \
    FOR(e, g[u]) if (RESIDUE(e->dst, u) > 0 && d[u] + 1 < d[e->dst])  \
      Q.push(e->dst), d[e->dst] = d[u] + 1; \
  } \
}
#define PUSH(u, v) { \
  Weight delta = min(excess[u], RESIDUE(u, v)); \
  flow[u][v] += delta; flow[v][w] -= delta; \
  excess[u] -= delta; excess[v] += delta; }
Weight maximumFlow(const Graph &g, int s, int t) {
  int n = g.size(), count = 0;
  Matrix flow(n, Array(n)), capacity(n, Array(n)); // adj. matrix
  REP(u,n) FOR(e,g[u]) capacity[e->src][e->dst] += e->weight;

  vector<Weight> excess(n); excess[s] = INF; // initialize step
  vector<int> d(n);
  GLOBAL_RELABELING();
  vector< queue<int> > B(n); B[ d[s] ].push( s );

  for (int b = d[s]; b >= 0; ) {
    if (B[b].empty()) { --b; continue; }
    int v = B[b].front(); B[b].pop();
    if (excess[v] == 0 || v == t) continue;

    FOR(e, g[v]) {
      int w = e->dst; // e is the current edge of v
      if (RESIDUE(v,w) > 0 && d[v] == d[w] + 1) { // (w,v) is admissible
        PUSH(v, w);
        if (excess[w] > 0 && w != t) B[d[w]].push( w );
      }
    }
    if (excess[v] == 0) continue;
    d[v] = n;
    FOR(e, g[v]) if (RESIDUE(v, e->dst) > 0)
      d[v] = min(d[v], d[e->dst] + 1);
    if (d[v] < n) B[b = d[v]].push(v);

    if (++count % n == 0) GLOBAL_RELABELING(); // !!HEURISTICS
  }
  return excess[t];
}

Verified