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++の実行ファイルで生成された出力ファイルを画像ファイルに変換するスクリプトを以下に示す。


続きを読む...

2012年3月31日土曜日

IPythonの埋め込みプロットが素晴らしい

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

先日、Tokyo.SciPy #3に参加して、SciPyの生みの親であるTravis Olipant氏のセッションを聞いた。その時に最近のIPythonでは埋め込みプロットができることを知ったので早速入れてみたところ、これがとてもクールでカッコよかったのでここで紹介したい。

埋め込みプロットはIPythonのバージョン0.11からできるようになっている。ただ、自分が使っているLinux環境では普通に持ってくるとバージョン0.10が入ってしまうので使えない。そこで、IPythonの公式サイトから最新版のバージョン0.12を持ってきて入れてみた。

ただ、埋め込みプロットができるのはターミナル版ではなくQt版のシェルなので、それに関連するライブラリなども一緒に入れなくてはならない。IPythonの起動時に何々のライブラリが足りないとかいろいろと文句を言われるが、足りないものをyumなりapt-getなりソースコードをコンパイルなりして順次入れていけば使えるようになると思う。

あと、IPythonの設定ファイルの構成が0.10から大きく変わった。以前は .ipython/ 内の pythonrc を利用していたのだが、それがなくなって、代わりに .ipython/profile_default/ 内の ipython_config.py と ipython_qtconsole_config.py が使われるようになった。デフォルトの設定ファイルは以下のコマンドで作ることができる。

$ ipython profile create

そして、埋め込みプロットを使うためのQt版IPythonの起動コマンドは以下の通り。

$ ipython qtconsole --pylab=inline

これでQt版シェルが立ち上がって埋め込みプロットを表示できるようになる。

せっかくなので実際に埋め込みプロットを使ってみる。ここでは自分のお気に入りのマンデルブロ集合をプロットしてみた。因みに、シェル上での複数行のソースコードの編集も改良されたので入力がかなり楽になった。冒頭にスクリーンショットを貼っておく。

rangex = arange(-2.0, 1.0, 0.01) rangey = arange(-1.0, 1.0, 0.01) (X, Y), Z = meshgrid(rangex, rangey), [] for y in rangey: Z_ = [] for x in rangex: c, z = complex(x, y), 0.0 for i in range(100): if abs(z) > 2.0: break z = z * z + c Z_.append(i) Z.append(Z_) contourf(X, Y, Z, arange(100.0))

続きを読む...

2012年2月19日日曜日

Android端末をサーバにしてHTML5を使ったお絵かきBBSを作成する

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

最近、HTML5に触れる機会があり、その良さが何となく伝わってきたので、何かしら簡単なコードを書いてみたくなった。そこで、ちっちゃいけどリッチ、というギャップを楽しむためにAndroid端末を使うことにした。具体的には、Android端末をSL4Aウェブサーバにして、HTML5をインターフェイスにしたお絵かきBBSをPythonで書いてみた。

BBSでは絵とテキストを扱うことができて、それらはAndroid端末上でSQLiteのデータベースで管理される。利用者の利便を考えて、名前などはクッキーで保存する。3G回線や無線LANなどで接続することを考慮して、IPアドレスも取得できるようにした。また、書き込み時にサーバのAndroid端末が振動して書き込みがあったことを知らせてくれる。因みに、NTTドコモの3G回線で使うためには、グローバルIPが振られるmopera Uなどのサービスが必要で、spモードでは利用できない。

実際に作ってみたプログラムのスクリーンショットを冒頭に入れてみた。必要最低限の機能しかないが、それでも文章と画像を保存できるちゃんとした掲示板だ。草の根BBSのころを考えると隔世の感がある。

以下に今回作成したソースコードを載せておく。こんな短いコードでもちゃんと機能するのが不思議な感じだ。

image_bbs.py

