2015年1月25日日曜日

Android端末をサーバにしてPyramidフレームワークを利用する

このブログ記事をはてなブックマークに追加

しばらく前に、PythonのWebフレームワークであるPyramidを利用した。これがなかなか良くできており、Android端末上でも動かしてみたくなったので載せてみた。

ところで、自分が利用しているキャリアはドコモなんだけど、spモードだとグローバルIPが割り振られないので外部から端末にアクセスできない。なので、spモードを契約せずに、mopera Uを利用している。mopera UであればグローバルIPが割り振られるのでアクセスすることができるからだ。このためだけに、spモードにせず、mopera Uにしていると言っても過言ではない。

閑話休題。まず、Pyramidを動作させるにはAndroid端末用のPython環境であるSL4Aを入れる必要がある。次にPyramidを入れるのだが、必要なモジュールなどが複数あるのでそれも一緒に入れる。一応、Hello Worldプログラムを動かすのに必要なものはすべて挙げたが、もし足りない場合は実行時に何が足りないかエラー表示が出るので、それを参考に入れて欲しい。また、Android SDKが導入されていることを前提にしている。


まず、Android端末にシェルで入ってプログラムを展開するディレクトリを準備する。環境に合わせてディレクトリは読み替えて欲しい。

> adb shell * daemon not running. starting it now on port 5037 * * daemon started successfully * $ cd /sdcard/sl4a/scripts $ mkdir pyramid $ exit

適当なディレクトリで以下のプログラム・モジュールをPC側で展開する。


PC側からAndroid端末にプログラム・モジュールをコピーする。

> adb push pyramid /sdcard/sl4a/scripts/pyramid/pyramid/ > adb push translationstring /sdcard/sl4a/scripts/pyramid/translationstring/ > adb push venusian /sdcard/sl4a/scripts/pyramid/venusian/ > adb push webob /sdcard/sl4a/scripts/pyramid/webob/ > adb push zope /sdcard/sl4a/scripts/pyramid/zope/ > adb push repoze /sdcard/sl4a/scripts/pyramid/repoze/

pkg_resources.pyについては、PCに入っているPythonのsite-packagesなどから適当に持ってきてコピーする。

> adb push pkg_resources.py /sdcard/sl4a/scripts/pyramid/

これでプログラム・モジュールについては準備完了だ。次に、Android端末側の設定をする。

まずはSL4Aを起動させる。pyramidディレクトリが作成されているはずなので、そこをタップして中を確認してみる。一通り、ファイルがあることを確認したら、メニューからViewをタップし、Interpretersを選択する。更にメニューを開き、Start Serverをタップして、Publicを選択する。これで、Android端末がサーバとして機能する。Androidの通知画面を開くと、SL4A Serviceという通知があるので、それをタップすることでServer情報を確認できる。そこに表示されたドメイン名でAndroid端末にアクセスできる。以前の記事でもサーバにする方法を書いたので、そちらも参考にしてほしい。

これで、Android端末をサーバとして機能させることができたので、こんどはHello Worldプログラムを動かして、動作確認をしてみる。

プログラムをhelloworld.pyとして保存した後、Android端末にコピーする。

> adb push helloworld.py /sdcard/sl4a/scripts/pyramid/

そして、コピー先のディレクトリにhelloworld.pyが表示されるので、それをタップして実行する。先ほど得られたドメイン名(ここでは 12345.mopera.net とする)にポート番号8080でアクセスする。


http://12345.mopera.net:8080/hello/world

これで以下の画面が表示されれば成功だ。


携帯の電波が届くところであればどこでも使えて、Pyramidが載っている手のひらWebサーバってなんか良い。

続きを読む...

2015年1月15日木曜日

IRCボットの作成

このブログ記事をはてなブックマークに追加

久々のブログ更新。書きたいときに書くというスタンスです。最近、会社でも技術ブログを立ち上げたので、そちらで書いても良いかも。

今回はIRCボットの作成について書いてみる。しばらく前から社内IRCで技術関連の雑談用チャンネルを利用している。技術系のメーリングリストもあるのだけど、わざわざメールで話題にするようなことでもない、軽い内容を気軽に議論できる場が欲しかったので有志で立ち上げた。IT企業にとって技術的な雑談をできる場はとても大事だと思う。

