Processing math: 100%
My Algorithm : kopricky アルゴリズムライブラリ

kopricky アルゴリズムライブラリ

Push Relabel

コードについての説明(個人的メモ)

最大流を求めるアルゴリズム.
最大流の値を高速に求めるために最終的な結果が依然として正しい flow となっていない可能性がある(excess を全て解消していない). なので valid な flow routing を求めたい場合は valiate() のコメントアウトを外す.
Stack(LIFO) を用いるかコメントアウトしている Queue(FIFO) を用いるかでケースによっては多少速度が変わることがある.
元論文は "A New Approach to the Maximum-Flow Problem" [Goldberg and Tarjan 1988]".

(追記)
容量が整数である場合に超過量を 2 べきでスケーリングしていく excess scaling 版の PushRelabel も 2 つ目の実装として置いた. 時間計算量は初期超過量の最大値 (maxe∈δ+(s)e.cap) を U としたとき O(nm+n2⌈log(U+1)⌉) となる. ケースによってはこちらの方が速い.

(関数)
add_edge(from,to,cap) : 容量 cap の辺 (from,to) を追加する
solve(s,t) : s, t 間の最大フローを計算する

時間計算量: O(n2m1/2) (ただしこれは m=Ω(n) かつ m=O(n2) の場合の計算量)

