ei1333の日記

ぺこい

最小費用流双対について 

わからない

厳密性のない議論をするのは、最高だな!厳密性のある議論を見たい人は、ごめんなさい

さいしょに

A. 最小費用流について

最小費用流って知っていますか.

ぼくは知っています. 知らなかったらごめんなさい.

アルゴリズムとかは理解してなくてもいいんですが、次の機能をもつプログラムを持っていると良いです. 持ってなかったら書いてください. 書きたくない人は書かなくていいです.

  • add_edge(u, v, cap, cost):= 頂点 u から頂点 v に最大流量 cap, コスト cost の辺を張る.
  • min_cost(S, T, F):= 頂点 S から頂点 T の流量 F の最小費用流を返す.

B. 一般化した最小費用流について

今後のために便利なので, 次のような問題を考えます(参考: 最小費用流の負辺除去 - あなたは嘘つきですかと聞かれたら「YES」と答えるブログ).

  • N 頂点 M 辺の有向グラフを考える.
  • 各頂点 v(1vN) について, 水の量 Dv が定まっている. Dv>0 のとき, 頂点 i に水が Dv あって, Dv<0 のとき水が Dv 不足していることを表している. 特に v=1NDv=0 である.
  • 各辺 e(1eM) について, 最大流量 cape と, 水を 1 流すごとにかかるコスト coste が定まっている.
  • 水の過不足を解消するための最小コストを求めよ.

この問題はさっきの A. の最小費用流を使って解くことができます. 具体的には, 次のようにすればよいです.

  • 2 つの超頂点を S,T を生やす.
  • v(1vN) ついて, Dv>0 のとき add_edge(S, v, Dv, 0), Dv<0 のとき add_edge(S, v, Dv, 0) をする.
  • min_cost_flow(S, T, v=1Nmax(0,Dv)) が答え.

これはなぜかというとこまっちゃうんですが, 水が溢れている頂点から水が不足している頂点に最小費用流を流しているイメージです.

実装例を以下に示しました.

template< typename flow_t, typename cost_t >
struct edge {
  int src, to;
  flow_t cap;
  cost_t cost;

  edge() = default;

  edge(int src, int to, flow_t cap, cost_t cost) : src(src), to(to), cap(cap), cost(cost) {}
};

template< typename flow_t, typename cost_t, template< typename, typename > class MF >
cost_t normalized_min_cost_flow(const vector< flow_t > &D, const vector< edge< flow_t, cost_t > > &E) {
  const int N = (int) D.size(), M = (int) E.size();
  MF< flow_t, cost_t > flow(N + 2);
  const int S = N, T = N + 1;
  
  for(auto &e : E) {
    flow.add_edge(e.src, e.to, e.cap, e.cost);
  }
  
  flow_t in = 0;
  for(int i = 0; i < N; i++) {
    if(D[i] > 0) {
      flow.add_edge(S, i, D[i], flow_t(0));
      in += D[i];
    } else if(D[i] < 0) {
      flow.add_edge(i, T, -D[i], flow_t(0));
    }
  }
  return flow.min_cost_flow(S, T, in);
}

S から T へ流量 F のフローを流す場合も, この関数を使って解くことができて, DS=F,DT=F で, 他の Di0 とすればよいです.

以下のコードは S=0,T=N1 の場合について, AOJ(http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=GRL_6_B&lang=jp)でverifyしています.

int main() {
  using flow_t = int;
  using cost_t = int;

  int N, M, F;
  cin >> N >> M >> F;
  vector< edge< flow_t, cost_t > > E(M);
  for(int i = 0; i < M; i++) cin >> E[i].src >> E[i].to >> E[i].cap >> E[i].cost;
  vector< flow_t > D(N);
  D[0] += F;
  D[N - 1] -= F;
  cout << normalized_min_cost_flow< flow_t, cost_t, PrimalDual >(D, E) << endl;
}

C. 最小流量制約付きの最小費用流について

最小流量制約付きの最小費用流について考えます.

  • add_edge(u, v, low, high, cost):= 頂点 u から頂点 v に流量 [low,high], コスト cost の辺を張る

普通(?) の最小費用流と何が違うかというと, それぞれの辺に最小流量の制約がついています. この場合は, ちょっと頑張ると最小流量が 0 の辺に帰着できます.

具体的には, 最初から流量 low だけ水が流れていたことにします(このときのコスト low×cost を最小費用流を求めた後に加算することとします). 流れていたことにするのはさっきの B. 一般化した最小費用流 を用いるとかなり簡単で, Du から low をひいて Dvlow を加算すればよいです.

このあと, 残りの highlow 分の流量についての辺を張れば良くて, add_edge(u, v, highlow, cost) をすれば終わりです.

D. 負辺のある最小費用流について