素のIRCの場合、発言内容がサーバに保存されず、途中から参加しても話題についていけないという問題が出たので、適当なボットで対応することにしたが、やはりエンジニアがメインの会社なのだからボットも自作が基本でしょ、ということでPythonで作ってみた。まあ、作ったといっても「Python でシンプルな IRC クライアントを作成する」のサイトを参考にさせてもらったのでスクラッチから作成したわけではない。

コマンド techlog: を入れると、その日の発言がボットからのtalkで返される。ただ、最初は1日分の発言しか取れない仕様だったので、日を跨ぐとその前の発言が取得できないし、一時的にチャットを抜けた分のログが欲しいだけでも強制的に1日分の発言を取得してしまう。

そこで、SQLite3を使ってユーザ毎にログイン・ログアウト時刻を管理して、ログアウトしてからログインするまでの発言を取得できるように変更することにした。新規に参加したメンバーでも自動的にデータベースに登録される。せっかくユーザの時刻情報があるのだから、techtime:, techtimeall: というコマンドを作って、それぞれ自分の滞在時間、ユーザ全員分の滞在時間も取得できるようにしてみた。

IRCサーバ上で以下のコマンドを実行すると、IRCにボットが常駐する。

$ ./irc_bot.py > irc_bot.log &

ソースコードは以下の通り。


そこそこ社内IRCも使えるようになったのだが、最近、別プロジェクトで利用を始めたSlackがかなり便利だったので、そこに技術系チャンネルを作ったところ、ほとんどのメンバーがSlackに流れてしまった…。次はSlackでボット作成かなぁ。

続きを読む...

2014年1月25日土曜日

スピログラフ

このブログ記事をはてなブックマークに追加

スピログラフに初めて触れたのは小学校に上がる前だった。当時はそれがスピログラフという名称であることを知らず、輪っかをぐるぐる回して綺麗な模様を描ける道具としてお気に入りの玩具になっていた。自分はマンデルブロ集合やジュリア集合のような数学的に表現できる図形や模様にともて魅力を感じるのだが、幼少のころのスピログラフがその切っ掛けだったのかもしれない。

長方形の定規に大きさの異なる3枚の歯車が付属しており、そして2つの円がくり抜かれている。さらに、矢印、五角形、半円などの型が子供心をくすぐったものだ。因みに、歯車と円に付いている歯数は刻印されており、3つの歯車はそれぞれ36、52、63、円の方は96、105の歯を持っている。

幼少のころ以来、スピログラフを見かけることはなかったのだが、最近になって100円ショップでも売られているという噂を聞き、少し探ってみた。近くの100円ショップでは扱っていなかったが、「デザイン定規」という名称でネットで78円で売っていたので思わず購入してしまった(当然だが送料のほうが高くついた…)。


スピログラフは一般的には内トロコイド(hypotrochoid; 内余擺線)といい、以下の式で表現できる(Wikipediaより)。


定円の半径 rc、動円の半径 rm、描画点の半径 rd とし、下図では、それぞれ、青、緑、赤の線で示している。そして、回転角を θとした軌跡がスピログラフの模様となる。


スピログラフでは半径よりも、円周の歯の数で示した方が分かりやすいだろう。円周と半径の比は変わらないのだから、歯の数の比がそのまま半径の比となる。例えば定円の歯が30で動円の歯が10ならば、rc:rm=3:1ということだ。

さて、幼少のころはこんなに綺麗な模様が数式で表現できるなどと夢にも思わなかったわけだが、今はそれを理解し、コンピュータを使って自由に描くことができる。コンピュータであれば、動円が定円の外側に来る外トコロイド(epitrochoid; 外余擺線)も描くことができるし、描画点を動円の外に出すことも思いのままだ。

せっかくコンピュータを使って簡単に表示できるようになったというのに、プログラムを組むことが難しくては意味がない。しかし、Python+matplotlibであれば簡単に描くことができる。上述の式をそのままコードに落とすだけだ。以下がIPythonで描いたスピログラフの模様だ。図と合わせて一画面に収まってしまう。


