なぜ，

最小費用流 @ C++ – 落書き、時々落学

をしたかというと，仇打ちです．

当初，最小費用循環流を実装していた．

→できた

→検証しよう

→UVa 10594 時間オーバー，別の問題（2部グラフの最小費用最大マッチング）ではOKだった．

→別のアルゴリズムで UVa 10594を解く．

という流れ．はじめは，速度とか気にしてなかったのが原因か．

で本題．out-of-kilterは最小費用循環流を解くアルゴリズム．

といっても，最小費用流と最小費用循環流は簡単に変換できるから，最小費用流のアルゴリズムといっても良い．

とりあえず私の中のout-of-kilterのイメージは相補性条件が成り立つように，

負閉路消したり，ポテンシャル調整したりしていく．

なかなか，なんでこれで最適解が求まるか，傍から見ると不思議だったりしたアルゴリズム．

以前から実装してみたかったアルゴリズム．

そんなに複雑なことをしないので，実装はそこまで難しくないのだが，関数内関数が使えないのに，使おうとして，

酷いマクロを召喚している． 

#include "graph.hpp"
#include <stack>
#include <queue>

#define RC(_t, _h) (c(_t, _h) + p[_t] - p[_h])  // reduced cost
#define R(_t, _h, _f) \
  (_f ?                                                                 \
   RC(_t, _h) > 0 ? l(_t, _h) - f(_t, _h): u(_t, _h) - f(_t, _h):       \
   RC(_t, _h) < 0 ? f(_t, _h) - u(_t, _h): f(_t, _h) - l(_t, _h)        \
   )
#define EXIST_EDGE(_t, _h, _f) \
  (_f ?                                                                 \
   (RC(_t, _h) > 0 && f(_t, _h) < l(_t, _h)) || (RC(_t, _h) <= 0 && f(_t, _h) < u(_t, _h)): \
   (RC(_t, _h) < 0 && f(_t, _h) > u(_t, _h)) || (RC(_t, _h) >= 0 && f(_t, _h) > l(_t, _h)) \
   )
#define IN_KILTER(_t, _h, _f) \
  (_f ?                                                                \
   RC(_t, _h) == 0 && l(_t, _h) <= f(_t, _h) && f(_t, _h) <  u(_t, _h): \
   RC(_t, _h) == 0 && l(_t, _h) <  f(_t, _h) && f(_t, _h) <= u(_t, _h)  \
   )
#define OUT_OF_KILTER(_t, _h, _f) (EXIST_EDGE(_t, _h, _f) && !IN_KILTER(_t, _h, _f))
#define OUT_OF_KILTER_(_x) OUT_OF_KILTER(_x.first.first, _x.first.second, _x.second)
#define FOR(x, xs) for (__typeof( (xs).begin() ) x = (xs).begin(); x != (xs).end(); x++)
#define ABS(x) ((x) >= 0? (x): -(x))

template <class T>
bool
circulation(const Graph &g,  // 単純有向グラフ
            const EdgeProperty<T> &c,  // 費用
            const EdgeProperty<T> &l,  // 下限容量
            const EdgeProperty<T> &u,  // 上限容量
            EdgeProperty<T> f  // 最小費用循環流が記録される
            ) {
  const int n = g.succ.size();
  
  T *p = new T[n]; fill_n(p, n, 0);  // initialize potential
  for (Vertex v = 0; v < n; ++v) FOR (w, g.succ[v]) f(v, *w) = 0;  // initialize flow
  typedef pair<Edge,bool> EB;
  stack<EB> out;  // edges in out-of-kilter
  for (Vertex v = 0; v < n; ++v) {  // intialize out-of-kilter
    FOR (w, g.succ[v]) {
      if (OUT_OF_KILTER(v, *w, true )) out.push(EB(Edge(v, *w), true));
      if (OUT_OF_KILTER(v, *w, false)) out.push(EB(Edge(v, *w), false));
    }
  }
  
  while (!out.empty()) {
    
    // choose an edge in out-of-kilter
    Edge e = out.top().first;
    bool forward = out.top().second;
    if(!OUT_OF_KILTER(e.first, e.second, forward)) {
        out.pop();
        continue;
      }
    
    // find circuit which contains edge a, BFS from head to tail
    bool *fw = new bool[n];
    T d = R(e.first, e.second, forward);
    Vertex t, h; // tail and head of Edge e. that is, tail -> head.
    if (forward) t = e.first, h = e.second;
    else t = e.second, h = e.first;
    queue<Vertex> q;
    Vertex* prev = new Vertex[n]; fill_n(prev, n, -1), prev[h] = h;
    q.push(h);
    
    while (!q.empty() && prev[t] < 0) {
      Vertex v = q.front(); q.pop();
      FOR (w, g.succ[v]) if (EXIST_EDGE(v, *w, true))
        q.push(*w), prev[*w] = v, fw[*w] = true, d = min(d, R(v, *w, fw[*w]));
      FOR (w, g.pred[v]) if (EXIST_EDGE(*w, v, false))
        q.push(*w), prev[*w] = v, fw[*w] = false, d = min(d, R(*w, v, fw[*w]));
    }
    
    if (prev[t] >= 0) {  // circuit is found
      for (Vertex v = t; v != prev[v]; v = prev[v]) // update flow
        if (fw[v]) f(prev[v], v) += d;
        else       f(v, prev[v]) -= d;
      f(e) += forward? d: -d;
      continue;
    }
    
    // v is reachable from h <=> prev[v] >= 0
    d = INF;
    stack<EB> stack;  // candidate edge of out-of-kilter
    for (Vertex v = 0; v < n; ++v)
      if (prev[v] >= 0) {
        FOR (w, g.succ[v]) if (RC(v, *w) > 0 && prev[*w] <  0 && f(v, *w) <= u(v, *w))
            d = min(d, RC(v, *w)), stack.push(EB(Edge(v, *w), false));
      } else {
        FOR (w, g.succ[v]) if (RC(v, *w) < 0 && prev[*w] >= 0 && l(v, *w) <= f(v, *w))
          d = min(d, -RC(v, *w)), stack.push(EB(Edge(v, *w), true));
      }
    if (stack.empty()) return false; // there is no feasible circulation
    // update potential
    for (Vertex v = 0; v < n; ++v)
      if (prev[v] < 0) p[v] += d;
    // update out-of-kilter edges
    while(!stack.empty()) {
      if(OUT_OF_KILTER_(stack.top()))
        out.push(stack.top());
      stack.pop();
    }
  }
  return true;
}


int main() {
  Graph g(2);
  Edge e[2] = {g.add_edge(0, 1), g.add_edge(1, 0)};
  EdgeProperty<int> c, l, u, f;
  c(e[0]) = 10;
  l(e[0]) = 0;
  u(e[0]) = 5;
  c(e[1]) = 10;
  l(e[1]) = 3;
  u(e[1]) = 5;
  circulation<int>(g, c, l, u, f);
  return 0;
}

このままだと，計算量が費用関数に依存するがスケーリングという手法を用いると

費用関数に依存しない，強多項式時間のアルゴリズムが得られる．

らしい．雰囲気は分かるんだけど，細かいところが良く分からん．