辺のコストに負のものが存在する場合の最小費用流について考えます. add_edge(u, v, low, high, cost) について, cost<0 の場合です.

最小流量の場合と同じノリでOKです.

この場合はできるだけ流すとオトクなので, 最小から流量 high だけ水が流れていたことにしましょう(C. と同様に high×cost を最後に加算することとします). Du から high をひいて Dvhigh を加算します.

流すのを highlow 分だけキャンセルできます. これは流量 highlow 分について逆辺を張れば良くて, add_edge(v, u, highlow, cost) をすれば終わりです.

ここまでの機能(最小流量, 負辺) を持つ最小費用流を解くための実装例を次に示します. NG はフローを流せなかった場合の戻り値です.

template< typename flow_t, typename cost_t >
struct edge {
  int src, to;
  flow_t low, high;
  cost_t cost;

  edge() = default;

  edge(int src, int to, flow_t high, cost_t cost) : src(src), to(to), low(0), high(high), cost(cost) {}

  edge(int src, int to, flow_t low, flow_t high, cost_t cost) : src(src), to(to), low(low), high(high), cost(cost) {}
};

template< typename flow_t, typename cost_t, template< typename, typename > class MF >
cost_t normalized_min_cost_flow(vector< flow_t > D, const vector< edge< flow_t, cost_t > > &E, cost_t NG = -1) {
  const int N = (int) D.size(), M = (int) E.size();
  MF< flow_t, cost_t > flow(N + 2);
  const int S = N, T = N + 1;

  cost_t sum = 0;
  for(auto &e : E) {
    if(e.cost < 0) {
      sum += e.cost * e.high;
      D[e.src] -= e.high;
      D[e.to] += e.high;
      flow.add_edge(e.to, e.src, e.high - e.low, -e.cost);
    } else {
      sum += e.cost * e.low;
      D[e.src] -= e.low;
      D[e.to] += e.low;
      flow.add_edge(e.src, e.to, e.high - e.low, e.cost);
    }
  }

  flow_t in = 0;
  for(int i = 0; i < N; i++) {
    if(D[i] > 0) {
      flow.add_edge(S, i, D[i], flow_t(0));
      in += D[i];
    } else if(D[i] < 0) {
      flow.add_edge(i, T, -D[i], flow_t(0));
    }
  }

  auto ret = flow.min_cost_flow(S, T, in);
  if(ret == -1) return NG;
  return ret + sum;
}

E. 双対について

双対って知っていますか

僕も知りません

整数性を仮定する整数性とかあるみたいですが, らてまるたがそのうち記事を書きます

F. 最小費用流の双対について

なんか問題の解き方としては真面目に双対をとってがんばる汎用的なやつと, なんかの問題の双対の形を覚えておいてそれに当てはめる方法があるみたいです. ここでは後者のうち最小費用流の双対について定式化して解くことにします.

う し た ぷ に き あ く ん 笑(最小費用流の双対 | beet’s soil) にも書かれている通り, 最小費用流の双対をとると, 次のような問題がでてきます.

max{eloweaehighebe)v=0N1Dvpv}
s.t.
(aebe)+(pvpu)costuv
ae0,be0,pv0

つまり条件が 2 変数の差に関するものだったら, 最小費用流です. 本当か?

各制約について, 頂点 u から頂点 v に容量 [low,high], コスト costuv の辺を張って, 各頂点に対する水の過不足の情報 Dv を与えた最小費用流を求めれば, 目的のものを得ることができます.

以降, 色付き文字(黒ではない)が変数です. これらの値を最小費用流によりうまく設定して目的関数を最大化します.

よくわからないのでここからは問題を解きます. ところで, 双対性 が詳しいので, こんなたぴちゃんみたいな記事よりもこっちを読むと幸せになれるかもしれません.

G. How to Create a Good Game

問題概要は次の通りです.

N 頂点のDAGが与えられて, 各 i について 0iN1 のパスが存在する. i 番目の辺は頂点 xi から頂点 yi に重み si で張られている. DAG の最大長を維持しながら辺に重みを足すとき, 足せる重みの和の最大値を求めよ.

与えられた DAG の最大長を D とします. また ai:=i に足す重み とします. このとき答えは,

maxi=1Mai
s.t.
DAG の最大長が D
ai0

となります. 最大長が D の言い換えが必要です. ここで dv:= 頂点 0 から v までの最大長とします. すると

dxi+si+aidyi(1iM) (a)
dN1d0D (b)

と書くことができます. 最小費用流の双対っぽい形をしていることに気づくと(そういう問題を選んだのでそれはそうなんですが), その形にあわせて書き直すと, 次の通りになります.

max{e=i,j(loweaehighebe)v=1NDvdv}
s.t.
Dv=0
(aibi)+(dxidyi)si(1iM) (a)
lowi=1,highi= (a)
(ajbj)+(dN1d0)D (b)
lowj=0,highj= (b)
ae0,be0,dv0