ここでは、Pythonのソースコードも挙げておこう。せっかくなので少しだけmatplotlibの調整をしてみた。必ず縦と横のスケールが同じになるようにし、グラフ表示も正方形にした。コードを見てもらえばやり方はすぐに解るだろう。以下のコマンドで実行できる。

$ spirograph.py rc rm rd

例えば次のように実行すれば冒頭の図が表示される。

$ spirograph.py 96 63 30


最後に参考にしたサイトを挙げておく。


続きを読む...

2014年1月1日水曜日

2013年と2014年の違い

このブログ記事をはてなブックマークに追加

2013年を終え、2014年を迎えた。さて、昨年と今年の違いは何か?

昨年のブログで2013を素因数分解すると「3×11×61」になることを示したが、2014は「2×19×53」となる。そして差を取ると、「2-3, 19-11, 53-61」、つまり、「-1, 8, -8」となり、1月、8月が昨年と今年の差であることがわかる。そして、負を示した1月は運気が下がり、正と負を持つ8月は運気の上下があるに違いない! 因みに2013+2014=4027が素数にあることから、昨年から今年にかけて行っていることは、これまでとは違ったユニークな内容となる!

…などと、MMRばりなことを書けるぐらいの余裕を持ちつつ、今年一年を頑張りたいと思いますので、本年もどうぞよろしくお願いいたします。

以下にPython+NZMATHを使った素因数分解を示す。

>>> import nzmath.factor.methods as methods >>> methods.factor(2013) [(3, 1), (11, 1), (61, 1)] >>> methods.factor(2014) [(2, 1), (19, 1), (53, 1)] >>> methods.factor(2013+2014) [(4027, 1)]

続きを読む...

2013年9月21日土曜日

クッキークリッカーとプログラマ

このブログ記事をはてなブックマークに追加

先週ぐらいからクッキークリッカー(Cookie Clicker)というJavaScriptを使ったブラウザゲームが流行っている。クリックするだけのゲームと聞いて、最初はあまり興味を持てなかったのだが、自分の周りであまりにもやっている人が多いので少し遊んでみることにした。

クッキークリッカーを簡単に説明すると、まずはクリックすることでクッキーを作り、作ったクッキーを使ってクッキーの生産性を高めるためのアップグレードやアイテムを購入し、たまに出現するゴールデンクッキー(Golden Cookie)をクリックすることでさらに多量のクッキーが得られるので、それらを駆使してできるだけたくさんのクッキーを作るというゲームだ。

このようにとてもシンプルなゲームなのだが、最初はちまちまとしか作れなかったクッキーが徐々に増えていき、様々なアイテムやイベントを通すことで、終盤では毎秒数千億クッキーを作れるというところに人々を惹きつける魅力があるらしい。しかし、自分は「人間の代わりにプログラムを働かせることでどこまで効率が良くなるか」という目的のために自動化プログラムを作成したくなった。コンピュータは人間を楽にさせるために存在するべきだ。

このゲームはJavaScriptで作られており、そのコードや変数を修正すれば簡単にチートを行うことができる。しかし、自分はそれについては全く興味がない。元のシステムには全く手を入れずに人間が操作する部分のみをプログラムで最適化したいのだ。

まず、最初に考えついたのは自動クリックだ。これは様々なツールが世に溢れているのでそれを使っても良いのだが、プログラマだったら、このぐらいは自分で作りたい。そこで、Pythonとpymouseを使って自動クリックプログラムを作成した。これで、毎秒100~200クリックすることができるようになった。これ以上の速度にするとブラウザがついて行けずプチフリーズを頻繁に起こし逆に生産性が落ちてしまう。因みに画面上のエフェクトなどは設定で消しておくほうが良いだろう。

次に、ゴールデンクッキーへの対処だ。ゴールデンクッキーはランダムな間隔でブラウザ上のどこかに出現する。そのクッキーをクリックすると溜まっているクッキーの1割分の枚数が貰えたり、一定時間現在の7倍の効率でクッキーを焼けるようになったりする。この効果はかなり大きく是非ともこれを自動化したい。

