Edmonds Karp と Push Relabel @ C++

最大流アルゴリズムを実装した．

Edmonds Karpは増大道に沿ってフローを更新（大域的操作）に対して，

Push RelabelはPush か Relabel（局所的更新）のみ．

だから，楽に実装できると思ったら，実はたいして変わらなかった．

Edmons Karp

// O(VE^2), BFS = O(E), iteration = O(VE)
template <class T>
pair< T, EdgeProperty<T> >
edmonds_karp(const Graph &g,  //多重枝無し
         const EdgeProperty<T> &c,  //capacity
         const Vertex &s,  //source
         const Vertex &t  //sink
         ) {
  const int n = g.succ.size();  // # of vertices
  T m = 0;  // maximum flow value
  EdgeProperty<T> f;  //flow
  for (Vertex u = 0; u < n; ++u) foreach (Vertex v, g.succ[u]) f[Edge(u, v)] = 0;  // Initialize flow
  while (true) {  // find augmenting path
    queue<Vertex> q;
    Vertex *p = new Vertex[n];  // previous vertex in path
    bool *r = new bool[n];  // r[i] <=> edge ?->i in path
    T d = INF;  // TODO: INF<T>
    q.push(s), fill_n(p, n, -1), p[s] = s;
    while (!q.empty() && p[t] < 0) {
      Vertex u = q.front(); q.pop();
      foreach (Vertex v, g.succ[u]) if (p[v] < 0 && f(u, v) < c(u, v))  // forward edge
        q.push(v), p[v] = u, r[v] = true, d = min(d, c(u, v) - f(u, v));
      foreach (Vertex v, g.pred[u]) if (p[v] < 0 && 0 < f(v, u))   // backward edge
        q.push(v), p[v] = u, r[v] = false, d = min(d, f(v, u));
    }
    if (p[t] < 0) return pair< T, EdgeProperty<T> >(m, f);  // No augmenting path
    for (Vertex v = t; v != p[v]; v = p[v])  // udpate flow
      if (r[v]) f[Edge(p[v], v)] += d;
      else f[Edge(v, p[v])] -= d;
    m += d;
  }
}

Push Relabel

// O(V^2E)
// find excess vertex = O(1)
// iteration (push or relabel) = O(V^2E)
template <class T>
pair< T, EdgeProperty<T> >
stack_push_relabel(const Graph &g,  //多重枝無し
       const EdgeProperty<T> &c,  //capacity
       const Vertex &s,  //source
       const Vertex &t  //sink
       ) {
  const int n = g.succ.size();  // # of vertices
  EdgeProperty<T> f;  //flow
  for (Vertex u = 0; u < n; ++u) foreach (Vertex v, g.succ[u]) f[Edge(u, v)] = 0;  // initialize preflow
  int *h = new int[n]; fill_n(h, n, 0), h[s] = n;  //  initialize height
  T *e = new T[n]; fill_n(e, n, 0);  //initialize excess
  stack<Vertex> stack;
  foreach(Vertex u, g.succ[s]) f[Edge(s, u)] = c(s, u), e[u] += f(s, u), stack.push(u);
  e[s] = INT_MAX>>1;
 LOOP: while (!stack.empty()) {
    Vertex u = stack.top(); stack.pop();
    if (e[u] == 0) continue;
    int m = 2*n;  // forall v in V, h[v] < 2*n
    foreach (Vertex v, g.succ[u])  // push forward
      if (f(u, v) < c(u, v)) {  // (u, v) in Ef
        if (h[u] > h[v]) {
          T d = min(e[u], c(u, v) - f(u, v));
          f[Edge(u, v)] += d, e[u] -= d, e[v] += d;
          if (e[u] > 0) stack.push(u);
          if (e[v] > 0 && v != t && v != s) stack.push(v);
          goto LOOP;  // reset loop
        } else { m = min (m, h[v]); }
      }
    foreach (Vertex v, g.pred[u])   // push backward
      if (0 < f(v, u)) {  // (v, u) in Ef
        if (h[u] > h[v]) {
          T d = min(e[u], f(v, u));
          f[Edge(v, u)] -= d, e[u] -= d, e[v] += d;
          if (e[u] > 0) stack.push(u);
          if (e[v] > 0 && v != t && v != s) stack.push(v);
          goto LOOP;  // reset loop
        } else { m = min (m, h[v]); }
      }
    h[u] = m + 1;  //relabel(u);
    stack.push(u);
  }
  return pair< T, EdgeProperty<T> >(e[t], f);
}

頂点から出て行く枝と頂点に入る枝を別々のリストで管理しているから，ソースが冗長になっている．

これをひとつのリストにまとめてしまえば，枝の向きを気にしなくて良くなるので，コードはもう少し簡潔になる．

しかし，有効グラフなのに枝の向きを無視するのはちょっと．

枝の容量とあわせて，向きが分かれば良いという立場もあるかもしれない．

しかし，グラフと容量は別ものだから（常に両方あるとは限らない），グラフだけでも，枝の向きが分かるほうが

精神的には落ち着く．

最後に priority queue を使って，ラベルの大きいものから処理していくコード．

多分，正しく実装できているはず（ラベル最大のものから処理しているはず）．

// hightest label
template <class T>
pair< T, EdgeProperty<T> >
queue_push_relabel(const Graph &g,  //多重枝無し
       const EdgeProperty<T> &c,  //capacity
       const Vertex &s,  //source
       const Vertex &t  //sink
       ) {
  const int n = g.succ.size();  // # of vertices
  EdgeProperty<T> f;  //flow
  for (Vertex u = 0; u < n; ++u) foreach (Vertex v, g.succ[u]) f[Edge(u, v)] = 0;  // initialize preflow
  int *h = new int[n]; fill_n(h, n, 0), h[s] = n;  //  initialize height
  T *e = new T[n]; fill_n(e, n, 0);  //initialize excess
  typedef pair<int, Vertex> Data;
  priority_queue<Data> que;
  foreach(Vertex u, g.succ[s]) f[Edge(s, u)] = c(s, u), e[u] += f(s, u), que.push(Data(0, u));
  e[s] = INT_MAX>>1;
 LOOP: while (!que.empty()) {
    Vertex u = que.top().second; que.pop();
    if (e[u] == 0) continue;
    int m = 2*n;  // forall v in V, h[v] < 2*n
    foreach (Vertex v, g.succ[u])  // push forward
      if (f(u, v) < c(u, v)) {  // (u, v) in Ef
        if (h[u] > h[v]) {
          T d = min(e[u], c(u, v) - f(u, v));
          f[Edge(u, v)] += d, e[u] -= d, e[v] += d;
          if (e[u] > 0) que.push(Data(h[u], u));
          if (e[v] > 0 && v != t && v != s) que.push(Data(h[v], v));
          goto LOOP;  // reset loop
        } else { m = min (m, h[v]); }
      }
    foreach (Vertex v, g.pred[u])   // push backward
      if (0 < f(v, u)) {  // (v, u) in Ef
        if (h[u] > h[v]) {
          T d = min(e[u], f(v, u));
          f[Edge(v, u)] -= d, e[u] -= d, e[v] += d;
          if (e[u] > 0) que.push(Data(h[u], u));
          if (e[v] > 0 && v != t && v != s) que.push(Data(h[v], v));
          goto LOOP;  // reset loop
        } else { m = min (m, h[v]); }
      }
    h[u] = m + 1;  //relabel(u);
    que.push(Data(h[u], u));
  }
  return pair< T, EdgeProperty<T> >(e[t], f);
}

# ラベル付き break ぐらいあっても良いのに．