コード

  1. class Stack {
  2. private:
  3.     const int N, H;
  4.     vector<int> node;
  5. public:
  6.     Stack(const int _N, const int _H) : N(_N), H(_H), node(N+H){ clear(); }
  7.     bool empty(const int h) const { return node[N+h] == N+h; }
  8.     int top(const int h) const { return node[N+h]; }
  9.     void pop(const int h){ node[N+h] = node[node[N+h]]; }
  10.     void push(const int h, const int u){ node[u] = node[N+h], node[N+h] = u; }
  11.     void clear(){ iota(node.begin() + N, node.end(), N); }
  12. };
  13.  
  14. // class Queue {
  15. // private:
  16. //     const int N, H;
  17. //     vector<int> node, last;
  18. // public:
  19. //     Queue(const int _N, const int _H) : N(_N), H(_H), node(N+H), last(H){ clear(); }
  20. //     bool empty(const int h) const { return node[N+h] == N+h; }
  21. //     int top(const int h) const { return node[N+h]; }
  22. //     void pop(const int h){ node[N+h] = node[node[N+h]]; if(empty(h)) last[h] = N+h; }
  23. //     void push(const int h, const int u){ node[node[last[h]] = u] = N+h, last[h] = u; }
  24. //     void clear(){ iota(node.begin() + N, node.end(), N); iota(last.begin(), last.end(), N); }
  25. // };
  26.  
  27. class List {
  28. public:
  29.     struct node {
  30.         int prev, next;
  31.     };
  32.     const int N, H;
  33.     vector<node> dat;
  34.     List(const int _N, const int _H) : N(_N), H(_H), dat(N+H){ clear(); }
  35.     bool empty(const int h) const { return (dat[N+h].next == N+h); }
  36.     bool more_one(const int h) const { return dat[N+h].prev != dat[N+h].next; }
  37.     void insert(const int h, const int u){
  38.         dat[u].prev = dat[N+h].prev, dat[u].next = N+h;
  39.         dat[dat[N+h].prev].next = u, dat[N+h].prev = u;
  40.     }
  41.     void erase(const int u){
  42.         dat[dat[u].prev].next = dat[u].next, dat[dat[u].next].prev = dat[u].prev;
  43.     }
  44.     void clear(){
  45.         for(int i = N; i < N+H; ++i) dat[i].prev = dat[i].next = i;
  46.     }
  47. };
  48.  
  49. template <typename T> class PushRelabel {
  50. public:
  51.     struct edge {
  52.         const int to, rev;
  53.         T cap;
  54.         edge(const int _to, const int _rev, const T _cap) : to(_to), rev(_rev), cap(_cap){}
  55.     };
  56. private:
  57.     const int V;
  58.     int s, t, pot_max, checker;
  59.     vector<T> excess;
  60.     vector<int> potential, cur_edge, que;
  61.     List all_ver;
  62.     Stack act_ver;
  63.     // Queue act_ver;
  64.     int calc_active(){
  65.         pot_max = -1;
  66.         for(int i = 0; i < V; ++i){
  67.             if(potential[i] < V){
  68.                 cur_edge[i] = 0;
  69.                 pot_max = max(potential[i], pot_max);
  70.                 all_ver.insert(potential[i], i);
  71.                 if(excess[i] > 0 && i != t) act_ver.push(potential[i], i);
  72.             }else{
  73.                 potential[i] = V+1;
  74.             }
  75.         }
  76.         return pot_max;
  77.     }
  78.     void bfs(){
  79.         for(int i = 0; i < V; ++i) potential[i] = max(potential[i], V);
  80.         potential[t] = 0;
  81.         int qh = 0, qt = 0;
  82.         for(que[qt++] = t; qh++ < qt;){
  83.             int u = que[qh - 1];
  84.             for(const edge& e : G[u]){
  85.                 if(potential[e.to] == V && G[e.to][e.rev].cap > 0){
  86.                     potential[e.to] = potential[u] + 1, que[qt++] = e.to;
  87.                 }
  88.             }
  89.         }
  90.     }
  91.     int init(){
  92.         potential[s] = V + 1;
  93.         bfs();
  94.         for(edge& e : G[s]){
  95.             if(potential[e.to] < V){
  96.                 G[e.to][e.rev].cap = e.cap, excess[s] -= e.cap, excess[e.to] += e.cap;
  97.                 e.cap = 0;
  98.             }
  99.         }
  100.         return calc_active();
  101.     }
  102.     int global_relabel(){
  103.         bfs();
  104.         all_ver.clear(), act_ver.clear();
  105.         return calc_active();
  106.     }
  107.     void gap_relabel(const int u){
  108.         for(int i = potential[u]; i <= pot_max; ++i){
  109.             for(int id = all_ver.dat[V+i].next; id < V; id = all_ver.dat[id].next){
  110.                 potential[id] = V + 1;
  111.             }
  112.             all_ver.dat[V+i].prev = all_ver.dat[V+i].next = V + i;
  113.         }
  114.     }
  115.     int discharge(const int u){
  116.         for(int& i = cur_edge[u]; i < (int)G[u].size(); ++i){
  117.             edge& e = G[u][i];
  118.             if(potential[u] == potential[e.to] + 1 && e.cap > 0){
  119.                 if(push(u, e)) return potential[u];
  120.             }
  121.         }
  122.         return relabel(u);
  123.     }
  124.     bool push(const int u, edge& e){
  125.         const T f = min(e.cap, excess[u]);
  126.         const int v = e.to;
  127.         e.cap -= f, excess[u] -= f;
  128.         G[v][e.rev].cap += f, excess[v] += f;
  129.         if(excess[v] == f && v != t) act_ver.push(potential[v], v);
  130.         return (excess[u] == 0);
  131.     }
  132.     int relabel(const int u){
  133.         ++checker;
  134.         int prv = potential[u], cur = V;
  135.         for(int i = 0; i < (int)G[u].size(); ++i){
  136.             const edge& e = G[u][i];
  137.             if(cur > potential[e.to] + 1 && e.cap > 0){
  138.                 cur_edge[u] = i;
  139.                 cur = potential[e.to] + 1;
  140.             }
  141.         }
  142.         if(all_ver.more_one(prv)){
  143.             all_ver.erase(u);
  144.             if((potential[u] = cur) == V) return potential[u] = V+1, prv;
  145.             act_ver.push(cur, u);
  146.             all_ver.insert(cur, u);
  147.             pot_max = max(pot_max, cur);
  148.         }else{
  149.             gap_relabel(u);
  150.             return pot_max = prv - 1;
  151.         }
  152.         return cur;
  153.     }
  154.     int validate_discharge(const int u){
  155.         for(int& i = cur_edge[u]; i < (int)G[u].size(); ++i){
  156.             edge& e = G[u][i];
  157.             if(potential[u] == potential[e.to] + 1 && e.cap > 0){
  158.                 if(validate_push(u, e)) return potential[u] - V;
  159.             }
  160.         }
  161.         return validate_relabel(u);
  162.     }
  163.     bool validate_push(const int u, edge& e){
  164.         T f = min(e.cap, excess[u]);
  165.         const int v = e.to;
  166.         e.cap -= f, excess[u] -= f;
  167.         G[v][e.rev].cap += f, excess[v] += f;
  168.         if(excess[v] == f) act_ver.push(potential[v] - V, v);
  169.         return (excess[u] == 0);
  170.     }
  171.     int validate_relabel(const int u){
  172.         int cur = 2 * V;
  173.         for(int i = 0; i < (int)G[u].size(); ++i){
  174.             const edge& e = G[u][i];
  175.             if(cur > potential[e.to] + 1 && e.cap > 0){
  176.                 cur_edge[u] = i;
  177.                 cur = potential[e.to] + 1;
  178.             }
  179.         }
  180.         potential[u] = cur;
  181.         act_ver.push(cur - V, u);
  182.         return cur - V;
  183.     }
  184.     void validate(){
  185.         for(int i = 0; i < V; ++i){
  186.             if(potential[i] >= V){
  187.                 cur_edge[i] = 0, potential[i] = V + 1;
  188.                 if(excess[i] > 0) act_ver.push(1, i);
  189.             }
  190.         }
  191.         potential[s] = V;
  192.         int level = 0;
  193.         while(level >= 0){
  194.             if(act_ver.empty(level)){
  195.                 --level;
  196.                 continue;
  197.             }
  198.             int u = act_ver.top(level);
  199.             act_ver.pop(level);
  200.             level = validate_discharge(u);
  201.         }
  202.     }
  203. public:
  204.     vector<vector<edge> > G;
  205.     PushRelabel(const int node_size)
  206.         : V(node_size), pot_max(-1), checker(0), excess(V, (T)0),
  207.             potential(V, 0), cur_edge(V), que(V), all_ver(V, V), act_ver(V, V), G(V){}
  208.     void add_edge(const int _from, const int _to, const T _cap){
  209.         G[_from].emplace_back(_to, (int)G[_to].size(), _cap);
  210.         G[_to].emplace_back(_from, (int)G[_from].size() - 1, 0);
  211.     }
  212.     T solve(const int source, const int sink){
  213.         s = source, t = sink;
  214.         int level = init();
  215.         while(level >= 0){
  216.             if(act_ver.empty(level)){
  217.                 --level;
  218.                 continue;
  219.             }
  220.             int u = act_ver.top(level);
  221.             act_ver.pop(level);
  222.             level = discharge(u);
  223.             if(checker >= V / 2){
  224.                 level = global_relabel();
  225.                 checker = 0;
  226.             }
  227.         }
  228.         // validate();
  229.         return excess[t];
  230.     }
  231. };