最初は一般的なテンプレートマッチングにより検出しようと思ったが、ゴールデンクッキーは360度ランダムに回転した状態で出現するので、それに対処しようとすると、PILで簡単に済ませるには処理時間が掛かってしまいそうだし、OpenCVを入れるのも面倒だし、クッキー1枚ごときの検出で手間を掛けたくない。そこで、ゴールデンクッキーだけに使われている色を検出することにした。まず、PILでデスクトップのキャプチャを行い、ブラウザの領域のみをクロップし、その領域内でゴールデンクッキーで使われている色を検出できれば、その位置がクリックする位置になる。検出されなければクリックを行わない。

以上の自動化でかなり効率よく稼げるようになり、最終的に作ったクッキーの約7割がクリックによる生産になった。1回目は10京個以上のクッキーを作ったあとにResetを行い、2回目でアップグレードが100%になるまでに要した時間は「リサーチ」の待ち時間を含めて9時間程だった。途中で一時的に中断したり、効率的にアップグレードしたわけではないので、その辺を考えてクッキーを焼けばもっと早く100%になったと思う。また、9時間と言っても自動化プログラムが動いている間は席を離れていても良いので、自分で操作したのはアイテム購入ぐらいだったが。

というわけで、自分はクッキークリッカーをそれなりに楽しんだ。楽しんだ部分は主に自動化プログラムを作成したところなので、一般にプレイされている方とは楽しみ方が違うのかもしれないけど。最後に作成したプログラムを以下に示しておく。


続きを読む...

2013年5月26日日曜日

Pythonを使って簡単にデータを視覚化する

このブログ記事をはてなブックマークに追加

世の中のことをもっと知るにはどうしたら良いだろうと思うときがある。世の中の多くの事柄はログやデータに落とされる。Googleなどの検索サイトは良い例だろう。さて、そのログやデータをどうすれば良いのか? 多くの場合、視覚化が有効な手段となる。

まずは身の回りの日常的なデータやログを何とかしたい。ただ、日常のデータを視覚化するのに数十行以上のコードは書きたくない。まるで息をするかのごとく自然に視覚化を行いたいのだ。そのためには1~2行、長くて数行で済ませることが必要だ。そこでPythonとmatplotlibを使う。加えて、IPythonがあればなお良い。IPythonの導入については以前のブログ記事であるIPythonの埋め込みプロットが素晴らしいを参考にして欲しい。

まずは事前にnumpyとmatplotlibをインポートしておく。できればscipyも。

>>> from numpy import * >>> from pylab import *

短いコードで視覚化を行うためには、Pythonの内包表記は必須だ。例えば、5, 2, 1, 5, 8をデータとするグラフを書きたいのならIPythonを使って以下の1行で実現できる。

>>> plot([5, 2, 1, 5, 8])



数値が行ごとにinput.txtというファイルに書かれていた場合は以下の通り。

>>> plot([int(x) for x in file("input.txt")])

map関数を使えばもっとスマートに書ける。ファイル内の数値を文字列として配列で読み込み、それらの文字列をmap関数によりintで整数に変換する。

>>> plot(map(int, file("input.txt")))

さて、ログが2次元で書かれていた場合はどうするのか。例えば以下のデータがinput.txtに書かれていたとする。

2, 2.5 5, 6.2 6, 3.6 7, 6.3 10, 1.9

この場合、以下のようにすればプロットできる。ファイルから行ごとの文字列を読み出して、それを","を区切りとしてリスト化する。文字列で格納されている数値をmap関数を使ってfloatで実数に変換している。最後にmap(list, zip(*data))を使って転置している。ここでdataはリストによる行列とする。リストによる行列の転置はイディオムとして覚えておくと便利だ。

>>> x, y = map(list, zip(*[map(float, xs.split(",")) for xs in file("input.txt")])) >>> plot(x, y)



ログファイルなどの場合、単に数値が並んでいることは少なく一緒にテキストが書き込まれていることが多い。例えば次のログ(input.log)について視覚化を行なってみる。

Data analysis log: Cycle 1: Calculating ... done! Factor = 0.265 Unit = 25.8 Cycle 2: Calculating ... done! Factor = 0.315 Unit = 28.2 Cycle 3: Calculating ... done! Factor = 0.415 Unit = 30.5 Cycle 4: Calculating ... done! Factor = 0.625 Unit = 40.6 Cycle 5: Calculating ... done! Factor = 0.895 Unit = 55.2