移項して, lowi,highi,Dv などを適切な係数に設定すると, 上のような形になります(be は常に 0 であってほしいので, highe の値を に設定します).

うくにきあちゃんなのでこの通りに辺を張って, フローを流して終了です.

int main() {
  using flow_t = int64;
  using cost_t = int64;

  int N, M;
  cin >> N >> M;
  vector< vector< int > > g(N);
  vector< int > X(M), Y(M), S(M);
  for(int i = 0; i < M; i++) {
    cin >> X[i] >> Y[i] >> S[i];
    g[X[i]].emplace_back(i);
  }
  auto dp = make_v< int >(N);
  for(int i = 0; i < N; i++) {
    for(auto &j : g[i]) chmax(dp[Y[j]], dp[i] + S[j]);
  }

  vector< edge< flow_t, cost_t > > es;
  vector< flow_t > D(N);
  for(int i = 0; i < M; i++) es.emplace_back(Y[i], X[i], 1, inf, -S[i]);
  es.emplace_back(0, N - 1, 0, inf, dp[N - 1]);
  cout << normalized_min_cost_flow< flow_t, cost_t, PrimalDual >(D, es) << endl;
}

H. 123パズル

beet-aizu.hatenablog.com

I. ラグランジュ双対

ラグランジュ双対ってなんですか 知りません

LPでは, maxf(x)s.t.g(x)0minλ0maxf(x)λg(x) が等しくなります(強双対性).

J. Longest Shortest Path

みんなだいすき Longest Shortest Path. ラグランジュ双対により解ける問題です.

N 頂点 M 辺の重み付き有向グラフが与えられる. i 番目の辺は頂点 ui から vi に重み di で結んでいて, コスト x×ce 支払うことで, 重みに x 加算できる. 予算 P の範囲内で, st 最短経路の長さを最長化せよ.

bi:=i に加算したコストと定義します. やりたいことは

max{st 最短距離 }
s.t
i=1McibiP
bi0

です. ラグランジュ双対をとると,

minλ0{max{st 最短距離 λ(i=1McibiP)}}
s.t.
bi0

となります. ところで, st 最短距離 は, そのまま LP双対で書き直せて

minλ0{max{ptpsλ(i=1McibiP)}}
s.t.
pvipuidi+bi
pv0,bi0

となります. で, もうちょっと目的関数の最小値の内側(最大値部分)を整理すると,

max{i=1Mλcibi+(ptps)}+λP

となります. 例のごとく最小費用流の双対っぽい形をしているので, ねねちゃんをすると,

max{e=i(loweaehighebe)v=0N1Dvpv}+λP
s.t.
Ds=1,Dt=1,Di=0(is,it)
lowi=0,highi=λci
(aibi)+(puipvi)di
ai0,bi0,pv0

となり終了. 最小値部分は三分探索をするとよいです. λ が小さすぎるとフローが存在しない場合があるので, ちょっと工夫してあります(コード参照). ところで, 三分探索をしなくても解けるらしいです.

int main() {
  using flow_t = double;
  using cost_t = double;

  int N, M, P, S, T;
  cin >> N >> M >> P >> S >> T;
  --S, --T;
  vector< int > V(M), U(M), D(M), C(M);
  for(int i = 0; i < M; i++) {
    cin >> U[i] >> V[i] >> D[i] >> C[i];
    --U[i], --V[i];
  }

  auto check = [&](double lambda) -> pair< bool, double > {
    vector< flow_t > X(N);
    X[S] = -1.0;
    X[T] = 1.0;
    vector< edge< flow_t, cost_t > > E;
    for(int i = 0; i < M; i++) E.emplace_back(V[i], U[i], 0.0, C[i] * lambda, D[i]);
    auto ret = normalized_min_cost_flow< flow_t, cost_t, PrimalDual >(X, E, inf);
    if(ret == inf) return {false, -1};
    else return {true, ret + lambda * P};
  };

  double left = 0.0, right = 1.0;
  for(int i = 0; i < 60; i++) {
    double a = (left * 2 + right) / 3;
    double b = (left + right * 2) / 3;
    auto L = check(a);
    auto R = check(b);
    if(!R.first) left = b;
    else if(!L.first) left = a;
    else if(L.second < R.second) right = b;
    else left = a;
  }
  cout << check(right).second << endl;
}

K. 飽きた

飽きました サヨナラー

L. う し た ぷ に き あ 王 国 笑 のすすめ

今, う し た ぷ に き あ 王 国 笑 に入ると, じょえくんからささやかなプレゼントが贈られます しらんけど

I. びーと

よるごはんをたべにいこうぜ

J. らて

またあそぼうね

K. たぷ

たぷちゃんって可愛くないですか 一番好きな天使です!

L. う