Hatena::Grouptopcoder

(iwi) { 反省します

TopCoder: [[iwi]] / Twitter: @iwiwi

|

2013-08-21

「ZDD+フロンティア法」はコンテストで役に立つか

13:44 | 「ZDD+フロンティア法」はコンテストで役に立つか - (iwi) { 反省します を含むブックマーク はてなブックマーク - 「ZDD+フロンティア法」はコンテストで役に立つか - (iwi) { 反省します 「ZDD+フロンティア法」はコンテストで役に立つか - (iwi) { 反省します のブックマークコメント

ZDD, フロンティア法とは

問題候補

関連しそうなものを思いつく限り列挙したが,残念ながら,これらの問題では ZDD を使う利点はあまりないように思う.その理由は以下.

(他に何かご存知でしたら教えてもらえると喜びます)

考察

ZDD を使うための必須条件は,まず

  • 扱うものが連結成分,木,森,パス,サイクルあたり(グラフ以外の問題への適用まで考えるとこの限りではないが)
  • メモリ使用量にかなり余裕がある

ということ.例えば,DP なら O(WH 3^W) 時間でも配列を使い回せば O(3^W) スペースでできる問題でも,ZDD であれば必ず O(WH 3^W) スペース(しかもそこそこの定数倍)が必要.

さらに,ZDD を使って有利になるための条件は,例えば

  • 1 つ見つけたいのではなく,最適なものを見つけたかったり全て数え上げたかったり,いわゆる DP 的処理(or列挙)がしたい
  • DP 中の状態遷移の中で書くと大変だが,構成後の ZDD に後から Family Algebra で変更を加える方針なら凄い簡単になるような制約がついている

というようなものが考えられる.

まず,DP は,複雑であればあるほど ZDD が有利になる.これは,普通の DP 解法では状態遷移をしながらそれに加えてさらに最適値の更新も行わなければならないのに対し,ZDD の構築部分では状態遷移だけを考えればよいから.しかし,実際にはこういう系の問題はそもそも状態遷移だけで凄い複雑なので,変な目的関数は滅多にでない.また,扱いやすいライブラリ化に繋がるかもしれないが,次の段落で述べるように,実際には毎回のように状態遷移に手を入れることになりそう.(あと,集合に含まれるか否かの情報しかない ZDD にしてしまってからでは計算できないような DP もある.)

次に,後で Family Algebra で解への制約を追加できるのは,実装面で楽になるかもしれない.ただ,上の問題たちを見ていると,実際にはこれはあまり上手く活用できなさそう.というのも,1 つの制約の追加に,ZDD 全体のサイズが計算量でかかってしまうから.例えば,「隣り合うマスを使わない」みたいなルールの場合,制約を ZDD 構成後に入れようとすると,入れる制約の数が O(WH) 有り,そんなに制約を後から入れると計算量がやばい(O(W^2H^2 3^W) とかになる).制約が定数個に収まっていないとこの方法は使えず,実際にはあまりこの方法を取れる機会は多くないかも.

2013-08-16

ICFP Programming Contest 2013 感想など

02:23 | ICFP Programming Contest 2013 感想など - (iwi) { 反省します を含むブックマーク はてなブックマーク - ICFP Programming Contest 2013 感想など - (iwi) { 反省します ICFP Programming Contest 2013 感想など - (iwi) { 反省します のブックマークコメント

2013/09/27 追記:

優勝しました!

こっちの記事を見て下さい → http://d.hatena.ne.jp/iwiwi/20130927/1380255924



チームについて

参加記を書こうかとも思ったけれど,だるいので,適当に書きたい事だけ書きます.

今年のチーム

順番は,僕視点で仲間になっていった順.チーム名は "Unagi: The Synthesis" です(何気にたぶん未公開情報).


チーム発足まで

このチームは,2 年前に準優勝(http://d.hatena.ne.jp/iwiwi/20110920/1316528846)したチーム "Unagi: The Gathering" をベースとしています.優勝者は 301 点,準優勝の我々は 300 点.とても悔しい思いをしました.

去年は,僕が VKCup(http://topcoder.g.hatena.ne.jp/iwiwi/20120715/1342360532)でロシアに行くのとかぶっていたこともあり崩壊.しかし主砲である wata は一人で出て Lightning 部門で優勝していた,流石すぎる(https://twitter.com/wata_orz/status/245500307205935104).

そして今年.去年がつまらなそうな問題だったので,ウーン,と思いつつ,今年はなんとなく Web ページ等から気合入ってそうなオーラを感じ取り,またやるか!という気分に.また本気でやりたいなあという事で wata と話しているところに imos からも声がかかったので,とりあえず 2 年前のチームメイトに声をかける.でも ir5, hos が日程の都合がつかないということで,一旦 4 人に.

2 年前のチームはあれはあれで1つの究極のチームだったと思います.例えば,TopCoder の平均レーティングがとにかくやばい.でも,過去の反省からするに,もっと幅広い力を持つほうが有利だということで,chokudai と sulume に声をかけてチーム結成.これは,予想以上に(特に,今年の問題には)完璧なチョイスでした.


問題について

このへんは,本当に人によって意見が異なる点だと思います.あくまで個人的な感想です.

今年の問題に関する自分の個人的な感想は,こんなところかなあと思います.

  • 重大な問題点はなし(Unable to decide equality は若干イラッ)
  • 運営は完璧(サーバ落ちない.問題文に欠陥がなく読みやすい.向こうのプログラムにバグがない.メールの反応もよい.)
  • そこそこ面白い,でもちょっとストイックすぎて寂しい

どうしても対比したくなるのは,3 年前の,車の出荷の問題ですね.こっちで計算をおこなって,サーバにある大量の問題を解いていくスタイルが類似.あれと対比すると,僕個人は,以下のように感じました.

  • あの時と比べて運営が段違いに良くてストレスがない
  • でもそれ抜きで考えれば,実は問題はあの時のほうが面白かったかも

問題の面白さなんかでいうと,やっぱり 2 年前も神がかり的に楽しかったと感じていて,今年はそれと比べてもやっぱりちょっと劣るかなと思います.

なんでそう感じるかというと,絡んでくる要素が少ない点じゃないかと思います.2 年前とか 3 年前の問題は,いろいろな面でやることがあって面白かった.そもそも,自分が 2 年前から ICFPC に本気参加しているのも,3 年前の優勝チームの人たち(http://d.hatena.ne.jp/tanakh/20100702)が色んな要素で強さを見せつけていてチョーカッケー!!俺もこんな総合力勝負っぽいコンテストでも勝ってみてー!!!って思ったから.なので,その点,今年はなんやかんや殆どが探索のチューニングだったなぁと思う(あと実行システム).まぁ,それも十分楽しかったんだけど.


ICFPC で勝つためには

色々事前に(チームを組む前から)考えていて,今年やってみて改めて色々感じて,色々あるけど,だるいからまた今度.

ウチのチームは,あとは皆のプログラミング言語が共通してないというのだけがとにかく重大な問題だと思うのだよなあ.2 年前みたいな,大きめの1つのプログラムを作らないといけないとき,一気に行動が制限される.この辺,皆 C++Python が書ける Googler で結成してるチームはかなり有利だなあと思う(まぁうちの Googler は PHP を書いてたんだけど).


結果について

wata (マジキチ) さんのツイートをご覧ください.

来年のためのメモ

2013-01-14

スキー

00:27 | スキー - (iwi) { 反省します を含むブックマーク はてなブックマーク - スキー - (iwi) { 反省します スキー - (iwi) { 反省します のブックマークコメント

皆様!!!!!スキー!!!!行きましょう!!!!

http://www.zusaar.com/event/489059


f:id:iwiwi:20130115002810g:image

2013-01-05

pow_mod, inverse の速度

14:22 | pow_mod, inverse の速度 - (iwi) { 反省します を含むブックマーク はてなブックマーク - pow_mod, inverse の速度 - (iwi) { 反省します pow_mod, inverse の速度 - (iwi) { 反省します のブックマークコメント

久々にコンテストの問題が解きたくなって AOJ で 問題 MinimumCostPath を解いた.個人的にはムチャクチャ面白いと思った.

ところで,別に TLE はしなかったんだが,mod の inverse がボトルネックだった.自分のライブラリの mod の inverse とからへん,コンテスト始めたころの大昔に作ってから変えてなくて,今やあまり気に入っていないと思いつつバグってはいないということで長らく放置していたが,この際ようやく弄ることにした.


ベンチマークコード

適当.g++ -O2.標準入力は適当に10とか入れてる.

const ll MOD = 1000000009;

const int N = 2000000;

int main() {
  int a;
  scanf("%d", &a);

  ll t = a;
  for (int k = 1; k <= N; ++k) {
    t = inverse(t, MOD) * k % MOD;
  }
  printf("%lld\n", t);

  return 0;
}

追記:手元の環境は

~/contests/tmp% cat /proc/cpuinfo [00:35:39]

processor : 0

vendor_id : GenuineIntel

cpu family : 6

model : 26

model name : Intel(R) Core(TM) i7 CPU 960 @ 3.20GHz

stepping : 5

cpu MHz : 1600.000

cache size : 8192 KB

~/contests/tmp% gcc -v [00:36:44]

Using built-in specs.

Target: x86_64-linux-gnu

Configured with: ../src/configure -v --with-pkgversion='Ubuntu 4.4.3-4ubuntu5.1' --with-bugurl=file:///usr/share/doc/gcc-4.4/README.Bugs --enable-languages=c,c++,fortran,objc,obj-c++ --prefix=/usr --enable-shared --enable-multiarch --enable-linker-build-id --with-system-zlib --libexecdir=/usr/lib --without-included-gettext --enable-threads=posix --with-gxx-include-dir=/usr/include/c++/4.4 --program-suffix=-4.4 --enable-nls --enable-clocale=gnu --enable-libstdcxx-debug --enable-plugin --enable-objc-gc --disable-werror --with-arch-32=i486 --with-tune=generic --enable-checking=release --build=x86_64-linux-gnu --host=x86_64-linux-gnu --target=x86_64-linux-gnu

Thread model: posix

gcc version 4.4.3 (Ubuntu 4.4.3-4ubuntu5.1)

今までのやつ

駆け出しだった頃の僕が,適当に bmerry か誰かのやつを参考にして書いた.中で再帰している上にいちいち丁寧に正にするために % を 2 回もとってて遅そう.(でもマイナスが帰るとそれが原因のバグをつくりそうで恐ろしいという気持ちは理解できる)

inline ll mod(ll a, ll m) { return (a % m + m) % m; }

ll inverse(ll a, ll m) {
  if ((a = mod(a, m)) == 1) return 1;
  return mod((1 - m * inverse(m % a, a)) / a, m);
}

~/contests/tmp% time ./a.out <<< 10 [14:10:35]

249011238

./a.out <<< 10 2.84s user 0.01s system 99% cpu 2.870 total

~/contests/tmp% time ./a.out <<< 10 [14:10:39]

249011238

./a.out <<< 10 2.85s user 0.00s system 99% cpu 2.870 total

~/contests/tmp% time ./a.out <<< 10 [14:10:43]

249011238

./a.out <<< 10 2.83s user 0.02s system 99% cpu 2.868 total


inline

再帰関数だからあんま意味分からんと思うが一応 inverse に inline つけてみた.

~/contests/tmp% g++ -O2 inverse.cpp [14:11:39]

~/contests/tmp% time ./a.out <<< 10 [14:11:41]

249011238

./a.out <<< 10 2.75s user 0.00s system 99% cpu 2.763 total

~/contests/tmp% time ./a.out <<< 10 [14:11:44]

249011238

./a.out <<< 10 2.75s user 0.00s system 99% cpu 2.762 total

~/contests/tmp% time ./a.out <<< 10 [14:11:49]

249011238

./a.out <<< 10 2.75s user 0.00s system 100% cpu 2.748 total

ちょっぴりはやくなったw



pow_mod 1

今まで使っていた pow_mod を使う.これも再帰していて丁寧に mod を呼んでいて見るからに遅そうである.

inline ll mod(ll a, ll m) { return (a % m + m) % m; }

ll pow_mod(ll a, ll k, ll m) {
  if (k == 0) return 1;
  ll t = pow_mod(a, k / 2, m);
  return mod(mod(t * t, m) * (k % 2 ? a : 1), m);
}

inline ll inverse(ll a, ll m) {
  return pow_mod(a, m - 2, m);
}

~/contests/tmp% time ./a.out <<< 10 [14:13:19]

249011238

./a.out <<< 10 3.39s user 0.00s system 99% cpu 3.407 total

~/contests/tmp% time ./a.out <<< 10 [14:13:23]

249011238

./a.out <<< 10 3.39s user 0.00s system 99% cpu 3.407 total

~/contests/tmp% time ./a.out <<< 10 [14:13:28]

249011238

./a.out <<< 10 3.37s user 0.02s system 99% cpu 3.411 total

遅くなったw



pow_mod 2

pow_mod を高速そうな実装にしてみる.

inline ll pow_mod(ll a, ll k, ll m) {
  ll r = 1;
  for (; k > 0; k >>= 1) {
    if (k & 1) (r *= a) %= m;
    (a *= a) %= m;
  }
  return r;
}

~/contests/tmp% time ./a.out <<< 10 [14:16:27]

249011238

./a.out <<< 10 0.45s user 0.01s system 99% cpu 0.463 total

~/contests/tmp% time ./a.out <<< 10 [14:16:28]

249011238

./a.out <<< 10 0.47s user 0.00s system 100% cpu 0.468 total

~/contests/tmp% time ./a.out <<< 10 [14:16:31]

249011238

./a.out <<< 10 0.45s user 0.01s system 99% cpu 0.463 total

かなりはやくなったw

少し予想はしていたが,今までのはやはり爆遅だった



非再帰 exgcd

Mi_Sawa さんに教えてもらった(ありがとうございます).確かに速そう.

inline ll exgcd(ll a, ll b, ll &x, ll &y){
  ll u[] = {a, 1, 0}, v[] = {b, 0, 1};
  while(*v){
    ll t = *u / *v;
    rep(i, 3) swap(u[i] -= t * v[i], v[i]);
  }
  x = u[1]; y = u[2];
  return u[0];
}
inline ll inverse(ll a, ll m){
  ll p, q;
  exgcd(a, m, p, q);
  return (p%m+m)%m;
}
手元環境

~/contests/tmp% g++ -O2 inverse.cpp [23:26:03]

~/contests/tmp% time ./a.out <<< 10 [23:26:08]

249011238

./a.out <<< 10 0.65s user 0.00s system 100% cpu 0.649 total

~/contests/tmp% time ./a.out <<< 10 [23:26:10]

249011238

./a.out <<< 10 0.64s user 0.01s system 100% cpu 0.644 total

~/contests/tmp% time ./a.out <<< 10 [23:26:16]

249011238

./a.out <<< 10 0.64s user 0.01s system 100% cpu 0.644 total

自分の環境だと僅かに pow_mod より遅い.

ideone

しかし,ideone に貼りつけたら,exgcd のほうが 2 倍ぐらい速い.

Codeforces

ideone と似たっけかになった.

pow_mod 0: 1.015625 sec

pow_mod 1: 1.014648 sec

pow_mod 2: 1.014649 sec

exgcd 0: 0.695312 sec

exgcd 1: 0.694336 sec

exgcd 2: 0.694336 sec

TopCoder

やっぱしexgcd のほうが 2 倍ぐらい速い.

pow_mod 0: 0.369712 sec

pow_mod 1: 0.351320 sec

pow_mod 2: 0.351127 sec

exgcd 0: 0.171135 sec

exgcd 1: 0.170425 sec

exgcd 2: 0.170366 sec

結論

exgcd のほうが速い環境がどうやら多い,pow_mod のほうが速い時も差は僅か

exgcd のほうがよさそうかも



非再帰 exgcd ループアンローリング

ループアンローリングしたらちょっと速くなった

http://ideone.com/Dz1LLI



整理して

inline ll inverse(ll a, ll m){
  ll b = m, u = 1, v = 0;
  while (b) {
    ll t = a / b;
    swap(a -= t * b, b);
    swap(u -= t * v, v);
  }
  return (u % m + m) % m;
}

Mi_SawaMi_Sawa2013/01/05 20:38自分の環境ではexgcdの方が速かったです.
(pow_modだとコンスタントにlog Mくらいかかりますが, exgcdだと速い時がそれなりにあるので, 良さそうでもある)

inline ll exgcd(ll a, ll b, ll &x, ll &y){
ll u[] = {a, 1, 0}, v[] = {b, 0, 1};
while(*v){
ll t = *u / *v;
rep(i, 3) swap(u[i] -= t * v[i], v[i]);
}
x = u[1]; y = u[2];
return u[0];
}
inline ll invMod(ll a, ll m){
ll p, q;
exgcd(a, m, p, q);
return (p%m+m)%m;
}

iwiwiiwiwi2013/01/06 00:23コメントありがとうございます!!
試してみました,詳しくは記事に追記します.自分の環境だと僅かに pow_mod のほうが高速でしたが,ideone で実行したら確かに貼ってくださった exgcd のほうが 2 倍とか速かったです.

Mi_SawaMi_Sawa2013/01/06 14:36もう一つよいのは, Mが素数である必要がなく, 逆元が無い時も,
ax = gcd(a,M)
なるxを見つけてくれる所ですね.
負の数にちゃんと対応していなかったり(普通に正規化していれば大丈夫ですが, 一度これで引っかかった), 逆元の存在のチェックをしていないので, ロバストにするなら,
ll u[] = {llabs(a), 1, 0}, v[] = {llabs(b), 0, 1};
にしたり,
assert(exgcd(a, m, p, q) == 1)
にしておくとよいです.

Mi_SawaMi_Sawa2013/01/06 16:22あ, ll u[] = {llabs(a), 1, 0}, v[] = {llabs(b), 0, 1};だとダメな気がします, コメント汚して申し訳ないです…
ループ終わったあと, u[0]が負ならrep(i, 3) u[i] = -u[i];とかすればよい気がします…

iwiwiiwiwi2013/01/07 11:35ありがとうございます.張り付けで使うならちょっと長くなっても周りのちょっとしたミスも吸収できるようにできるだけロバストにしておくのが良さそうですよね.

2012-12-19

SMT Solver (z3) で探索問題を解く

23:01 | SMT Solver (z3) で探索問題を解く - (iwi) { 反省します を含むブックマーク はてなブックマーク - SMT Solver (z3) で探索問題を解く - (iwi) { 反省します SMT Solver (z3) で探索問題を解く - (iwi) { 反省します のブックマークコメント

ちょっと触る機会があったのでメモ程度に書いておく.Python は普段使わないのでコードはあまりよくないかもしれない.


はじめに

SMT Solver?

SMT Solver とは SAT Solver より色々できるひとたち. SMT って,SAT に比べてあまり名前を聞かないという人も多いかと思うが,SAT Solver の界隈では,実は SAT と SMT は並列にタイトルに入るぐらいには話題沸騰なのである.


z3?

最強の SMT Solver の 1 つっぽい. Microsoft Research で開発されている. Python から簡単に呼べる z3py というのがあって,ウェブからすぐ使ってみることのできるデモまである.すごい.これを見れば大体 SMT Solver でどんなことができるかの雰囲気がよく分かると思う.

最新版を git から clone してビルドしてインストールしたらスグに z3 が python から使えた.ビルドの際,dos2unix がないと駄々をこねていたので,ln -s todos unix2dos してやった.



書いてみたコード

Captain Q's Treasure (ICPC福岡2011 G)

公式データセットに対して 4 秒で走る.とはいえそもそも IP だから SMT Solver で解く意義は薄い.とりあえず定式化が簡単そうなのでやってみた.

from z3 import *

while 1:
    #
    # Input
    #
    h, w = map(int, raw_input().split())
    if h == 0 and w == 0:
        break
    field = [raw_input() for i in range(h)]

    #
    # Constraints
    #
    var_x = [ [ Int("x_%s_%s" % (i + 1, j + 1)) for j in range(w) ] for i in range(h) ]

    solver = Solver()
    solver.add([var_x[y][x] == 0 if field[y][x] == '.' else And(0 <= var_x[y][x], var_x[y][x] <= 1)
                for x in range(w) for y in range(h) ])

    for y in range(h):
        for x in range(w):
            if '0' <= field[y][x] and field[y][x] <= '9':
                xs = []
                for dy in [-1, 0, 1]:
                    for dx in [-1, 0, 1]:
                        cx = x + dx
                        cy = y + dy
                        if cx >= 0 and cx < w and cy >= 0 and cy < h:
                            xs += [var_x[cy][cx]]
                solver.add(Sum(xs) == int(field[y][x]))

    xs = [var_x[y][x] for y in range(h) for x in range(w)]

    #
    # Solve
    #
    low = -1
    high = h * w
    while high - low > 1:
        mid = (low + high) / 2
        solver.push()
        solver.add(Sum(xs) <= mid)
        if solver.check() == sat:
            high = mid
        else:
            solver.pop()
            low = mid
    print high

15 パズル

次に挑戦するものとしてはタフすぎた気配がするw 定式化大変だった.

で,性能だが,3x3 の 8 パズルだとわりと何とかなるぽいが,4x4 の 15 パズルだと 10 手以下のものは解けるけど増えるにしたがってやばくなっていく.

しょんぼりやな.

from z3 import *

god = 20  # limit...

def solve(field, goal):
    #goal = ((1,2,3),(4,5,6),(7,8,0))
    h = len(field)
    w = len(field[0])

    print [h, w]

    x0, y0 = -1, -1
    for x in range(w):
        for y in range(h):
            if field[y][x] == 0:
                x0 = x
                y0 = y

    print [x0, y0]

    var_f = [ [ [ Int("f_%s_%s_%s" % (g, y + 1, x + 1)) for x in range(w) ] for y in range(h) ] for g in range(god + 1) ]
    var_d = [ Int("d_%s" % (g)) for g in range(god) ]
    var_x = [ Int("x_%s" % (g)) for g in range(god + 1) ]
    var_y = [ Int("y_%s" % (g)) for g in range(god + 1) ]

    #
    # Placement
    #
    constraint_init = []

    for g in range(god):
        constraint_init += [var_f[0][y][x] == field[y][x] for x in range(w) for y in range(h) ]
        constraint_init += [var_x[0] == x0]
        constraint_init += [var_y[0] == y0]
        constraint_init += [var_f[god - 1][y][x] == goal[y][x] for x in range(w) for y in range(h) ]

    #
    # Movement
    #
    constraint_move = []
    dx = [1, 0, -1, 0, 0]
    dy = [0, 1, 0, -1, 0]

    for g in range(god):
        cx = var_x[g]
        cy = var_y[g]
        d = var_d[g]
        nx = var_x[g + 1]
        ny = var_y[g + 1]
        constraint_move += [Or([And(d == i, nx == cx + dx[i], ny == cy + dy[i]) for i in range(5)])]
        constraint_move += [0 <= nx, nx < w, 0 <= ny, ny < h]

    #
    # Board
    #
    constraint_board = []

    def gen_move(cf, nf, x, y, d, r, i = 0):
        if i == 4:
            return nf[y][x] == cf[y][x]
        else:
            j = (4 if i == 4 else (i + 2) % 4) if r else i
            tx = min(max(x + dx[j], 0), w - 1)
            ty = min(max(y + dy[j], 0), h - 1)
            return If(d == i, nf[ty][tx] == cf[y][x], gen_move(cf, nf, x, y, d, r, i + 1))

    for g in range(god):
        cx = var_x[g]
        cy = var_y[g]
        cf = var_f[g]
        nx = var_x[g + 1]
        ny = var_y[g + 1]
        nf = var_f[g + 1]
        d = var_d[g]
        for x in range(w):
            for y in range(h):
                constraint_board += [
                    If(And(x == cx, y == cy), gen_move(cf, nf, x, y, d, False),
                       If(And(x == nx, y == ny), gen_move(cf, nf, x, y, d, True),
                          nf[y][x] == cf[y][x]))]

    #
    # Solve!
    #
    solver = Solver()
    solver.add(constraint_init + constraint_move + constraint_board)
    for target in reversed(xrange(god)):
        solver.add(var_d[target] == 4)
        print [target, solver.check()]

solve(((1,2,3,4),(9,5,6,7),(13,10,11,8),(0,14,15,12)),
      ((1,2,3,4),(5,6,7,8),(9,10,11,12),(13,14,15,0)))
|