Unitの値をそのままプロットするなら以下のようにすれば良い。予めreモジュールをインポートしておき、re.matchで先頭の文字列と一致したものだけ選び出している。

>>> plot([float(l.split("=")[1]) for l in file("input.log") if re.match("Unit =", l)])



Unitを横軸、Factorを縦軸にしてドットでプロットしたい場合は以下のようにする。

>>> x = [float(l.split("=")[1]) for l in file("input.log") if re.match("Unit =", l)] >>> y = [float(l.split("=")[1]) for l in file("input.log") if re.match("Factor =", l)] >>> plot(x, y, "o")



最後に少し複雑な例を挙げておく。あるデータに対してスペクトラルクラスタリングを行いたいとする。まずは、numpyとmatplotlibの他に、scipy関連の固有値問題・k近傍法・kd木のモジュールをインポートしておく。

>>> from scipy.linalg import eig >>> from scipy.cluster.vq import kmeans2 >>> from scipy.spatial.kdtree import KDTree

そして、クラスタリングしたいデータをcircles.datとすると、以下のコードで視覚化できる。

>>> ps = array([map(float,l.split()) for l in file("circles.dat")]) >>> kt = KDTree(ps) >>> knn = [[((lambda a,b:exp((-(sqrt(dot(a-b,(a-b).conj()))**2))/1.5))(p,ps[nb]),nb) for nb in kt.query(p,16)[1] if i!=nb] for i,p in enumerate(ps)] >>> W = zeros([len(knn)]*2) >>> _ = [W[p].__setitem__(nb,d) for p,nn in enumerate(knn) for d,nb in nn if p in zip(*knn[nb])[1]] >>> w,v = map(real,eig(diag([sum(Wi) for Wi in W])-W)) >>> Y = array([e[1] for e in sorted(zip(w,v.T))[:3]]).T >>> res,idx = kmeans2(Y,3,minit='random') >>> _ = [plot(p[0],p[1],('ro','go','bo')[i]) for p,i in zip(ps,idx)]



ここではforループを行うために内包表記を使っており、また、内包表記のリストが必要ない場合は、_(アンダースコア)に入れておく。これはシェル上で余計なデータを表示させないためである。また、内包表記上では配列への代入ができないので__setitem__を利用している。さらに、Pythonのリストの代わりにnumpyのarrayを利用すれば転置などもmap(list, zip(*data))の代わりにarray(data).Tとすれば良いので分かりやすい。

Pythonを使えばたったこれだけだ。これだけのことで、これまで見えなかったものが見えてくるのは楽しい。

続きを読む...

2013年3月10日日曜日

Python: timeitモジュールを使ってお手軽に実行時間計測

このブログ記事をはてなブックマークに追加

Pythonではtimeitモジュールを使って簡単に実行時間を計測できる。以下に示すコードはPythonでは一般的なコードだと思う。たらい回し関数(竹内関数)と階乗を求める関数だ。これを実行した時の時間計測、さらにコードに含まれている関数の実行時間を簡単に測るにはどうすればよいか。


ここでtimeitモジュールが利用できる。まず、Pythonシェル環境では以下の通り。

main関数の測定は次のように行う。stmtで時間計測したい関数を指定して、setupで最初の1回だけ実行する文を指定する。numberは実行する回数となっている。正確を期するなら複数回を指定するべきだろう。

>>> import timeit >>> timeit.timeit(stmt='test_timeit.main()', setup='import test_timeit', number=1) 12 3628800 2.663659573618423

たらい回し関数(tarai)のみの測定。

>>> timeit.timeit(stmt='test_timeit.tarai(12, 6, 0)', setup='import test_timeit', number=1) 2.6390463109480637

階乗関数(factorial)のみの測定。100万回実行。

>>> timeit.timeit(stmt='test_timeit.factorial(10)', setup='import test_timeit', number=1000000) 2.3606330638673825

さらに、コマンドラインでの計測も可能だ。