# -*- coding: utf-8 -*- import sys,os,cgi,sqlite3,datetime import socket,fcntl from wsgiref.simple_server import make_server import android droid=android.Android() LIMIT=10 # 最大表示記事数. DB_FILE='/sdcard/image_bbs.sqlite' P=8080 con=sqlite3.connect(DB_FILE) cur=con.cursor() cur.execute('CREATE TABLE IF NOT EXISTS bbs (id INTEGER NOT NULL PRIMARY KEY AUTOINCREMENT, user TEXT, datetime TEXT, image TEXT, text TEXT)') INSERT_DB='INSERT INTO bbs VALUES(NULL,?,?,?,?)' def ipconfig(): s=socket.socket(socket.AF_INET,socket.SOCK_DGRAM) s.connect(('gmail.com',80)) return s.getsockname()[0] def post(user,image,text): if image=='' or text=='': return if user=='': user='匿名' dt=datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S') cur.execute(INSERT_DB,(cgi.escape(user.decode('utf-8')),dt,cgi.escape(image.decode('utf-8')).replace('\n','<br />'),cgi.escape(text.decode('utf-8')).replace('\n','<br />'))) con.commit() droid.vibrate() #droid.notify(user,text) def bbs(environ,start_response): global IP if environ['PATH_INFO']=='/': user='' if environ.has_key('HTTP_COOKIE'): cookie={} for d in environ['HTTP_COOKIE'].strip(';').split(';'): d=d.split('=') cookie[d[0]]=d[1] if cookie.has_key('BBSUSER'): user=cookie['BBSUSER'] if environ['REQUEST_METHOD']=='POST': fs=cgi.FieldStorage(fp=environ['wsgi.input'],environ=environ,keep_blank_values=1) user=fs.getfirst('user','').strip() image=fs.getfirst('bbsImage','') text=fs.getfirst('text','').strip() post(user,image,text) data=u"""<!DOCTYPE html> <html><head> <meta http-equiv="Set-Cookie" content="BBSUSER=%s"> <meta http-equiv="Content-Type" content="text/html; charset=utf-8" /> <title>Image BBS by Android</title></head><body> <p><a href="http://%s:%d">http://%s:%d</a></p> <form action="/" method="post" id="form"> <span>名前:<input type="text" name="user" size="20" maxlength="30" value="%s" /></span> <div><canvas id="bbsCanvas" width="400" height="200" style="border:1px solid;"></canvas></div> <textarea name="text" cols="40" rows="5"></textarea><br /> <input type="hidden" name="bbsImage" id="bbsImage" /> <span><input type="button" onclick="clearCanvas();" value="画像消去" /></span> <span><input type="button" value="更新" onclick="location.reload(true);" /></span> <span><input type="button" value="送信" onclick="setImage();this.form.submit();" /></span> </form> <script type="text/javascript"> var mouseDownFlag = false; window.onload = function() { draw(); setImage(); }; function draw() { var canvas = document.getElementById('bbsCanvas'); if (!canvas || !canvas.getContext) return false; var ctx = canvas.getContext('2d'); canvas.onmousemove = function(e) { if (mouseDownFlag) { var rect = e.target.getBoundingClientRect(); ctx.beginPath(); ctx.arc(e.clientX - rect.left, e.clientY - rect.top, 3, 0, Math.PI * 2, false); ctx.fill(); } } canvas.onmousedown = function(e) { mouseDownFlag = true; } canvas.onmouseup = function(e) { mouseDownFlag = false; } canvas.onmouseout = function(e) { mouseDownFlag = false; } } function setImage() { var canvas = document.getElementById('bbsCanvas'); if (!canvas || !canvas.getContext) return false; document.getElementById("bbsImage").value = canvas.toDataURL(); } function clearCanvas() { var canvas = document.getElementById('bbsCanvas'); if (!canvas || !canvas.getContext) return false; var ctx = canvas.getContext('2d'); ctx.clearRect(0, 0, 400, 200); setImage(); } </script>""" % (user,IP,P,IP,P,user) cur.execute('SELECT * FROM bbs ORDER BY id DESC') for i, row in enumerate(cur): if i>=LIMIT: break data+=('<div><p>%d <b>%s</b> %s</p><img src="%s" alt="Image" width="400" height="200" style="border:1px solid;" /><p>%s</p></div>' % row) data+='</body></html>' start_response('200 OK',[('Content-type','text/html;charset=utf-8')]) return [data.encode('utf-8')] IP=ipconfig() httpd=make_server('',P,bbs) httpd.serve_forever()

続きを読む...