コード(excess scaling 版)

  1. class Stack {
  2. private:
  3.     const int N, H;
  4.     int sz;
  5.     vector<int> node;
  6. public:
  7.     Stack(const int _N, const int _H) : N(_N), H(_H), node(N+H){ clear(); }
  8.     bool empty() const { return sz == 0; }
  9.     bool empty(const int h) const { return node[N+h] == N+h; }
  10.     int top(const int h) const { return node[N+h]; }
  11.     void pop(const int h){ --sz, node[N+h] = node[node[N+h]]; }
  12.     void push(const int h, const int u){ ++sz, node[u] = node[N+h], node[N+h] = u; }
  13.     void clear(){ sz = 0, iota(node.begin() + N, node.end(), N); }
  14. };
  15.  
  16. // class Queue {
  17. // private:
  18. //     const int N, H;
  19. // 	int sz;
  20. //     vector<int> node, last;
  21. // public:
  22. //     Queue(const int _N, const int _H) : N(_N), H(_H), node(N+H), last(H){ clear(); }
  23. //     bool empty() const { return sz == 0; }
  24. //     bool empty(const int h) const { return node[N+h] == N+h; }
  25. //     int top(const int h) const { return node[N+h]; }
  26. //     void pop(const int h){ --sz, node[N+h] = node[node[N+h]]; if(empty(h)) last[h] = N+h; }
  27. //     void push(const int h, const int u){ ++sz, node[node[last[h]] = u] = N+h, last[h] = u; }
  28. //     void clear(){ sz = 0, iota(node.begin() + N, node.end(), N), iota(last.begin(), last.end(), N); }
  29. // };
  30.  
  31. template <typename T> class PushRelabel {
  32. public:
  33.     static_assert(std::is_integral<T>::value, "Integral required.");
  34.     struct edge {
  35.         const int to, rev;
  36.         T cap;
  37.         edge(const int _to, const int _rev, const T _cap) : to(_to), rev(_rev), cap(_cap){}
  38.     };
  39. private:
  40.     const int V;
  41.     int s, t, checker;
  42.     vector<T> excess;
  43.     vector<int> potential, cur_edge, que;
  44.     Stack act_ver;
  45.     // Queue act_ver;
  46.     static unsigned long long ceil2(unsigned long long v){
  47.         --v;
  48.         v = v | (v >> 1), v = v | (v >> 2);
  49.         v = v | (v >> 4), v = v | (v >> 8);
  50.         v = v | (v >> 16), v = v | (v >> 32);
  51.         return ++v;
  52.     }
  53.     int calc_active(const T delta){
  54.         int pot_min = V;
  55.         for(int i = 0; i < V; ++i){
  56.             if(potential[i] < V){
  57.                 if(excess[i] >= delta && i != t){
  58.                     act_ver.push(potential[i], i);
  59.                     pot_min = min(pot_min, potential[i]);
  60.                 }
  61.             }else{
  62.                 potential[i] = V + 1;
  63.             }
  64.         }
  65.         return pot_min;
  66.     }
  67.     void bfs(){
  68.         for(int i = 0; i < V; ++i) potential[i] = max(potential[i], V);
  69.         potential[t] = 0;
  70.         int qh = 0, qt = 0;
  71.         for(que[qt++] = t; qh++ < qt;){
  72.             int u = que[qh - 1];
  73.             for(const edge& e : G[u]){
  74.                 if(potential[e.to] == V && G[e.to][e.rev].cap > 0){
  75.                     potential[e.to] = potential[u] + 1, que[qt++] = e.to;
  76.                 }
  77.             }
  78.         }
  79.     }
  80.     T init(){
  81.         T mx = 0;
  82.         potential[s] = V + 1;
  83.         bfs();
  84.         for(edge& e : G[s]){
  85.             if(potential[e.to] < V){
  86.                 G[e.to][e.rev].cap = e.cap, excess[s] -= e.cap, excess[e.to] += e.cap;
  87.                 mx = max(mx, e.cap), e.cap = 0;
  88.             }
  89.         }
  90.         return mx;
  91.     }
  92.     int global_relabel(const T delta){
  93.         bfs();
  94.         act_ver.clear();
  95.         return calc_active(delta);
  96.     }
  97.     int discharge(const int u, const T delta){
  98.         for(int& i = cur_edge[u]; i < (int)G[u].size(); ++i){
  99.             edge& e = G[u][i];
  100.             if(potential[u] == potential[e.to] + 1 && e.cap > 0){
  101.                 return push(u, e, delta) ? potential[e.to] : potential[u];
  102.             }
  103.         }
  104.         return relabel(u);
  105.     }
  106.     bool push(const int u, edge& e, const T delta){
  107.         const int v = e.to;
  108.         T f = min(e.cap, excess[u]);
  109.         if(v != t) f = min(f, 2 * delta - 1 - excess[v]);
  110.         e.cap -= f, excess[u] -= f;
  111.         G[v][e.rev].cap += f, excess[v] += f;
  112.         if(excess[u] >= delta) act_ver.push(potential[u], u);
  113.         if(excess[v] >= delta && v != t){
  114.             act_ver.push(potential[v], v);
  115.             return true;
  116.         }else{
  117.             return false;
  118.         }
  119.     }
  120.     int relabel(const int u){
  121.         ++checker;
  122.         int prv = potential[u], cur = V;
  123.         for(int i = 0; i < (int)G[u].size(); ++i){
  124.             const edge& e = G[u][i];
  125.             if(cur > potential[e.to] + 1 && e.cap > 0){
  126.                 cur_edge[u] = i, cur = potential[e.to] + 1;
  127.             }
  128.         }
  129.         if((potential[u] = cur) == V) return potential[u] = V + 1, prv;
  130.         act_ver.push(cur, u);
  131.         return prv;
  132.     }
  133.     int validate_calc_active(const T delta){
  134.         int pot_min = 2 * V;
  135.         for(int i = 0; i < V; ++i){
  136.             if(excess[i] >= delta && i != t){
  137.                 act_ver.push(potential[i] - V, i);
  138.                 pot_min = min(pot_min, potential[i]);
  139.             }
  140.         }
  141.         return pot_min;
  142.     }
  143.     int validate_discharge(const int u, const T delta){
  144.         for(int& i = cur_edge[u]; i < (int)G[u].size(); ++i){
  145.             edge& e = G[u][i];
  146.             if(potential[u] == potential[e.to] + 1 && e.cap > 0){
  147.                 return validate_push(u, e, delta) ? potential[e.to] : potential[u];
  148.             }
  149.         }
  150.         return validate_relabel(u);
  151.     }
  152.     bool validate_push(const int u, edge& e, const T delta){
  153.         const int v = e.to;
  154.         T f = min({e.cap, excess[u], 2 * delta - 1 - excess[v]});
  155.         e.cap -= f, excess[u] -= f;
  156.         G[v][e.rev].cap += f, excess[v] += f;
  157.         if(excess[u] >= delta) act_ver.push(potential[u] - V, u);
  158.         if(excess[v] >= delta){
  159.             act_ver.push(potential[v] - V, v);
  160.             return true;
  161.         }else{
  162.             return false;
  163.         }
  164.     }
  165.     int validate_relabel(const int u){
  166.         int prv = potential[u], cur = 2 * V;
  167.         for(int i = 0; i < (int)G[u].size(); ++i){
  168.             const edge& e = G[u][i];
  169.             if(cur > potential[e.to] + 1 && e.cap > 0){
  170.                 cur_edge[u] = i, cur = potential[e.to] + 1;
  171.             }
  172.         }
  173.         potential[u] = cur, act_ver.push(cur - V, u);
  174.         return prv;
  175.     }
  176.     void validate(){
  177.         T mx = 0;
  178.         for(int i = 0; i < V; ++i){
  179.             if(potential[i] >= V){
  180.                 cur_edge[i] = 0, potential[i] = V + 1;
  181.                 mx = max(mx, excess[i]);
  182.             }
  183.         }
  184.         potential[s] = V;
  185.         T delta = (ceil2(mx + 1) >> 1);
  186.         int level = validate_calc_active(delta) - V;
  187.         while(delta > 0){
  188.             while(!act_ver.empty()){
  189.                 if(act_ver.empty(level)){
  190.                     ++level;
  191.                     continue;
  192.                 }
  193.                 const int u = act_ver.top(level);
  194.                 act_ver.pop(level);
  195.                 level = validate_discharge(u, delta) - V;
  196.             }
  197.             if(delta == 1) break;
  198.             delta >>= 1, level = validate_calc_active(delta) - V;
  199.         }
  200.     }
  201. public:
  202.     vector<vector<edge> > G;
  203.     PushRelabel(const int node_size)
  204.         : V(node_size), checker(0), excess(V, (T)0),
  205.             potential(V, 0), cur_edge(V), que(V), act_ver(V, V), G(V){}
  206.     void add_edge(const int _from, const int _to, const T _cap){
  207.         G[_from].emplace_back(_to, (int)G[_to].size(), _cap);
  208.         G[_to].emplace_back(_from, (int)G[_from].size() - 1, 0);
  209.     }
  210.     T solve(const int source, const int sink){
  211.         s = source, t = sink;
  212.         T delta = (ceil2(init() + 1) >> 1);
  213.         int level = calc_active(delta);
  214.         while(delta > 0){
  215.             while(!act_ver.empty()){
  216.                 if(act_ver.empty(level)){
  217.                     ++level;
  218.                     continue;
  219.                 }
  220.                 const int u = act_ver.top(level);
  221.                 act_ver.pop(level);
  222.                 level = discharge(u, delta);
  223.                 if(checker >= V / 2){
  224.                     level = global_relabel(delta);
  225.                     checker = 0;
  226.                 }
  227.             }
  228.             if(delta == 1) break;
  229.             delta >>= 1, level = calc_active(delta);
  230.         }
  231.         // validate();
  232.         return excess[t];
  233.     }
  234. };

verify 用の問題

AOJ : Maximum Flow 提出コード 提出コード(excess scaling 版)