main関数の測定。コマンドオプションの-sは最初に1回だけ実行する文、-nは実行する文の回数、-rはタイマのリピート数になっている。

$ python -m timeit -n 1 -r 1 -s "import test_timeit" "test_timeit.main()" 12 3628800 1 loops, best of 1: 2.63 sec per loop

たらい回し関数(tarai)のみの測定。

$ python -m timeit -n 1 -r 1 -s "import test_timeit" "test_timeit.tarai(12, 6, 0)" 1 loops, best of 1: 2.63 sec per loop

階乗関数(factorial)のみの測定。100万回実行。

$ python -m timeit -n 1000000 -r 1 -s "import test_timeit" "test_timeit.factorial(10)" 1000000 loops, best of 1: 2.33 usec per loop

ところで、上に示した階乗関数については「車輪の再発明」なので、通常はmath.factorialを使った方が良いだろう。

続きを読む...

2013年1月1日火曜日

2013年と素因数分解

このブログ記事をはてなブックマークに追加

あけましておめでとうございます。

新年を迎えて、ふと、今年の西暦である2013を素因数分解したらどんな数に分解されるのか気になったので考えてみた。桁をすべて足し合わせると 2+0+1+3=6 と3の倍数になり、3で割り切れることは明白だ。3で割ると671となる。671に対し桁を一つ飛ばしにして足し合わせた数の差が (6+1)-7=0 で11の倍数(0を含む)なので11で割り切れることも簡単にわかる。671を11で割ると61の素数となる。つまり、2013の素因数分解は 3×11×61 となる。

大きな素数同士による合成数の素因数分解は非常に困難であることが知られている。この性質を利用して多くの暗号アルゴリズムで大きな素数による合成数が利用されている。現在のところ、素因数分解アルゴリズムとしては、楕円曲線法(ECM; Elliptic Curve Method)や複数多項式二次ふるい法(MPQS; Multiple Polynomial Quadratic Sieve)がよく使われているらしい。で、手軽にこれらを利用できるライブラリがないかと軽く調べたところNZMATH(ニジマス)という数論のためのPythonモジュールがあることを知ったので、早速インストールして使ってみた。

ところで、大きな素数といえばメルセンヌ素数が有名だが、マラン・メルセンヌは、2n-1 に対してnが257以下のとき、n = 2, 3, 5, 7, 13, 17, 19, 31, 67, 127, 257のときにのみ素数となると主張した。しかし、一部に間違いがあり、n = 61, 89, 107のときも素数であり、また、n = 67, 257は素数ではなく合成数であった。

そこで、誤って素数としていたn = 67のときの 147573952589676412927 について楕円曲線法(ECM)で解かせると、すぐに 193707721×761838257287 と結果が出力された。

>>> import nzmath.factor.methods as methods >>> methods.factor(147573952589676412927, method="ecm") [(193707721L, 1), (761838257287L, 1)]

更に、メルセンヌが見落としていた9番目と10番目のメルセンヌ素数である 2305843009213693951 (261-1) と 618970019642690137449562111 (289-1) を掛け合わせた46桁の合成数 1427247692705959880439315947500961989719490561 を複数多項式二次ふるい法(MPQS)で解いてみたところ、30分ほどで正しく2つの素数に分解できた。

>>> methods.factor(1427247692705959880439315947500961989719490561, method="mpqs") [(2305843009213693951L, 1), (618970019642690137449562111L, 1)]

ソースコードも簡単に確認できるし、NZMATHはなかなか楽しい。

それでは、本年もよろしくお願い申し上げます。

続きを読む...

2012年9月15日土曜日

「フカシギの数え方」の問題を解いてみた

このブログ記事をはてなブックマークに追加

先日、「『フカシギの数え方』 おねえさんといっしょ! みんなで数えてみよう!」という動画を見た。格子状のマスの左上から右下までの経路が何通りあるのかを調べて、格子が多くなればなるほど組み合わせの数が爆発的に増えることを教えてくれる動画だ。これは自己回避歩行(Self-avoiding walk)と呼ばれている問題らしい。

これだけ聞いてもそれほどインパクトはないのだが、動画に出てくるおねえさんの経路を調べあげる執念がもの凄く、ネット上でも結構な話題になっている。執念と言うよりも狂気に近い。しかし、話題になった割には動画内で言及されている高速なアルゴリズムを実装したという話を聞かなかったので、自分で確かめることにした。


動画のおねえさんは深さ優先探索によるプログラムを使っていると思われるが、それだとスパコンを使っても10×10マスの格子を解くのに25万年も掛かってしまう。そこで、高速化のためにゼロサプレス型二分決定グラフ(ZDD; Zero-Suppressed Binary Decision Diagram)と呼ばれるアルゴリズムを利用することにした。このアルゴリズムを開発したのは北大の湊先生で、ZDDによりすべての経路を見つけ出すアルゴリズムとしてクヌース先生のSIMPATHを使った。ZDDについてはクヌース先生も強い関心を持っていて、 The Art of Computer Programming Volume 4, Fascicle 1(TAOCP 4-1)ではBDD/ZDDの詳細な解説を読むことができる。演習問題の解説だけで書籍の半分を使っていることからしても気合の入れようがわかるだろう。

実際に自分のノートPCでZDDアルゴリズムを使ったコードを走らせたら、ほんの10秒程度で10×10マスの問題を解いてしまった。おねえさんがスパコンで25万年かかった問題をノートPCでたった10秒である。約8千億倍の高速化だ。これだけ劇的に変わるとやっぱり楽しい。そして、アルゴリズムの重要性を再認識させられた。

さて、以下におねえさんが利用したであろう自作の深さ優先探索(DFS)プログラムと高速な解法であるZDDアルゴリズムの両方を載せておく。ZDDについては4つのプログラムを1つのPythonスクリプトで統合している。これは、クヌース先生のSIMPATHおよびSIMPATH-REDUCEを利用しており、また、SGB (Stanford GraphBase)ライブラリでグラフを作成し、経路の数え上げにGMPを使ったC++で自作コードを組んだためである。最初はSIMPATHを参考に実装して一つのソースコードに統一しようかと思ったが、解説もちゃんと書かれているクヌース先生のコードをそのまま利用することにした。これらのプログラムの簡単な解説と使い方を載せておく。

まずはDFSによる実装を以下に示す。Node構造体から経路を辿って数え上げるだけの単純なプログラムである。


実行方法は以下の通り。引数の7は7×7のノードを表しており、マスで表現すれば6×6マスとなる。そして計算結果として、575,780,564通りの経路が探索されたことを示している。因みに7×7ノード(6×6マス)で約3分ほどの時間が掛かっており、それ以上の大きさだと現実的な時間で解くことが難しくなってくる。

% ./count_lattice_path_dfs 7 Count all paths from (0, 0) to (6, 6) in the 7x7 lattice graph: Count: 575780564

次に、ZDDを利用したプログラムのビルド方法を示す。

ビルドをする前にSGBおよびGMPを準備しておく。また、SIMPATHおよびSIMPATH-REDUCEについてはCWEBを利用しているので、それも入れておく。

最初に、SGBライブラリを用いて格子状のグラフを作成する。

% gcc gen_lattice.c -O3 -lgb -o gen_lattice

ソースコード内でグラフを下記のように作成している。board関数はチェスのような格子状のグラフを作成することができる。他にもグラフを作成する便利な関数が多数あるので、興味のある方はSGBライブラリのgb_basic.wあたりを参照して欲しい。

Graph* g = board(N, M, 0L, 0L, 1L, 0L, 0L);


次に、SIMPATHおよびSIMPATH-REDUCEを以下のようにビルドする。

% wget http://www-cs-faculty.stanford.edu/~uno/programs/simpath.w % ctangle simpath.w % gcc simpath.c -O3 -lgb -o simpath

% wget http://www-cs-faculty.stanford.edu/~uno/programs/simpath-reduce.w % ctangle simpath-reduce.w % gcc simpath-reduce.c -O3 -lgb -o simpath-reduce

因みに、ctangleでC言語への変換、cweaveでTeX形式への変換を行うことができるので、以下のようにTeX形式のドキュメントも一緒に作っておくことをお勧めする。また、これらのアルゴリズムを理解するために前述のTAOCP Vol. 4を読んでおくことが望ましい。

% cweave simpath.w % tex simpath.tex % dvipdf simpath.dvi # PDFファイルが必要なら. % cweave simpath-reduce.w % tex simpath-reduce.tex % dvipdf simpath-reduce.dvi # PDFファイルが必要なら.

最後に、ZDDから経路を数え上げるプログラムをコンパイルする。

% g++ count_path_mp.cpp -O3 -lgmp -o count_path_mp

経路を一つ一つ数えていては高速に数え上げるという点で意味がなくなってしまうので、ちゃんとメモ化再帰で数える。ZDDは高度に圧縮されているのでメモ化再帰でほとんど瞬時に数え上げることができる。


そして、上記のプログラムを統括するPythonスクリプトが以下となる。グルー言語としてもPythonは優秀だ。


実行方法は以下の通り。

% ./count_lattice_path.py 8 Count all paths from (0, 0) to (7, 7) in the 8x8 lattice graph: Count: 789360053252

DFSのプログラム同様、引数の8は8×8ノード(7×7マス)を表している。そして計算結果として、789,360,053,252通りの経路が探索されたことを示している。また、中間ファイルとして以下が出力される。

lattice_08.gb # SGBによるグラフデータ. lattice_08.out # SIMPATHの出力: not-necessarily-reduced BDD lattice_08.out.r # SIMPATH-REDUCEの出力: ZDD lattice_08.out.rc # SIMPATH-REDUCEの出力を再番付. lattice_08.log # 一連の実行ファイルのログ.

因みに、DFSで7×7ノード(6×6マス)を解くと3分ほどかかったが、ZDDでは0.1秒未満で完了した。また、ZDDで14×14ノード(13×13マス)は数分で終わることを確認したが、それ以上は途中で搭載メモリ以上に必要なメモリが大きくなるので実行していない。

続きを読む...

2012年5月14日月曜日

テトレーション・フラクタル

このブログ記事をはてなブックマークに追加

5年ほど前、このブログで虚数のテトレーションという記事を書いたことがある。テトレーション(tetration)とは、自らのべき乗を指定された回数反復する演算のことで、na と表現する。35 の場合、555 = 1.911×102185 となる。Pythonの関数で表現すれば以下のようになる。


以前の記事では、i が 0.43828+0.36059i に収束することを見つけたのだが、今回はそれに関連したフラクタルについて紹介したい。

テトレーション naa には複素数を指定することができる。このとき、a = x + yi として、それを複素平面に置く。ここで正の整数nを大きくしていき、発散と判定されたnに対応した色で平面を色分けする。発散しない場合、予め決めておいたnの上限値を使う。これで得られる図をテトレーション・フラクタル(tetration fractal)と呼ぶ。

n(x + yi)=a + biとしたとき、n+1(x + yi)=a' + b'i は以下のように計算できる。




Pythonでの実装はActiveSate CodeTetration Fractalとして載っている。今回はこれを少し修正したコード、C++で書いたコード、更にTBBで並列化したコードを用意して、それらの実行速度のベンチマークを取ってみた。

Python C++ C++ with TBB
258.732秒 16.999秒 0.976秒

ここで、複素平面の領域は(-1.5, 0)-(-0.75, 0.75)であり、最大繰り返し数は256としている。画像の大きさは1,024×1,024とした。実行環境はIntel Xeon X5650の2CPUで、12コア・24スレッドとなっている。

以下に、ソースコードを示す。詳細についてはコードを読んでいただきたい。

Pythonによる実装:

% ./tetration.py tetration.png -1.5 0.0 0.75 0.75 1024 1024



C++による実装:

% g++ tetration.cpp -o tetration -O3
% ./tetration tetration.out -1.5 0.0 0.75 0.75 1024 1024
% ./image_tetration.py tetration.out tetration.png



TBBを利用したC++による並列化実装:

% g++ tetration_tbb.cpp -o tetration_tbb -ltbb -O3
% ./tetration_tbb tetration.out -1.5 0.0 0.75 0.75 1024 1024
% ./image_tetration.py tetration.out tetration.png



C++の実行ファイルで生成された出力ファイルを画像ファイルに変換するスクリプトを以下に示す。


続きを読む...