tag:blogger.com,1999:blog-61533010783353369922024-03-05T17:34:33.577+09:00良いもの。悪いもの。Good things. Bad things.noxhttp://www.blogger.com/profile/01142433117618538300noreply@blogger.comBlogger311125tag:blogger.com,1999:blog-6153301078335336992.post-66354847707266306252020-08-05T00:23:00.001+09:002020-08-05T00:28:07.655+09:00Thunderbirdで標準搭載されたPGPによる暗号化を利用する電子メールクライアントのMozilla Thunderbird 78でやっとOpenPGPの機能が標準搭載された。<br /><br />
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiWkfdvFYmyeRZmhymGFTYk0rYluzEz4VJITBnjQwicXO70Pths-RI60KhZQKW3IAEnPuufj0MmRpUaYG67SMLPrhfRugRVk6j3LiIDvU0MqqLz-xEWa0ZOG5fextbF2iwvmhtZoz13sAy5/s512/Mozilla_Thunderbird_logo.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="512" data-original-width="512" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiWkfdvFYmyeRZmhymGFTYk0rYluzEz4VJITBnjQwicXO70Pths-RI60KhZQKW3IAEnPuufj0MmRpUaYG67SMLPrhfRugRVk6j3LiIDvU0MqqLz-xEWa0ZOG5fextbF2iwvmhtZoz13sAy5/s0/Mozilla_Thunderbird_logo.png" /></a></div><br />
とはいっても、78.0.1ではまだ不完全で、78.1.0でほぼ完成したようだが、デフォルトではOpenPGPは使えない。今後リリース予定の78.2からデフォルトで動作するようになるらしい。<br /><br />
78.1.0で使うためには、「オプション」の「一般」の一番下の「設定エディター...」を開き、 <code>mail.openpgp.enable</code> を <code>true</code> にする必要がある。
ただし、設定をいじると動作保証対象外になるので、自己責任で対応すること。
今回の動作検証ではWindows 10上で実行させた<a href="https://ftp.mozilla.org/pub/thunderbird/releases/78.1.0/win64/ja/">64ビット版Thunderbird 78.1.0</a>を使っている。<br /><br />
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjYNWNtnZvlzIYmAhS0a4D7BnfRgT2nv_an8lu0u_n9yzANIXhCFM2OagZXd6c7ykLUIVke2-hnrY0jSNwc4uUrA_8QONP5T2p814N95IFKmRFvWQxng-KnPN0MwZCLvUA3ZYa7X1OEKXJb/s796/thunderbird_pgp_01.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="587" data-original-width="796" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjYNWNtnZvlzIYmAhS0a4D7BnfRgT2nv_an8lu0u_n9yzANIXhCFM2OagZXd6c7ykLUIVke2-hnrY0jSNwc4uUrA_8QONP5T2p814N95IFKmRFvWQxng-KnPN0MwZCLvUA3ZYa7X1OEKXJb/s640/thunderbird_pgp_01.png" width="640" /></a></div><br />
設定を済ませたら、「エンドツーエンド暗号化」のOpenPGPの「Add Key...」で秘密鍵をインポートもしくは作成する。
<a href="https://gnupg.org/download/">GnuPG</a>などですでに鍵を作っているのであれば、それをエクスポートして読み込ませるだけだ。<br /><br />
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhMg1h7_fuSoe_1tNRmsHnLF2BaWFe3OXk99yqZz3oBdj_vbQfJppV7eKpdNeQBjKwJrXgdjWT5yZCmsK_oW7waR7HCv8ci4RxVHrOeMbyEloNyVDyTzX6hd-lJETVlG3UsNmBPA6BRfWNM/s1191/thunderbird_pgp_02.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="679" data-original-width="1191" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhMg1h7_fuSoe_1tNRmsHnLF2BaWFe3OXk99yqZz3oBdj_vbQfJppV7eKpdNeQBjKwJrXgdjWT5yZCmsK_oW7waR7HCv8ci4RxVHrOeMbyEloNyVDyTzX6hd-lJETVlG3UsNmBPA6BRfWNM/s640/thunderbird_pgp_02.png" width="640" /></a></div><br />
公開鍵をインポートする場合、「OpenPGP Key Manager」の「Import Public Key From File」で対応可能だ。<br /><br />
早速、メールを作成し、「セキュリティ」から「暗号化が必要」と「このメッセージにデジタル署名」を選択して送信する。<br /><br />
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhlwwyyhxAka8wAXQWiCm3WAy3Jl3lBABn3rL8zEoaxVtLfSdlX8NXHq_XL0dGN8apkCqSWsx5afThZ935otQGH4E4zCXH4RwfEgH-EFEdVRsYbxpe3rFYLoorCS2tfqvBUC9idf7RYxoM-/s861/thunderbird_pgp_03.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="831" data-original-width="861" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhlwwyyhxAka8wAXQWiCm3WAy3Jl3lBABn3rL8zEoaxVtLfSdlX8NXHq_XL0dGN8apkCqSWsx5afThZ935otQGH4E4zCXH4RwfEgH-EFEdVRsYbxpe3rFYLoorCS2tfqvBUC9idf7RYxoM-/s640/thunderbird_pgp_03.png" width="640" /></a></div><br />
メッセージは自分自身に送っているのであまり意味はないが、それぞれの機能がきちんと動作していることを確認することはできる。<br /><br />
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjQYilJzcvlUsbxSQe96ipcUhmKnp5JThzjHv0I6irigW9E8o57sDAzbCnbvluvDMChgfoilg9bJkQfwgnFBT40Y9aEVt4LZPDsa1eSXCxb-dtHFEYRq3Y6-L9x32qb60JJYXNPW8VaQ_Rr/s1754/thunderbird_pgp_04.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="863" data-original-width="1754" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjQYilJzcvlUsbxSQe96ipcUhmKnp5JThzjHv0I6irigW9E8o57sDAzbCnbvluvDMChgfoilg9bJkQfwgnFBT40Y9aEVt4LZPDsa1eSXCxb-dtHFEYRq3Y6-L9x32qb60JJYXNPW8VaQ_Rr/s640/thunderbird_pgp_04.png" width="640" /></a></div><br />
この通り、暗号化とデジタル署名ができていることを確認できた。<br /><br />
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjnXA6NLSvH534Hf8YPCtWpu7v_eF3kgbhd79zYEYpgs8rJlkXpwWTrF2N5HGgQJlgW0P0-aryiSWCheygRYnfZdkboIVdRmJ-LJUawWZjH536jcexNSiyKoVmJErBmoFyrv7QHs5oHC1n9/s602/thunderbird_pgp_05.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="392" data-original-width="602" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjnXA6NLSvH534Hf8YPCtWpu7v_eF3kgbhd79zYEYpgs8rJlkXpwWTrF2N5HGgQJlgW0P0-aryiSWCheygRYnfZdkboIVdRmJ-LJUawWZjH536jcexNSiyKoVmJErBmoFyrv7QHs5oHC1n9/s0/thunderbird_pgp_05.png" /></a></div><br />
PGPは昔からあるが、自分の周りでは使われているのをほとんど見たことがない。
Thunderbirdに標準搭載されたのを機に、広まってくれないだろうか。<br />
noxhttp://www.blogger.com/profile/01142433117618538300noreply@blogger.com0tag:blogger.com,1999:blog-6153301078335336992.post-20934169004463926192020-08-02T20:45:00.000+09:002020-08-02T20:45:19.023+09:00古典的で未来的な人にやさしい顔映像匿名化現代では簡単に写真や動画を撮ることができる。それは、個人情報の流出と直結する。
特に動画に映った人物について、匿名化することはなかなかに大変だ。
目線やモザイク、塗りつぶしなどで顔を隠すケースをよく見かけるが、個人的にはちょっと微妙。犯罪臭がするし、淫靡な印象を与えてしまう。
ほかにもアイコンやシンボルで顔を隠したりすることができるツールなどもあるが、ちょっと自分好みではないし、何よりオリジナリティに欠ける。<br /><br />
そこで、古典的でもあり、未来的でもある、そして新しい、そんな顔映像匿名化を考えてみた。<br /><br />
匿名化するにあたって、まずは顔、特に目と鼻を検出できなくてはならない。
幸いなことに、<a href="https://pytorch.org/docs/stable/torchvision/models.html#keypoint-r-cnn">Keypoint R-CNN</a>モデルを使えばそれが可能になる。
利用するだけならば、PyTorchのチュートリアルレベルの知識を持っていれば十分。<br /><br />
そして、コードはできるだけシンプルで分かりやすく、その場ですぐに書けてしまえることが重要だ。
これは大事なことで、たとえ面白そうな技術でも、それを体験できるまでに時間がかかれば多くの人は興味を失ってしまう。
体験して理解さえすれば、達成感と面白さを実感できるし、次を目指す動機づけになる。<br /><br />
閑話休題。さて、今回作るのは、動画の顔画像の匿名化だ。そして、古典的で未来的、何のことだと思われるかもしれないが、使っている技術は新しく、できたものは懐かしいということ。
具体的には、顔画像に細工をして「へのへのもへじ」にしてしまうのだ。とはいっても、目と鼻を置き換えるだけなので「のもの」を使う。
江戸時代中期にはよく知られていたという「へのへのもへじ」は、まさに古典的で懐かしいと言えるだろう。<br /><br />
作成したコードは <a href="https://github.com/foota/nomono-anonymizer">foota / nomono-anonymizer</a> にある。使い方は次の通り。<br /><br />
<code>$ python nomono-anonymizer.py [入力動画ファイル名] [出力動画ファイル名]</code><br /><br />
動画はフレームごとに画像データとして処理され、それを事前学習した<a href="https://pytorch.org/docs/stable/torchvision/models.html#keypoint-r-cnn">Keypoint R-CNN</a>のモデルに推論をさせる。
モデルは画像データに対し人物領域を検出し、それに対して目、鼻、耳、肩、肘など、17のキーポイントを推定する。
今回は、右目・左目・鼻の3つのキーポイントを利用し、取得した座標に対して、画像データ上に「のもの」の文字を書き込むだけだ。
あとは、それを出力画像として動画としてエンコードすれば完了となる。<br /><br />
出来上がった動画は以下の通り。使った入力動画はPixabayの「<a href="https://pixabay.com/videos/id-6582/">イタリアの人々</a>」(1280x720/19秒)と「<a href="https://pixabay.com/videos/id-34418/">食事中の客</a>」(1280x720/13秒)の2つだ。<br /><br />
<strong>イタリアの人々</strong><br />
<div align="center"><iframe width="560" height="315" src="https://www.youtube.com/embed/4xwN4jtAcJw" frameborder="0" allow="accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe></div><br />
<strong>食事中の客</strong><br />
<div align="center"><iframe width="560" height="315" src="https://www.youtube.com/embed/vuQ1YnijqoQ" frameborder="0" allow="accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe></div><br />
見てわかる通り、これだけ御託を並べて結局やったことは動画に落書きしただけ。まあ、楽しかったので良いのだけど。<br /><br />
最後にCPU/GPUによる実行時間の比較を載せておく。GPUを使わないとなかなか厳しい処理速度だ。<br /><br />
<div align="center">
<table border="1" width="60%">
<caption>
実行時間の比較
</caption>
<tr align="center">
<th>動画</th>
<th>CPU / GPU</th>
<th>実行時間 (秒)</th>
</tr>
<tr>
<td align="center" rowspan="2">イタリアの人々<br /><a href="https://pixabay.com/videos/id-6582/">入力</a> | <a href="https://www.youtube.com/embed/4xwN4jtAcJw">出力</a></td>
<td align="center">CPU</td>
<td align="right">1646.350</td>
</tr>
<tr>
<td align="center">GPU</td>
<td align="right">75.423</td>
</tr>
<tr>
<td align="center" rowspan="2">食事中の客<br /><a href="https://pixabay.com/videos/id-34418/">入力</a> | <a href="https://www.youtube.com/embed/vuQ1YnijqoQ">出力</a></td>
<td align="center">CPU</td>
<td align="right">1325.606</td>
</tr>
<tr>
<td align="center">GPU</td>
<td align="right">60.313</td>
</tr>
</table>
</div><br />
noxhttp://www.blogger.com/profile/01142433117618538300noreply@blogger.com0tag:blogger.com,1999:blog-6153301078335336992.post-89856813077603683542020-07-23T15:51:00.006+09:002020-07-29T05:19:19.636+09:00人工知能が見ている世界今の世の中、そこら中に人工知能(AI)と呼ばれるシステムがあふれている。特にディープラーニングによるニューラルネットワークモデルを使ったシステムは高度な処理が可能であると広く考えられているようだ。では、それを使ったAIはどのようにこの世界を見ているのだろうか。<br /><br />
ディープラーニングで得られる特徴を含んだ値(特徴マップ)がAIの視覚と呼べるだろう。AIについて扱った教材では、ディープラーニングにより学習されたニューラルネットワークモデルの特徴マップを例示していることがあるが、特定の画像に対する特定の特徴の一部を載せているだけなので、AIがどのように世界を見ているのか、それだけで判断するのは難しい。<br /><br />
そこで、任意の動画に対してニューラルネットワークモデルを使って特徴マップを映像化することにした。今回、対象とするネットワークアーキテクチャは、<a href="https://pytorch.org/hub/pytorch_vision_alexnet/">AlexNet</a>と<a href="https://pytorch.org/hub/pytorch_vision_vgg/">VGG19</a>だ。それぞれ2012年と2014に公表されており、どちらのアーキテクチャも<a href="http://www.image-net.org/">ImageNet</a>の分類できわめて高い性能を誇っている。<br /><br />
特徴マップの映像化を行うにあたって、まずは動画処理のために<a href="https://github.com/kkroening/ffmpeg-python">ffmpeg-python</a>を利用することにした。これは、名前からもわかるようにPythonのFFmpegラッパーで、動画を扱うのにはかなり便利なライブラリとなっている。次に、AlexNetとVGG19を利用するために<a href="https://pytorch.org/">PyTorch</a>とそのパッケージである<a href="https://pytorch.org/docs/stable/torchvision/index.html">torchvision</a>を利用した。torchvisionにはAlexNetとVGG19の学習済みネットワークモデルも用意されている。使い方は簡単で、例えばAlexNetなら、以下のコードでネットワークモデルをすぐに使うことができる。<br /><br />
<code>model = torchvision.models.alexnet(pretrained=True)</code><br /><br />
今回は学習をせず用意されているモデルを使うだけなので、推論モードにして勾配計算をしないように <code>model.eval()</code> と <code>torch.set_grad_enabled(False)</code> を記述しておく。<br /><br />
処理の流れは次の通り。入力データとして適当な動画ファイルを用意する。動画ファイルが読み込まれると、ffmpegによりフレーム画像を一枚ずつ読み出し、それをnumpyで配列に変換。変換されたデータをPyTorchのモデルに渡す。前処理としてtransposeやunsqueezeでモデルが処理できるような形に変換しておく。モデル内での計算は特徴のみを抽出するため、 <code>model.features()</code> を呼び出す。続けて <code>model.avgpool()</code> を呼び出してもいいが、プーリングによりデータが圧縮されて出力画像の解像度が落ちる。<br /><br />
モデルによる出力データとしてチャネル数分の特徴マップが出てくるので、それを出力動画の画面に並べることにする。AlexNetであればチャネル数は256なので16×16に並べる(VGG19なら512なので32×16)。まず、モデルからの出力を(Ho×Wo,
256)の形に変換して(Ho,
Woはそれぞれ出力された特徴マップの高さと幅を表す)、人に見やすく視覚化するために、特徴マップごとに正規化(min-max
scaling;
0-255)を行っている。そして、それを縦横16×16枚の特徴マップの画像として並ぶようにデータを変換し、画像処理ライブラリである<a href="https://pillow.readthedocs.io/en/stable/">PIL (pillow)</a>を用いて元のサイズにスケーリングする。さらに、オリジナルの画像を透過度(アルファ値)0.3にして重ね合わせて、動画フレームのための画像データを生成する。最後にffmpegにより動画ファイルにエンコードして完了となる。<br /><br />
コードはGitHubリポジトリの<a
href="https://github.com/foota/neural-networks-view"
>foota / neural-networks-view</a
>に置いてある。使い方は以下の通り。<br /><br />
<code>python nnview.py 入力動画ファイル名 出力動画ファイル名</code><br /><br />
以下に示す動画がAlexNetが見ている世界になる。256個もの眼を持って、矯めつ眇めつそれを眺め、何であるのかを認識している。人が見て何となく認識できる映像もあれば、ただの点滅にしか見えないものまで様々だ。<br /><br />
<div align="center">
<iframe width="560" height="315" src="https://www.youtube.com/embed/LBSlxE_b38s" frameborder="0" allow="accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe></div><br />
同様に、VGG19の見ている世界が以下の映像だ。AlexNetの倍の眼を持ち、512個もの異なった視点を全結合層で統合して判断を下していると考えるとなんだか不思議な気分になる。入力映像として<a href="https://pixabay.com/videos/">Pixabay</a>の<a href="https://pixabay.com/videos/id-39264/">ハトの動画</a> [1280x720]を使わせてもらった。<br /><br />
<div align="center">
<iframe width="560" height="315" src="https://www.youtube.com/embed/-sbFQaa77W4" frameborder="0" allow="accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe></div><br />
最後に、手元の実行環境(CPU: Intel Core i9 10900K @3.70GHz / GPU: NVIDIA GeForce RTX 2080Ti)での速度比較を示す。CPUとGPUの速度差は、AlexNetで4.4倍、VGG19で24.8倍であった。この違いが生じている理由は、AlexNetと比較して計算量が多いVGG19の方が効果的にGPUを利用できているためだ。<br /><br />
<div align="center">
<table border="1" width="80%">
<caption>
実行時間の比較 (入力データ: <a href="https://pixabay.com/videos/id-39264/">ハトの動画</a> [1280x720])
</caption>
<tr align="center">
<th>アーキテクチャ</th>
<th>CPU / GPU</th>
<th>実行時間 (秒)</th>
</tr>
<tr>
<td align="center" rowspan="2">AlexNet</td>
<td align="center">CPU</td>
<td align="right">44.204</td>
</tr>
<tr>
<td align="center">GPU</td>
<td align="right">9.998</td>
</tr>
<tr>
<td align="center" rowspan="2">VGG19</td>
<td align="center">CPU</td>
<td align="right">550.690</td>
</tr>
<tr>
<td align="center">GPU</td>
<td align="right">22.204</td>
</tr>
</table>
</div><br />
今回はAlexNetとVGG19の2つのアーキテクチャについて動画を作成してみたが、特徴マップを抽出できればどのアーキテクチャにも適用できると思うので、興味があれば試して欲しい。また、いろいろな動画で試してみることで、ニューラルネットワークしか知り得ない、面白い特徴が見つかるかもしれない。<br /><br />
※ 入力フレーム画像の正規化が入っていなかったのでコード修正・動画差し替え (2020-07-29)<br />
noxhttp://www.blogger.com/profile/01142433117618538300noreply@blogger.com0tag:blogger.com,1999:blog-6153301078335336992.post-73076462622916758142020-07-12T13:55:00.003+09:002020-07-24T05:37:20.959+09:0012年前の非接触型ICカードリーダPaSoRiで遊ぶ最近ずっとリモートワークだったので、これを機に書斎の整理整頓を行った。段ボールに無造作に放り込まれたPC関連機器や電子ガジェットを漁っていると、2004年10月に発売されたソニーの非接触型ICカードリーダPaSoRi RC-S320が出てきた。2008年3月に購入したので実に12年以上前のモノだ。<br /><br />
<div class="separator" style="clear: both; text-align: center;"><a target="_blank" href="https://www.amazon.co.jp/gp/product/B0009YVAW4/ref=as_li_tl?ie=UTF8&camp=247&creative=1211&creativeASIN=B0009YVAW4&linkCode=as2&tag=footawwwservi-22&linkId=2f1ee5ca19cb45df5c9cc32a09282d0c"><img border="0" src="//ws-fe.amazon-adsystem.com/widgets/q?_encoding=UTF8&MarketPlace=JP&ASIN=B0009YVAW4&ServiceVersion=20070822&ID=AsinImage&WS=1&Format=_SL250_&tag=footawwwservi-22" ></a><img src="//ir-jp.amazon-adsystem.com/e/ir?t=footawwwservi-22&l=am2&o=9&a=B0009YVAW4" width="1" height="1" border="0" alt="" style="border:none !important; margin:0px !important;" /></div><br /><br />
このまま捨て置くのは勿体ないし、ちょっと興味をそそられたので、<a href=" http://handasse.blogspot.com/2008/03/python-pasorisuica.html">以前のブログ記事</a>に上げた<a href=" http://handasse.blogspot.com/2008/04/python-pasorisuica.html">コード</a>がちゃんと使えるか、最近新調したWindows 10搭載PCで確認してみた。<br /><br />
当然、ドライバが入っていないのでそのままでは動かない。早速最新のドライバを入れようと思ったが、実はRC-S320はとっくに販売終了していて、新しいPaSoRiとはドライバが違う。まあ、そりゃそうか。調べたところ古いドライバ(<a href=" https://www.sony.co.jp/Products/felica/consumer/download/old2_felicaportsoftware.html">NFCポートソフトウェア Version 5.4.8.6</a>)を落としてくれば動くことが分かった。動作確認環境はWindows 8.1までみたいだけど、まあ、大丈夫でしょ。<br /><br />
ドライバを入れたところ、デバイスマネージャーでちゃんと認識された。よし、これで準備OKだ。動かしてみよう! …いやまてよ? 今時Python2は許されない。ということでPython3のコードにちょこちょこっと修正して、これで良し。PaSoRiの上にモバイルSuica機能を搭載したAndroid端末を置いて、さっそく実行 <code class="black">python read_felica.py</code> (カタカタカタッターン!)<br /><br />
<code>OSError: [WinError 193] %1 は有効な Win32 アプリケーションではありません。</code><br /><br />
あ、これよく見るやつだ。実行パス名にスペースとか入ってちゃんと認識しないやつだ。と思ったけど調べてみるとパス名は問題ないみたい。あとは、そうそう、64ビット環境で32ビットのライブラリ使ったとか。それだ! 実行環境はPython3の64ビット版なので読み込むDLLも64ビットじゃないと。でも、利用させてもらった <a href="http://felicalib.tmurakam.org/">felicalib</a> は32ビット版しかないみたい。<a href="https://ja.osdn.net/projects/felicalib/wiki/x64">0.4.3から64ビット版にも対応する</a>みたいだけど、最新版の0.4.2が2008年6月リリースでそれ以降出ていないということはもう出ない可能性が高い。ちなみにソニーのドライバは当時64ビット版が存在しなかったみたいだけど、今は入っているようだ。<br /><br />
それならとりあえず32ビット版Pythonで動くか確認してみよう。ちなみに、Pythonの64ビット版と32ビット版の共存は簡単。<a href="https://www.anaconda.com/products/individual#Downloads">Anaconda (64ビット版)</a>を入れているので、以下のように環境を作ってしまえば良い。<br /><br />
<code>> set CONDA_FORCE_32BIT=1
> conda create -n py37_win32 python=3.7
> set CONDA_FORCE_32BIT=</code><br /><br />
利用するときは <code class="black">conda activate py32_win32</code> とすれば32ビット版環境が立ち上がる。<br /><br />
さっそく実行してみると問題なく動いた。めでたしめでたし。…とは思えないのがエンジニアの性(さが)。普段64ビット環境を利用していて、RC-S320用の64ビット版ドライバが用意されていて、32ビット版環境で妥協するのは許されない感じ。なので、まずはfelicalib以外で対応されていないか確認してみた。<br /><br />
PythonでNFCタグを読み書きするのに<a href=" https://nfcpy.readthedocs.io/en/latest/index.html">nfcpy</a>というライブラリが近年ではよく使われているらしい。ただ、前準備も必要で、<a href="">Zadig</a>を使ってNFC Port/PaSoRiをWinUSBとして再インストールし、<a href=" https://libusb.info/">libusb</a>を入れて初めて使えるようになる。<br /><br />
それでnfcpyを使ってみたけど、認識しない。何故だとよくよく調べてみると、古いPaSoRiであるRC-S320には対応していないらしい。うーん、無駄骨か。動作したとしても、一般ユーザにプログラムを利用してもらう場合、libusbを入れること自体がそもそもハードルが高い。libusbで使えるようにするとNFCポートソフトウェア付属のNFCポート自己診断も当然ながら失敗するし。結局、「ドライバーを元に戻す」でlibusbは使わないことにした。<br /><br />
なら、felicalibを64ビットネイティブでビルドするしかない。幸いソースコードとVS2008用のプロジェクトファイルは付属しているので、それをVS2019で開いてみる。適当に変換されたのでビルドしてみたらコンパイルエラー。<br /><br />
<code>RC1015 cannot open include file 'afxres.h'</code><br /><br />
どうやらMFC用ヘッダファイルの afxres.h がないということらしい。ならこれを winres.h に変更してリビルド。普通に通った。次にプラットフォームのWin32からx64を追加してビルド。問題なく64ビット用DLLが生成された。念のためdumpbinコマンドでドライバのDLLの中身を確認したけど、64ビット版と32ビット版で構成は同一だった。<br /><br />
64ビット用DLLが用意できたので、64ビットPython環境でread_felica.pyを実行してみる。そうすると型が合わないよというエラーが出て実行できない。なるほど、ctypesを使っているので、ちゃんと argtypesで指定しないとね。きちんと指定したところ、実行することができた。初めから64ビット環境用にfelicalibをビルドすればよかった。<br /><br />
というわけで、12年前のPaSoRiでも現役で利用できた。実行するにあたって、64ビット版felicalibをパスの通った場所に置く。そして、サイバネ駅コードが入ったCSVファイル(stationcode.csv)を実行ディレクトリに置く。stationcode.csvは、<a href=" http://www014.upp.so-net.ne.jp/SFCardFan/">IC SFCard Fan</a>の<a href=" http://www.denno.net/SFCardFan/index.php">路線・駅コード一覧・登録</a>からダウンロードしてきたExcelファイルをCSVフォーマット(UTF-8)に変換したファイルだけど、現在のところダウンロード時にタイムアウトが出てデータが落とせなくなっているようだ。南無三。仕方がないので、Webスクレイピング用のコードを書いて対処した。<br /><br />
最後に、作成したコードを以下のリポジトリに上げておく。興味があればいじってみてほしい。<br /><br />
<a href="https://github.com/foota/felica-pasori">https://github.com/foota/felica-pasori</a><br />noxhttp://www.blogger.com/profile/01142433117618538300noreply@blogger.com0tag:blogger.com,1999:blog-6153301078335336992.post-7749543293202762022017-06-08T19:11:00.000+09:002017-06-08T19:14:59.716+09:00PythonでSlackのボットを作成する職場でSlackを利用している。なかなか使いやすいし便利だ。しかし、技術者としては与えられた機能を利用するだけではやはり満足できない。自分好みの道具を作りたい。<br /><br />
というわけで、単なるチャットツールとして使うSlackでは物足りなくなったので、Slackのボット(bot)をPythonで作ってみることにした。Hubotのようなボット用フレームワークも利用可能だが、SlackのAPIをガシガシ制御する感じにしたかったので、<a href="https://api.slack.com/rtm">Slack Real Time Messaging API</a>をPythonから制御してSlackと連携した。今回は、書き込んだ数字を素因数分解するボットを作ってみた。<br />
<span id="fullpost"></span>
<div id="fullpost"><br />
まずは、Slack用のボットを利用する準備をする。<a href="https://api.slack.com/">Slack API</a>にアクセスして、<a href="https://api.slack.com/bot-users">Bot users</a>をクリックし、Custom bot usersにある<a href="https://my.slack.com/services/new/bot">creating a new bot user</a>のリンクからBotsのページを開く。
名前(ここでは @prime_decomposition)を入力したらAdd bot integrationボタンをクリックする。<br /><br />
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiw2SEXi0IO82vwJX_Mxo49tgi_5ZKyuq3dAPzX3WHvKA0UQe-l3qRwlh8Zm8uMj2OaB1O9Gs0zUoSWhXE07QB18yIbKPblmFeDaBfvz5WfZ3UjMMwssiRbRyVVxaGDpU1iYKxpwCIAHRrx/s1600/slack_bots_01f.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiw2SEXi0IO82vwJX_Mxo49tgi_5ZKyuq3dAPzX3WHvKA0UQe-l3qRwlh8Zm8uMj2OaB1O9Gs0zUoSWhXE07QB18yIbKPblmFeDaBfvz5WfZ3UjMMwssiRbRyVVxaGDpU1iYKxpwCIAHRrx/s400/slack_bots_01f.png" width="400" height="229" data-original-width="1154" data-original-height="661" /></a></div><br /><br />
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEghpPt8y8IiM5jXXKCz7HybVDGng4vpNS4-9kTWISjzyoXs9jNqugodf5iaUaKwvUw4dbB4TfJZdbpD6209U9AzlYB-HFn2RGHSNef2TYLLx9zL4eeObe7yuo_1wnQwUZKQoTWN_TpoLsmZ/s1600/slack_bots_02f.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEghpPt8y8IiM5jXXKCz7HybVDGng4vpNS4-9kTWISjzyoXs9jNqugodf5iaUaKwvUw4dbB4TfJZdbpD6209U9AzlYB-HFn2RGHSNef2TYLLx9zL4eeObe7yuo_1wnQwUZKQoTWN_TpoLsmZ/s400/slack_bots_02f.png" width="400" height="230" data-original-width="1157" data-original-height="665" /></a></div><br /><br />
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhzxuTpKcgwLYJHWZd2lTymwfdCO2OBFA8QvLBTDXVnJZoFdM3_qKB-8eIjtZJpwThs3swcpWeyYzH1nzs25ANuNYkQeBoCU54yM-Aiw0dK9uPi_QoTe4NI53px3iWu9ca4_xMgoHRHWRbY/s1600/slack_bots_03f.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhzxuTpKcgwLYJHWZd2lTymwfdCO2OBFA8QvLBTDXVnJZoFdM3_qKB-8eIjtZJpwThs3swcpWeyYzH1nzs25ANuNYkQeBoCU54yM-Aiw0dK9uPi_QoTe4NI53px3iWu9ca4_xMgoHRHWRbY/s400/slack_bots_03f.png" width="400" height="233" data-original-width="1155" data-original-height="674" /></a></div><br /><br />
Integration Settingsの画面ではAPI Tokenが表示されており、ボットのスクリプトでこれを利用する。あとは好みに応じてアイコンの変更やボットの説明を記述しておく。今回のボットでは、デフォルトのフェイスマークに用意してあるロボットのアイコンを使い、説明文は「素因数分解をするよ!」とした。あとは、Save Integrationのボタンをクリックして設定を保存する。<br /><br />
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhHfGV7IqKAXOqOYuYvnTGng2Y1gl0IZxO7NuoCK2-vs5ck4obOE9yi_341p49hpKtNG1yDXJmZPm8QYxL50mfOinEOsvRg2uYRk5SWwM4jgYIlb_1EwOvQ-hOQBKIJvzr4CPgPlaKfddpE/s1600/slack_bots_04a.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhHfGV7IqKAXOqOYuYvnTGng2Y1gl0IZxO7NuoCK2-vs5ck4obOE9yi_341p49hpKtNG1yDXJmZPm8QYxL50mfOinEOsvRg2uYRk5SWwM4jgYIlb_1EwOvQ-hOQBKIJvzr4CPgPlaKfddpE/s400/slack_bots_04a.png" width="400" height="233" data-original-width="1155" data-original-height="673" /></a></div><br /><br />
ボットの準備ができたら、それをSlackへメンバーとして加える。ボットを参加させたいチャンネル(ここの画像では #general だが #random に参加させたい場合は #random を選択すること)のChannel Settings(歯車のアイコン)から、Invite team members to join...を選択する。先ほど設定したボットがメンバーとして表示されているのでそれを選ぶ。これで下準備は完了だ。<br /><br />
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhnLyUw22osyVO-U6m8Hpr8MC9NZd-3DQAs9XLz1g8r1xMqW6fwFEXaSZ4GkhXkUnZd06k7mubWVqcaYxPSucrnzOay5S7HwstFjbQtP9GMUzJfu9YZ8YD-U7CttsUP6TMw8uJiR6WdtSGb/s1600/slack_bot_slide_05b.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhnLyUw22osyVO-U6m8Hpr8MC9NZd-3DQAs9XLz1g8r1xMqW6fwFEXaSZ4GkhXkUnZd06k7mubWVqcaYxPSucrnzOay5S7HwstFjbQtP9GMUzJfu9YZ8YD-U7CttsUP6TMw8uJiR6WdtSGb/s400/slack_bot_slide_05b.png" width="400" height="133" data-original-width="1221" data-original-height="406" /></a></div><br /><br />
次にPythonを使ってボット用スクリプトを作成する。目的は、Slack上で誰かが数字を書き込んだ際に、それを取得して素因数分解をすること。そのためには、Slackに書き込まれた内容を取得しなければならないが、そのためにslackclientというPythonのモジュールを利用する。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;"># pip install slackclient</span><br /><br />
pipであれば上記のコマンドでインストールが終わる。手動で導入する場合は、<a href="https://github.com/slackhq/python-slackclient">python-slackclient</a>から入れる。<br /><br />
また、素因数分解プログラムをスクラッチから作るのも面倒なので、数論用ライブラリ<a href="http://tnt.math.se.tmu.ac.jp/nzmath/">NZMATH(ニジマス)</a>を利用する。これもpipでインストールできる。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;"># pip install nzmath</span><br /><br />
さて、ここまでできたらあとはコードを書くだけだ。Slackからデータを取得する場合は SlackClient.rtm_read() を、書き込みには SlackClient.rtm_send_message([channel, message]) を使う。これらを使うためには SlackClient.rtm_connect() で通信を確立するだけで良い。<br /><br />
以下に簡単な例を書いておく。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">sc = SlackClient(TOKEN)
if sc.rtm_connect():
r = sc.rtm_read()
print r
sc.rtm_send_message(u"random", u"こんにちは")</span><br /><br />
こんな感じ。簡単。<br /><br />
今回のボットの場合、数字を取得してそれを素因数分解する。数字を取得すること自体は上に書いたとおり簡単なのだが、連続して数字が書かれた場合はそれを漏らさないためにキューに保存しておく必要がある。また、一定時間で解けない分解が難しい合成数などを取得した場合、タイムアウトを設定しなければならない。これらを実現するためにはスレッドを使って並列処理をする必要がある。<br /><br />
数字を取得したらキューに入れ、そのキューから別スレッドのworkerスレッドが処理を行う。workerスレッドはキューから数字を一つ取り出し、それの素因数分解を試みる。この時にタイムアウトを考慮するので、この素因数分解を行う関数についてもスレッドを利用する。そのスレッド内で時間を測り、計算が終了しないまま一定時間が過ぎたら処理を強制終了する。<br /><br />
実行方法はネットにつながっているPCから以下のコマンドを実行するだけ。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;"> $ ./prime_decomposition_bot.py</span><br /><br />
そして、実際のSlack上の画面は以下のようになる。<br /><br />
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiTAxoLHC8aq3dgMZDIa-sWK2R6vx_PCPr7FkZ5E345Nb1tIjG_CdMlv27bJCuTH0yvI76EYl82t6V13cQr__dLR0ET9fiWLjgFqtrCIv9fFAPZUmf2vkQfIW9X6oRFBAhgL8aaD11m_LL7/s1600/slack_bots_07.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiTAxoLHC8aq3dgMZDIa-sWK2R6vx_PCPr7FkZ5E345Nb1tIjG_CdMlv27bJCuTH0yvI76EYl82t6V13cQr__dLR0ET9fiWLjgFqtrCIv9fFAPZUmf2vkQfIW9X6oRFBAhgL8aaD11m_LL7/s400/slack_bots_07.png" width="400" height="228" data-original-width="935" data-original-height="534" /></a></div><br /><br />
今回書いたコードは以下の通り。今回は素因数分解をするボットだったが、これをベースに作り変えれば大抵のことはできると思う。<br /><br />
<script src="https://gist.github.com/foota/a91e53739244415466e875405dedb695.js"></script><br />
</div>
noxhttp://www.blogger.com/profile/01142433117618538300noreply@blogger.com3tag:blogger.com,1999:blog-6153301078335336992.post-71663521978410782222015-01-25T23:55:00.000+09:002015-01-25T23:59:03.025+09:00Android端末をサーバにしてPyramidフレームワークを利用するしばらく前に、PythonのWebフレームワークである<a href="http://www.pylonsproject.org/">Pyramid</a>を利用した。これがなかなか良くできており、Android端末上でも動かしてみたくなったので載せてみた。<br /><br />
ところで、自分が利用しているキャリアはドコモなんだけど、<a href="https://www.nttdocomo.co.jp/service/provider/spmode/">spモード</a>だとグローバルIPが割り振られないので外部から端末にアクセスできない。なので、spモードを契約せずに、<a href="https://www.mopera.net/">mopera U</a>を利用している。mopera UであればグローバルIPが割り振られるのでアクセスすることができるからだ。このためだけに、spモードにせず、mopera Uにしていると言っても過言ではない。<br /><br />
閑話休題。まず、Pyramidを動作させるにはAndroid端末用のPython環境である<a href="https://code.google.com/p/android-scripting/">SL4A</a>を入れる必要がある。次にPyramidを入れるのだが、必要なモジュールなどが複数あるのでそれも一緒に入れる。一応、Hello Worldプログラムを動かすのに必要なものはすべて挙げたが、もし足りない場合は実行時に何が足りないかエラー表示が出るので、それを参考に入れて欲しい。また、<a href="http://developer.android.com/intl/ja/sdk/index.html">Android SDK</a>が導入されていることを前提にしている。<br />
<span id="fullpost"></span>
<div id="fullpost"><br />
まず、Android端末にシェルで入ってプログラムを展開するディレクトリを準備する。環境に合わせてディレクトリは読み替えて欲しい。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">> adb shell
* daemon not running. starting it now on port 5037 *
* daemon started successfully *
$ cd /sdcard/sl4a/scripts
$ mkdir pyramid
$ exit</span><br /><br />
適当なディレクトリで以下のプログラム・モジュールをPC側で展開する。<br /><br />
<ul>
<li><a href="https://github.com/Pylons/pyramid">pyramid-master.zip</a></li>
<li><a href="https://pypi.python.org/pypi/translationstring">translationstring-1.3.tar.gz</a></li>
<li><a href="https://pypi.python.org/pypi/venusian">venusian-1.0.tar.gz</a></li>
<li><a href="https://pypi.python.org/pypi/WebOb/1.4">WebOb-1.4.tar.gz</a></li>
<li><a href="https://pypi.python.org/pypi/zope.deprecation">zope.deprecation-4.1.2.tar.gz</a></li>
<li><a href="https://pypi.python.org/pypi/zope.interface/4.1.2">zope.interface-4.1.2.tar.gz</a></li>
<li><a href="https://pypi.python.org/pypi/repoze.lru">repoze.lru-0.6.tar.gz</a></li>
</ul><br />
PC側からAndroid端末にプログラム・モジュールをコピーする。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">> 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/</span><br /><br />
pkg_resources.pyについては、PCに入っているPythonのsite-packagesなどから適当に持ってきてコピーする。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">> adb push pkg_resources.py /sdcard/sl4a/scripts/pyramid/</span><br /><br />
これでプログラム・モジュールについては準備完了だ。次に、Android端末側の設定をする。<br /><br />
まずはSL4Aを起動させる。pyramidディレクトリが作成されているはずなので、そこをタップして中を確認してみる。一通り、ファイルがあることを確認したら、メニューからViewをタップし、Interpretersを選択する。更にメニューを開き、Start Serverをタップして、Publicを選択する。これで、Android端末がサーバとして機能する。Androidの通知画面を開くと、SL4A Serviceという通知があるので、それをタップすることでServer情報を確認できる。そこに表示されたドメイン名でAndroid端末にアクセスできる。<a href="http://handasse.blogspot.com/2010/09/pythonandroid5.html">以前の記事</a>でもサーバにする方法を書いたので、そちらも参考にしてほしい。<br /><br />
これで、Android端末をサーバとして機能させることができたので、こんどは<a href="http://pyramid-doc-ja.readthedocs.org/en/latest/narr/firstapp.html">Hello World</a>プログラムを動かして、動作確認をしてみる。<br /><br />
プログラムをhelloworld.pyとして保存した後、Android端末にコピーする。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">> adb push helloworld.py /sdcard/sl4a/scripts/pyramid/</span><br /><br />
そして、コピー先のディレクトリにhelloworld.pyが表示されるので、それをタップして実行する。先ほど得られたドメイン名(ここでは 12345.mopera.net とする)にポート番号8080でアクセスする。<br /><br />
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjH-iOrhf0dFcbN3jnhVgY65QX6cTafeO-aRdvpMrGok_e6LwMp64B72TzKU1vbHptR3l08T5AaADyZBPDVlTTcTMdgDScF8U6C0u7Ce5xfqaJDAym0hYmK1k1OvICESxf39hyphenhyphenG9YHF7vbW/s1600/android_pyramid_20150125c.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjH-iOrhf0dFcbN3jnhVgY65QX6cTafeO-aRdvpMrGok_e6LwMp64B72TzKU1vbHptR3l08T5AaADyZBPDVlTTcTMdgDScF8U6C0u7Ce5xfqaJDAym0hYmK1k1OvICESxf39hyphenhyphenG9YHF7vbW/s400/android_pyramid_20150125c.png" /></a></div><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">http://12345.mopera.net:8080/hello/world</span><br /><br />
これで以下の画面が表示されれば成功だ。<br /><br />
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjj-6jLNIN35O4t8XSwUyccpXdYaGmAudDp4n1ti5vcOCQB16PQh-cxpZXLYWkjMS7Oh7hihN6RLAAmzw0QIin3tZoruoabUZ4jVpmq53w9gggNFQv19sjQTq3iO0Hih6zvNOUK6BwhCDIb/s1600/android_pyramid_20150125b.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjj-6jLNIN35O4t8XSwUyccpXdYaGmAudDp4n1ti5vcOCQB16PQh-cxpZXLYWkjMS7Oh7hihN6RLAAmzw0QIin3tZoruoabUZ4jVpmq53w9gggNFQv19sjQTq3iO0Hih6zvNOUK6BwhCDIb/s400/android_pyramid_20150125b.png" /></a></div><br />
携帯の電波が届くところであればどこでも使えて、Pyramidが載っている手のひらWebサーバってなんか良い。<br />
</div>noxhttp://www.blogger.com/profile/01142433117618538300noreply@blogger.com0tag:blogger.com,1999:blog-6153301078335336992.post-31184985704029510042015-01-15T20:57:00.000+09:002015-01-15T21:09:32.418+09:00IRCボットの作成久々のブログ更新。書きたいときに書くというスタンスです。最近、会社でも技術ブログを立ち上げたので、そちらで書いても良いかも。<br /><br />
今回はIRCボットの作成について書いてみる。しばらく前から社内IRCで技術関連の雑談用チャンネルを利用している。技術系のメーリングリストもあるのだけど、わざわざメールで話題にするようなことでもない、軽い内容を気軽に議論できる場が欲しかったので有志で立ち上げた。IT企業にとって技術的な雑談をできる場はとても大事だと思う。<br /><br />
素のIRCの場合、発言内容がサーバに保存されず、途中から参加しても話題についていけないという問題が出たので、適当なボットで対応することにしたが、やはりエンジニアがメインの会社なのだからボットも自作が基本でしょ、ということでPythonで作ってみた。まあ、作ったといっても「<a href="http://www.codereading.com/codereading/python/python-irc-client.html">Python でシンプルな IRC クライアントを作成する</a>」のサイトを参考にさせてもらったのでスクラッチから作成したわけではない。<br /><br />
コマンド <span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">techlog:</span> を入れると、その日の発言がボットからのtalkで返される。ただ、最初は1日分の発言しか取れない仕様だったので、日を跨ぐとその前の発言が取得できないし、一時的にチャットを抜けた分のログが欲しいだけでも強制的に1日分の発言を取得してしまう。<br /><br />
そこで、SQLite3を使ってユーザ毎にログイン・ログアウト時刻を管理して、ログアウトしてからログインするまでの発言を取得できるように変更することにした。新規に参加したメンバーでも自動的にデータベースに登録される。せっかくユーザの時刻情報があるのだから、<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">techtime:</span>, <span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">techtimeall:</span> というコマンドを作って、それぞれ自分の滞在時間、ユーザ全員分の滞在時間も取得できるようにしてみた。<br /><br />
IRCサーバ上で以下のコマンドを実行すると、IRCにボットが常駐する。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">$ ./irc_bot.py > irc_bot.log &</span><br /><br />
ソースコードは以下の通り。<br />
<span id="fullpost"></span>
<div id="fullpost">
<script src="https://gist.github.com/foota/6365c35001bc7d411fea.js"></script><br />
そこそこ社内IRCも使えるようになったのだが、最近、別プロジェクトで利用を始めた<a href="https://slack.com/">Slack</a>がかなり便利だったので、そこに技術系チャンネルを作ったところ、ほとんどのメンバーがSlackに流れてしまった…。次は<a href="https://api.slack.com/">Slackでボット作成</a>かなぁ。<br />
</div>
noxhttp://www.blogger.com/profile/01142433117618538300noreply@blogger.com0tag:blogger.com,1999:blog-6153301078335336992.post-73139744367543484202014-01-25T23:53:00.000+09:002014-01-26T01:29:46.523+09:00スピログラフ<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgzQFFPc36u-j-ed5miNaC9HkiYEnxTT8q70G2yE5_QwzbkU4vGgZn82OFy8Cr-N_oeWEPaqYBKj7cpGJyBL9XJyXHL2gqRuyBcwGfvN4lltavMVdMpJ-DtN0V07KVkNBOvx1iLhEDoTx4H/s1600/spirograph_plot_01.png" imageanchor="1" style="clear: left; float: left; margin-bottom: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgzQFFPc36u-j-ed5miNaC9HkiYEnxTT8q70G2yE5_QwzbkU4vGgZn82OFy8Cr-N_oeWEPaqYBKj7cpGJyBL9XJyXHL2gqRuyBcwGfvN4lltavMVdMpJ-DtN0V07KVkNBOvx1iLhEDoTx4H/s320/spirograph_plot_01.png" /></a></div>
スピログラフに初めて触れたのは小学校に上がる前だった。当時はそれがスピログラフという名称であることを知らず、輪っかをぐるぐる回して綺麗な模様を描ける道具としてお気に入りの玩具になっていた。自分はマンデルブロ集合やジュリア集合のような数学的に表現できる図形や模様にともて魅力を感じるのだが、幼少のころのスピログラフがその切っ掛けだったのかもしれない。<br /><br />
長方形の定規に大きさの異なる3枚の歯車が付属しており、そして2つの円がくり抜かれている。さらに、矢印、五角形、半円などの型が子供心をくすぐったものだ。因みに、歯車と円に付いている歯数は刻印されており、3つの歯車はそれぞれ36、52、63、円の方は96、105の歯を持っている。<br /><br />
幼少のころ以来、スピログラフを見かけることはなかったのだが、最近になって100円ショップでも売られているという噂を聞き、少し探ってみた。近くの100円ショップでは扱っていなかったが、「デザイン定規」という名称でネットで78円で売っていたので思わず購入してしまった(当然だが送料のほうが高くついた…)。<br /><br />
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgxb13O1C6k-N9Yr6IlFhd-JHoENelXG4CF4crkiEfBrT78Zo9zD1v2rFroFKhVJvSo2-gwg4t0qVgTiL-aeaI9L7e7kLNT0gzz09O1Igb4B6IDHfl3LLykCF4vbNX_r7TyWS2y8_CWIxms/s1600/spirograph_photo_02.jpg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgxb13O1C6k-N9Yr6IlFhd-JHoENelXG4CF4crkiEfBrT78Zo9zD1v2rFroFKhVJvSo2-gwg4t0qVgTiL-aeaI9L7e7kLNT0gzz09O1Igb4B6IDHfl3LLykCF4vbNX_r7TyWS2y8_CWIxms/s400/spirograph_photo_02.jpg" /></a></div><br />
スピログラフは一般的には内トロコイド(hypotrochoid; 内余擺線)といい、以下の式で表現できる(Wikipediaより)。<br /><br />
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj-u-8RA0Bve29QwLFZjSGfT2B2gqWeqNaBFx6tYM36taPZqyuu1sJfH1ovYC-iLLGUfFlfd6AZtsUp7pSxzF2Qqw8Tyk-Bo98BA1cuCntI8TWE9JCecXmXp6rqTu46Avmmqj0zXPx4-BmD/s1600/spirograph_formula_01.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj-u-8RA0Bve29QwLFZjSGfT2B2gqWeqNaBFx6tYM36taPZqyuu1sJfH1ovYC-iLLGUfFlfd6AZtsUp7pSxzF2Qqw8Tyk-Bo98BA1cuCntI8TWE9JCecXmXp6rqTu46Avmmqj0zXPx4-BmD/s400/spirograph_formula_01.png" /></a></div><br />
定円の半径 <i>r<sub>c</sub></i>、動円の半径 <i>r<sub>m</sub></i>、描画点の半径 <i>r<sub>d</sub></i> とし、下図では、それぞれ、青、緑、赤の線で示している。そして、回転角を θとした軌跡がスピログラフの模様となる。<br />
<span id="fullpost"><br />
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgjm215GpMu8WRO_6YUMDGLp-9T7VPktf0kpy8ehfCst1qebzHSuLqc9_3KqWKNwjRiyvzfjiXewqPbo70E55yEjJ-s-mlbmmeuo391nmGFGQ2zzOX0QYu1mCEy_P0f_6Wj0VWK2u02Nt90/s1600/spirograph_circles_01.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgjm215GpMu8WRO_6YUMDGLp-9T7VPktf0kpy8ehfCst1qebzHSuLqc9_3KqWKNwjRiyvzfjiXewqPbo70E55yEjJ-s-mlbmmeuo391nmGFGQ2zzOX0QYu1mCEy_P0f_6Wj0VWK2u02Nt90/s400/spirograph_circles_01.png" /></a></div><br />
スピログラフでは半径よりも、円周の歯の数で示した方が分かりやすいだろう。円周と半径の比は変わらないのだから、歯の数の比がそのまま半径の比となる。例えば定円の歯が30で動円の歯が10ならば、<i>r<sub>c</sub></i>:<i>r<sub>m</sub></i>=3:1ということだ。<br /><br />
さて、幼少のころはこんなに綺麗な模様が数式で表現できるなどと夢にも思わなかったわけだが、今はそれを理解し、コンピュータを使って自由に描くことができる。コンピュータであれば、動円が定円の外側に来る外トコロイド(epitrochoid; 外余擺線)も描くことができるし、描画点を動円の外に出すことも思いのままだ。<br /><br />
せっかくコンピュータを使って簡単に表示できるようになったというのに、プログラムを組むことが難しくては意味がない。しかし、Python+matplotlibであれば簡単に描くことができる。上述の式をそのままコードに落とすだけだ。以下がIPythonで描いたスピログラフの模様だ。図と合わせて一画面に収まってしまう。<br /><br />
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh6jynJ0iBgM6_Ebc1n6FBCswhP_tVMemET4gnWLUd9c9vKJ4Sp_dcjd7g1GDfrb1VSFUSjJOmtoaUIfnqFMcwReZh6GghjWeDHOphv_vNGB0LalJenyWey6VrkKiQpH6T5oVJLO5sG00lo/s1600/spirograph_ipython_01.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh6jynJ0iBgM6_Ebc1n6FBCswhP_tVMemET4gnWLUd9c9vKJ4Sp_dcjd7g1GDfrb1VSFUSjJOmtoaUIfnqFMcwReZh6GghjWeDHOphv_vNGB0LalJenyWey6VrkKiQpH6T5oVJLO5sG00lo/s400/spirograph_ipython_01.png" /></a></div><br />
ここでは、Pythonのソースコードも挙げておこう。せっかくなので少しだけmatplotlibの調整をしてみた。必ず縦と横のスケールが同じになるようにし、グラフ表示も正方形にした。コードを見てもらえばやり方はすぐに解るだろう。以下のコマンドで実行できる。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">$ spirograph.py <i>r<sub>c</sub></i> <i>r<sub>m</sub></i> <i>r<sub>d</sub></i></span><br /><br />
例えば次のように実行すれば冒頭の図が表示される。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">$ spirograph.py 96 63 30</span><br /><br />
<script src="https://gist.github.com/foota/8616573.js"></script><br />
最後に参考にしたサイトを挙げておく。<br /><br />
<ul>
<li><a href="http://ja.wikipedia.org/wiki/%E3%82%B9%E3%83%94%E3%83%AD%E3%82%B0%E3%83%A9%E3%83%95">スピログラフ (Wikipedia)</a></li>
<li><a href="http://ja.wikipedia.org/wiki/%E3%83%88%E3%83%AD%E3%82%B3%E3%82%A4%E3%83%89">トロコイド (Wikipedia)</a></li>
<li><a href="http://stillbe.web.fc2.com/compendium/basic/spiro/body03.html">スピログラフ 第3回 トロコイド (彷徨の神殿)</a></li>
</ul><br />
</span>noxhttp://www.blogger.com/profile/01142433117618538300noreply@blogger.com1tag:blogger.com,1999:blog-6153301078335336992.post-42451579190220295872014-01-01T16:46:00.000+09:002014-01-01T16:46:06.130+09:002013年と2014年の違い2013年を終え、2014年を迎えた。さて、昨年と今年の違いは何か?<br /><br />
昨年のブログで2013を素因数分解すると「3×11×61」になることを示したが、2014は「2×19×53」となる。そして差を取ると、「2-3, 19-11, 53-61」、つまり、「-1, 8, -8」となり、1月、8月が昨年と今年の差であることがわかる。そして、負を示した1月は運気が下がり、正と負を持つ8月は運気の上下があるに違いない! 因みに2013+2014=4027が素数にあることから、昨年から今年にかけて行っていることは、これまでとは違ったユニークな内容となる!<br /><br />
…などと、MMRばりなことを書けるぐらいの余裕を持ちつつ、今年一年を頑張りたいと思いますので、本年もどうぞよろしくお願いいたします。<br />
<span id="fullpost"><br />
以下にPython+<a href="http://tnt.math.se.tmu.ac.jp/nzmath/">NZMATH</a>を使った素因数分解を示す。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">>>> 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)]</span><br />
</span>
noxhttp://www.blogger.com/profile/01142433117618538300noreply@blogger.com1tag:blogger.com,1999:blog-6153301078335336992.post-44382330004109414462013-09-21T21:40:00.000+09:002013-09-22T01:13:13.373+09:00クッキークリッカーとプログラマ<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgpXjhoIa9bL3qMoiPoi0TUfzSWTXTNvjlje2fit5Ve8zu3h3XE59oVQ2cQ6LJWz0edlaBtt0gqp5ozBxfq4jpeFQWyHw21cR8soCU5FOJ7Gz2UfINi2Mcs-wM_PKBwTWtpMfooCIHrsTED/s1600/cookie_clicker_01.png" imageanchor="1" style="clear: left; float: left; margin-bottom: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgpXjhoIa9bL3qMoiPoi0TUfzSWTXTNvjlje2fit5Ve8zu3h3XE59oVQ2cQ6LJWz0edlaBtt0gqp5ozBxfq4jpeFQWyHw21cR8soCU5FOJ7Gz2UfINi2Mcs-wM_PKBwTWtpMfooCIHrsTED/s320/cookie_clicker_01.png" /></a></div>
先週ぐらいから<a href="http://orteil.dashnet.org/cookieclicker/">クッキークリッカー(Cookie Clicker)</a>というJavaScriptを使ったブラウザゲームが流行っている。クリックするだけのゲームと聞いて、最初はあまり興味を持てなかったのだが、自分の周りであまりにもやっている人が多いので少し遊んでみることにした。<br /><br />
クッキークリッカーを簡単に説明すると、まずはクリックすることでクッキーを作り、作ったクッキーを使ってクッキーの生産性を高めるためのアップグレードやアイテムを購入し、たまに出現するゴールデンクッキー(Golden Cookie)をクリックすることでさらに多量のクッキーが得られるので、それらを駆使してできるだけたくさんのクッキーを作るというゲームだ。<br /><br />
このようにとてもシンプルなゲームなのだが、最初はちまちまとしか作れなかったクッキーが徐々に増えていき、様々なアイテムやイベントを通すことで、終盤では毎秒数千億クッキーを作れるというところに人々を惹きつける魅力があるらしい。しかし、自分は「人間の代わりにプログラムを働かせることでどこまで効率が良くなるか」という目的のために自動化プログラムを作成したくなった。コンピュータは人間を楽にさせるために存在するべきだ。<br />
<span id="fullpost"><br />
このゲームはJavaScriptで作られており、そのコードや変数を修正すれば簡単にチートを行うことができる。しかし、自分はそれについては全く興味がない。元のシステムには全く手を入れずに人間が操作する部分のみをプログラムで最適化したいのだ。<br /><br />
まず、最初に考えついたのは自動クリックだ。これは様々なツールが世に溢れているのでそれを使っても良いのだが、プログラマだったら、このぐらいは自分で作りたい。そこで、Pythonと<a href="https://github.com/SavinaRoja/PyUserInput">pymouse</a>を使って自動クリックプログラムを作成した。これで、毎秒100~200クリックすることができるようになった。これ以上の速度にするとブラウザがついて行けずプチフリーズを頻繁に起こし逆に生産性が落ちてしまう。因みに画面上のエフェクトなどは設定で消しておくほうが良いだろう。<br /><br />
次に、ゴールデンクッキーへの対処だ。ゴールデンクッキーはランダムな間隔でブラウザ上のどこかに出現する。そのクッキーをクリックすると溜まっているクッキーの1割分の枚数が貰えたり、一定時間現在の7倍の効率でクッキーを焼けるようになったりする。この効果はかなり大きく是非ともこれを自動化したい。<br /><br />
最初は一般的なテンプレートマッチングにより検出しようと思ったが、ゴールデンクッキーは360度ランダムに回転した状態で出現するので、それに対処しようとすると、<a href="http://www.pythonware.com/products/pil/">PIL</a>で簡単に済ませるには処理時間が掛かってしまいそうだし、<a href="http://opencv.org/">OpenCV</a>を入れるのも面倒だし、クッキー1枚ごときの検出で手間を掛けたくない。そこで、ゴールデンクッキーだけに使われている色を検出することにした。まず、PILでデスクトップのキャプチャを行い、ブラウザの領域のみをクロップし、その領域内でゴールデンクッキーで使われている色を検出できれば、その位置がクリックする位置になる。検出されなければクリックを行わない。<br /><br />
以上の自動化でかなり効率よく稼げるようになり、最終的に作ったクッキーの約7割がクリックによる生産になった。1回目は10京個以上のクッキーを作ったあとにResetを行い、2回目でアップグレードが100%になるまでに要した時間は「リサーチ」の待ち時間を含めて9時間程だった。途中で一時的に中断したり、効率的にアップグレードしたわけではないので、その辺を考えてクッキーを焼けばもっと早く100%になったと思う。また、9時間と言っても自動化プログラムが動いている間は席を離れていても良いので、自分で操作したのはアイテム購入ぐらいだったが。<br /><br />
というわけで、自分はクッキークリッカーをそれなりに楽しんだ。楽しんだ部分は主に自動化プログラムを作成したところなので、一般にプレイされている方とは楽しみ方が違うのかもしれないけど。最後に作成したプログラムを以下に示しておく。<br /><br />
<script src="https://gist.github.com/foota/6648188.js"></script><br />
</span>
noxhttp://www.blogger.com/profile/01142433117618538300noreply@blogger.com2tag:blogger.com,1999:blog-6153301078335336992.post-41568563109423229292013-05-26T21:45:00.000+09:002013-05-30T01:53:48.329+09:00Pythonを使って簡単にデータを視覚化する<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiBe7QQmZO34-g-TF9OCcwpj0P93JZUy2LoBuPz1UkVyOzRmNDtiTdVcAfcn22e5ovlXhqQHkiW-gtf2VDb6E1c5BvRhd0ojrXQEplZrybHzC2VwP9Sn6A9HR8lmHHfMXdn-z2YIV4EuRTf/s1600/ipython_spectral_clustering_05.png" imageanchor="1" style="float:left;margin:0 10px 10px 0;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiBe7QQmZO34-g-TF9OCcwpj0P93JZUy2LoBuPz1UkVyOzRmNDtiTdVcAfcn22e5ovlXhqQHkiW-gtf2VDb6E1c5BvRhd0ojrXQEplZrybHzC2VwP9Sn6A9HR8lmHHfMXdn-z2YIV4EuRTf/s200/ipython_spectral_clustering_05.png" /></a>
世の中のことをもっと知るにはどうしたら良いだろうと思うときがある。世の中の多くの事柄はログやデータに落とされる。Googleなどの検索サイトは良い例だろう。さて、そのログやデータをどうすれば良いのか? 多くの場合、視覚化が有効な手段となる。<br /><br />
まずは身の回りの日常的なデータやログを何とかしたい。ただ、日常のデータを視覚化するのに数十行以上のコードは書きたくない。まるで息をするかのごとく自然に視覚化を行いたいのだ。そのためには1~2行、長くて数行で済ませることが必要だ。そこでPythonとmatplotlibを使う。加えて、IPythonがあればなお良い。IPythonの導入については以前のブログ記事である<a href="http://handasse.blogspot.com/2012/03/ipython.html">IPythonの埋め込みプロットが素晴らしい</a>を参考にして欲しい。<br />
<span id="fullpost"><br />
まずは事前にnumpyとmatplotlibをインポートしておく。できればscipyも。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">>>> from numpy import *
>>> from pylab import *</span><br /><br />
短いコードで視覚化を行うためには、Pythonの内包表記は必須だ。例えば、5, 2, 1, 5, 8をデータとするグラフを書きたいのならIPythonを使って以下の1行で実現できる。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">>>> plot([5, 2, 1, 5, 8])</span><br /><br />
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjNUnuXPAsFegNiM7cTfxVpTWgWf8nxYpYxG3prUU0PURETq3ERNmCcvTvMYjeGv2DN7zQOx6jfcnHUgWkutxAM8_EE4MDapHLIXmenVzJdql3PoKawHdM-HvPgr8uKXHEf19B16YB0LS5m/s1600/ipython_plot_01.png" imageanchor="1" ><img border="0" style="text-align:center;" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjNUnuXPAsFegNiM7cTfxVpTWgWf8nxYpYxG3prUU0PURETq3ERNmCcvTvMYjeGv2DN7zQOx6jfcnHUgWkutxAM8_EE4MDapHLIXmenVzJdql3PoKawHdM-HvPgr8uKXHEf19B16YB0LS5m/s320/ipython_plot_01.png" /></a><br /><br />
数値が行ごとにinput.txtというファイルに書かれていた場合は以下の通り。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">>>> plot([int(x) for x in file("input.txt")])</span><br /><br />
map関数を使えばもっとスマートに書ける。ファイル内の数値を文字列として配列で読み込み、それらの文字列をmap関数によりintで整数に変換する。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">>>> plot(map(int, file("input.txt")))</span><br /><br />
さて、ログが2次元で書かれていた場合はどうするのか。例えば以下のデータがinput.txtに書かれていたとする。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">2, 2.5
5, 6.2
6, 3.6
7, 6.3
10, 1.9</span><br /><br />
この場合、以下のようにすればプロットできる。ファイルから行ごとの文字列を読み出して、それを","を区切りとしてリスト化する。文字列で格納されている数値をmap関数を使ってfloatで実数に変換している。最後にmap(list, zip(*data))を使って転置している。ここでdataはリストによる行列とする。リストによる行列の転置はイディオムとして覚えておくと便利だ。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">>>> x, y = map(list, zip(*[map(float, xs.split(",")) for xs in file("input.txt")]))
>>> plot(x, y)</span><br /><br />
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhF1QZkpo_QM-uZ_jJQZQ1IzQe2eADKz4qi4ac9th1Ay4hlhPPClHKXvS7X4hZQOUj87uLSf8FsVd44I6QLyXMW7AyPNTkTqjXU4L_1nS8sd0ZmybecyP7mhX2haTzTNfZbj0ET1iClT5dx/s1600/ipython_plot_02.png" imageanchor="1" ><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhF1QZkpo_QM-uZ_jJQZQ1IzQe2eADKz4qi4ac9th1Ay4hlhPPClHKXvS7X4hZQOUj87uLSf8FsVd44I6QLyXMW7AyPNTkTqjXU4L_1nS8sd0ZmybecyP7mhX2haTzTNfZbj0ET1iClT5dx/s320/ipython_plot_02.png" /></a><br /><br />
ログファイルなどの場合、単に数値が並んでいることは少なく一緒にテキストが書き込まれていることが多い。例えば次のログ(input.log)について視覚化を行なってみる。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">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</span><br /><br />
Unitの値をそのままプロットするなら以下のようにすれば良い。予めreモジュールをインポートしておき、re.matchで先頭の文字列と一致したものだけ選び出している。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">>>> plot([float(l.split("=")[1]) for l in file("input.log") if re.match("Unit =", l)])</span><br /><br />
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgRUH-7gJMTQG-_p5YOc2WKa7b_y1jD9aa_338lnF3ASx_kET_SCxq18E4JosrJ00cUDzAd9iXc_1iymgrTg-WZkEgFtjGev16vhL89MnKlqpfLl3KutUl3kkmDDgliWVscqOZrBLvVm1cm/s1600/ipython_plot_log_01.png" imageanchor="1" ><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgRUH-7gJMTQG-_p5YOc2WKa7b_y1jD9aa_338lnF3ASx_kET_SCxq18E4JosrJ00cUDzAd9iXc_1iymgrTg-WZkEgFtjGev16vhL89MnKlqpfLl3KutUl3kkmDDgliWVscqOZrBLvVm1cm/s320/ipython_plot_log_01.png" /></a><br /><br />
Unitを横軸、Factorを縦軸にしてドットでプロットしたい場合は以下のようにする。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">>>> 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")</span><br /><br />
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg89sRS_7aTmwHwHwkg7a5GPtS0WVc42NkfAgVcfFmInmppnDml0vVDaakcfbFmqW2c-aV2jz2JdgNh0u8jFUsOVPaYb74ktFEx6RPNoGQFL2voPQw6xBM_v8chvSXuR58j7_v-caDNfAlZ/s1600/ipython_plot_log_02.png" imageanchor="1" ><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg89sRS_7aTmwHwHwkg7a5GPtS0WVc42NkfAgVcfFmInmppnDml0vVDaakcfbFmqW2c-aV2jz2JdgNh0u8jFUsOVPaYb74ktFEx6RPNoGQFL2voPQw6xBM_v8chvSXuR58j7_v-caDNfAlZ/s320/ipython_plot_log_02.png" /></a><br /><br />
最後に少し複雑な例を挙げておく。あるデータに対して<a href="http://en.wikipedia.org/wiki/Spectral_clustering">スペクトラルクラスタリング</a>を行いたいとする。まずは、numpyとmatplotlibの他に、scipy関連の固有値問題・k近傍法・kd木のモジュールをインポートしておく。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">>>> from scipy.linalg import eig
>>> from scipy.cluster.vq import kmeans2
>>> from scipy.spatial.kdtree import KDTree</span><br /><br />
そして、クラスタリングしたいデータを<a href="https://gist.github.com/foota/91b8939413baa3f44a28">circles.dat</a>とすると、以下のコードで視覚化できる。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">>>> 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)]</span><br /><br />
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiBe7QQmZO34-g-TF9OCcwpj0P93JZUy2LoBuPz1UkVyOzRmNDtiTdVcAfcn22e5ovlXhqQHkiW-gtf2VDb6E1c5BvRhd0ojrXQEplZrybHzC2VwP9Sn6A9HR8lmHHfMXdn-z2YIV4EuRTf/s1600/ipython_spectral_clustering_05.png" imageanchor="1" ><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiBe7QQmZO34-g-TF9OCcwpj0P93JZUy2LoBuPz1UkVyOzRmNDtiTdVcAfcn22e5ovlXhqQHkiW-gtf2VDb6E1c5BvRhd0ojrXQEplZrybHzC2VwP9Sn6A9HR8lmHHfMXdn-z2YIV4EuRTf/s320/ipython_spectral_clustering_05.png" /></a><br /><br />
ここではforループを行うために内包表記を使っており、また、内包表記のリストが必要ない場合は、_(アンダースコア)に入れておく。これはシェル上で余計なデータを表示させないためである。また、内包表記上では配列への代入ができないので__setitem__を利用している。さらに、Pythonのリストの代わりにnumpyのarrayを利用すれば転置などもmap(list, zip(*data))の代わりにarray(data).Tとすれば良いので分かりやすい。<br /><br />
Pythonを使えばたったこれだけだ。これだけのことで、これまで見えなかったものが見えてくるのは楽しい。<br />
</span>noxhttp://www.blogger.com/profile/01142433117618538300noreply@blogger.com0tag:blogger.com,1999:blog-6153301078335336992.post-17621296460469154162013-03-10T19:50:00.000+09:002013-03-10T19:53:15.765+09:00Python: timeitモジュールを使ってお手軽に実行時間計測Pythonでは<a href="http://docs.python.jp/2/library/timeit.html">timeitモジュール</a>を使って簡単に実行時間を計測できる。以下に示すコードはPythonでは一般的なコードだと思う。<a href="http://ja.wikipedia.org/wiki/%E7%AB%B9%E5%86%85%E9%96%A2%E6%95%B0">たらい回し関数(竹内関数)</a>と階乗を求める関数だ。これを実行した時の時間計測、さらにコードに含まれている関数の実行時間を簡単に測るにはどうすればよいか。<br /><br />
<script src="https://gist.github.com/foota/5127840.js"></script><br />
ここでtimeitモジュールが利用できる。まず、Pythonシェル環境では以下の通り。<br /><br />
main関数の測定は次のように行う。stmtで時間計測したい関数を指定して、setupで最初の1回だけ実行する文を指定する。numberは実行する回数となっている。正確を期するなら複数回を指定するべきだろう。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">>>> import timeit
>>> timeit.timeit(stmt='test_timeit.main()', setup='import test_timeit', number=1)
12
3628800
2.663659573618423</span><br /><br />
たらい回し関数(tarai)のみの測定。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">>>> timeit.timeit(stmt='test_timeit.tarai(12, 6, 0)', setup='import test_timeit', number=1)
2.6390463109480637</span><br /><br />
階乗関数(factorial)のみの測定。100万回実行。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">>>> timeit.timeit(stmt='test_timeit.factorial(10)', setup='import test_timeit', number=1000000)
2.3606330638673825</span><br /><br />
さらに、コマンドラインでの計測も可能だ。<br />
<span id="fullpost"><br />
main関数の測定。コマンドオプションの-sは最初に1回だけ実行する文、-nは実行する文の回数、-rはタイマのリピート数になっている。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">$ 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</span><br /><br />
たらい回し関数(tarai)のみの測定。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">$ 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</span><br /><br />
階乗関数(factorial)のみの測定。100万回実行。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">$ 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</span><br /><br />
ところで、上に示した階乗関数については「車輪の再発明」なので、通常はmath.factorialを使った方が良いだろう。<br />
</span>noxhttp://www.blogger.com/profile/01142433117618538300noreply@blogger.com0tag:blogger.com,1999:blog-6153301078335336992.post-89264164258014634762013-01-01T20:00:00.000+09:002013-01-01T20:00:30.959+09:002013年と素因数分解あけましておめでとうございます。<br /><br />
新年を迎えて、ふと、今年の西暦である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 となる。<br /><br />
大きな素数同士による合成数の素因数分解は非常に困難であることが知られている。この性質を利用して多くの暗号アルゴリズムで大きな素数による合成数が利用されている。現在のところ、素因数分解アルゴリズムとしては、楕円曲線法(ECM; Elliptic Curve Method)や複数多項式二次ふるい法(MPQS; Multiple Polynomial Quadratic Sieve)がよく使われているらしい。で、手軽にこれらを利用できるライブラリがないかと軽く調べたところ<a href="http://tnt.math.se.tmu.ac.jp/nzmath/index_j.html">NZMATH(ニジマス)</a>という数論のためのPythonモジュールがあることを知ったので、早速インストールして使ってみた。<br />
<span id="fullpost"><br />
ところで、大きな素数といえばメルセンヌ素数が有名だが、<a href="http://ja.wikipedia.org/wiki/%E3%83%9E%E3%83%A9%E3%83%B3%E3%83%BB%E3%83%A1%E3%83%AB%E3%82%BB%E3%83%B3%E3%83%8C">マラン・メルセンヌ</a>は、2<sup>n</sup>-1 に対してnが257以下のとき、n = 2, 3, 5, 7, 13, 17, 19, 31, 67, 127, 257のときにのみ素数となると主張した。しかし、一部に間違いがあり、n = 61, 89, 107のときも素数であり、また、n = 67, 257は素数ではなく合成数であった。<br /><br />
そこで、誤って素数としていたn = 67のときの 147573952589676412927 について楕円曲線法(ECM)で解かせると、すぐに 193707721×761838257287 と結果が出力された。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">>>> import nzmath.factor.methods as methods
>>> methods.factor(147573952589676412927, method="ecm")
[(193707721L, 1), (761838257287L, 1)]</span><br /><br />
更に、メルセンヌが見落としていた9番目と10番目のメルセンヌ素数である 2305843009213693951 (2<sup>61</sup>-1) と 618970019642690137449562111 (2<sup>89</sup>-1) を掛け合わせた46桁の合成数 1427247692705959880439315947500961989719490561 を複数多項式二次ふるい法(MPQS)で解いてみたところ、30分ほどで正しく2つの素数に分解できた。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">>>> methods.factor(1427247692705959880439315947500961989719490561, method="mpqs")
[(2305843009213693951L, 1), (618970019642690137449562111L, 1)]</span><br /><br />
ソースコードも簡単に確認できるし、NZMATHはなかなか楽しい。<br /><br />
それでは、本年もよろしくお願い申し上げます。<br />
</span>noxhttp://www.blogger.com/profile/01142433117618538300noreply@blogger.com0tag:blogger.com,1999:blog-6153301078335336992.post-37680787513747557812012-09-15T23:09:00.000+09:002012-09-24T00:44:00.154+09:00「フカシギの数え方」の問題を解いてみた先日、「<a href="http://youtu.be/Q4gTV4r0zRs">『フカシギの数え方』 おねえさんといっしょ! みんなで数えてみよう!</a>」という動画を見た。格子状のマスの左上から右下までの経路が何通りあるのかを調べて、格子が多くなればなるほど組み合わせの数が爆発的に増えることを教えてくれる動画だ。これは<a href="http://en.wikipedia.org/wiki/Self-avoiding_walk">自己回避歩行(Self-avoiding walk)</a>と呼ばれている問題らしい。<br /><br />
これだけ聞いてもそれほどインパクトはないのだが、動画に出てくるおねえさんの経路を調べあげる執念がもの凄く、ネット上でも結構な話題になっている。執念と言うよりも狂気に近い。しかし、話題になった割には動画内で言及されている高速なアルゴリズムを実装したという話を聞かなかったので、自分で確かめることにした。<br /><br />
<div align="center"><iframe width="560" height="315" src="http://www.youtube.com/embed/Q4gTV4r0zRs" frameborder="0" allowfullscreen></iframe></div><br />
<iframe src="http://rcm-jp.amazon.co.jp/e/cm?t=footawwwservi-22&o=9&p=8&l=as1&asins=4048687409&nou=1&ref=qf_sp_asin_til&fc1=000000&IS2=1<1=_blank&m=amazon&lc1=0000FF&bc1=FFFFFF&bg1=FFFFFF&f=ifr" style="width:120px;height:240px;float:left;margin:0 10px 10px 0;" scrolling="no" marginwidth="0" marginheight="0" frameborder="0"></iframe>
動画のおねえさんは深さ優先探索によるプログラムを使っていると思われるが、それだとスパコンを使っても10×10マスの格子を解くのに25万年も掛かってしまう。そこで、高速化のためにゼロサプレス型二分決定グラフ(ZDD; Zero-Suppressed Binary Decision Diagram)と呼ばれるアルゴリズムを利用することにした。このアルゴリズムを開発したのは北大の<a href="http://www-alg.ist.hokudai.ac.jp/~minato/index-j.html">湊先生</a>で、ZDDによりすべての経路を見つけ出すアルゴリズムとして<a href="http://www-cs-faculty.stanford.edu/~uno/">クヌース先生</a>のSIMPATHを使った。ZDDについてはクヌース先生も強い関心を持っていて、
<a href="http://www.amazon.co.jp/gp/product/4048687409/ref=as_li_qf_sp_asin_tl?ie=UTF8&camp=247&creative=1211&creativeASIN=4048687409&linkCode=as2&tag=footawwwservi-22">The Art of Computer Programming Volume 4, Fascicle 1</a><img src="http://www.assoc-amazon.jp/e/ir?t=footawwwservi-22&l=as2&o=9&a=4048687409" width="1" height="1" border="0" alt="" style="border:none !important; margin:0px !important;" />(TAOCP 4-1)ではBDD/ZDDの詳細な解説を読むことができる。演習問題の解説だけで書籍の半分を使っていることからしても気合の入れようがわかるだろう。<br /><br />
実際に自分のノートPCでZDDアルゴリズムを使ったコードを走らせたら、ほんの10秒程度で10×10マスの問題を解いてしまった。おねえさんがスパコンで25万年かかった問題をノートPCでたった10秒である。約8千億倍の高速化だ。これだけ劇的に変わるとやっぱり楽しい。そして、アルゴリズムの重要性を再認識させられた。<br />
<span id="fullpost"><br />
さて、以下におねえさんが利用したであろう自作の深さ優先探索(DFS)プログラムと高速な解法であるZDDアルゴリズムの両方を載せておく。ZDDについては4つのプログラムを1つのPythonスクリプトで統合している。これは、クヌース先生の<a href="http://www-cs-faculty.stanford.edu/~uno/programs/simpath.w">SIMPATH</a>および<a href="http://www-cs-faculty.stanford.edu/~uno/programs/simpath-reduce.w">SIMPATH-REDUCE</a>を利用しており、また、<a href="http://www-cs-faculty.stanford.edu/~uno/sgb.html">SGB (Stanford GraphBase)ライブラリ</a>でグラフを作成し、経路の数え上げに<a href="http://gmplib.org/">GMP</a>を使ったC++で自作コードを組んだためである。最初はSIMPATHを参考に実装して一つのソースコードに統一しようかと思ったが、解説もちゃんと書かれているクヌース先生のコードをそのまま利用することにした。これらのプログラムの簡単な解説と使い方を載せておく。<br /><br />
まずはDFSによる実装を以下に示す。Node構造体から経路を辿って数え上げるだけの単純なプログラムである。<br /><br />
<script src="https://gist.github.com/3721400.js?file=count_lattice_path_dfs.cpp"></script><br />
実行方法は以下の通り。引数の7は7×7のノードを表しており、マスで表現すれば6×6マスとなる。そして計算結果として、575,780,564通りの経路が探索されたことを示している。因みに7×7ノード(6×6マス)で約3分ほどの時間が掛かっており、それ以上の大きさだと現実的な時間で解くことが難しくなってくる。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">% ./count_lattice_path_dfs 7
Count all paths from (0, 0) to (6, 6) in the 7x7 lattice graph:
Count: 575780564</span><br /><br />
次に、ZDDを利用したプログラムのビルド方法を示す。<br /><br />
ビルドをする前に<a href="http://www-cs-faculty.stanford.edu/~uno/sgb.html">SGB</a>および<a href="http://gmplib.org/">GMP</a>を準備しておく。また、<a href="http://www-cs-faculty.stanford.edu/~uno/programs/simpath.w">SIMPATH</a>および<a href="http://www-cs-faculty.stanford.edu/~uno/programs/simpath-reduce.w">SIMPATH-REDUCE</a>については<a href="http://www-cs-faculty.stanford.edu/~uno/cweb.html">CWEB</a>を利用しているので、それも入れておく。<br /><br />
最初に、SGBライブラリを用いて格子状のグラフを作成する。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">% gcc gen_lattice.c -O3 -lgb -o gen_lattice</span><br /><br />
ソースコード内でグラフを下記のように作成している。board関数はチェスのような格子状のグラフを作成することができる。他にもグラフを作成する便利な関数が多数あるので、興味のある方はSGBライブラリのgb_basic.wあたりを参照して欲しい。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">Graph* g = board(N, M, 0L, 0L, 1L, 0L, 0L);</span><br /><br />
<script src="https://gist.github.com/3721375.js?file=gen_lattice.c"></script><br />
次に、SIMPATHおよびSIMPATH-REDUCEを以下のようにビルドする。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">% wget http://www-cs-faculty.stanford.edu/~uno/programs/simpath.w
% ctangle simpath.w
% gcc simpath.c -O3 -lgb -o simpath</span><br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">% 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</span><br /><br />
因みに、ctangleでC言語への変換、cweaveでTeX形式への変換を行うことができるので、以下のようにTeX形式のドキュメントも一緒に作っておくことをお勧めする。また、これらのアルゴリズムを理解するために前述のTAOCP Vol. 4を読んでおくことが望ましい。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">% cweave simpath.w
% tex simpath.tex
% dvipdf simpath.dvi # PDFファイルが必要なら.
% cweave simpath-reduce.w
% tex simpath-reduce.tex
% dvipdf simpath-reduce.dvi # PDFファイルが必要なら.</span><br /><br />
最後に、ZDDから経路を数え上げるプログラムをコンパイルする。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">% g++ count_path_mp.cpp -O3 -lgmp -o count_path_mp</span><br /><br />
経路を一つ一つ数えていては高速に数え上げるという点で意味がなくなってしまうので、ちゃんとメモ化再帰で数える。ZDDは高度に圧縮されているのでメモ化再帰でほとんど瞬時に数え上げることができる。<br /><br />
<script src="https://gist.github.com/3721384.js?file=count_path_mp.cpp"></script><br />
そして、上記のプログラムを統括するPythonスクリプトが以下となる。グルー言語としてもPythonは優秀だ。<br /><br />
<script src="https://gist.github.com/3721389.js?file=count_lattice_path.py"></script><br />
実行方法は以下の通り。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">% ./count_lattice_path.py 8
Count all paths from (0, 0) to (7, 7) in the 8x8 lattice graph:
Count: 789360053252</span><br /><br />
DFSのプログラム同様、引数の8は8×8ノード(7×7マス)を表している。そして計算結果として、789,360,053,252通りの経路が探索されたことを示している。また、中間ファイルとして以下が出力される。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">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 # 一連の実行ファイルのログ.</span><br /><br />
因みに、DFSで7×7ノード(6×6マス)を解くと3分ほどかかったが、ZDDでは0.1秒未満で完了した。また、ZDDで14×14ノード(13×13マス)は数分で終わることを確認したが、それ以上は途中で搭載メモリ以上に必要なメモリが大きくなるので実行していない。<br />
</span>
noxhttp://www.blogger.com/profile/01142433117618538300noreply@blogger.com0tag:blogger.com,1999:blog-6153301078335336992.post-5240296150189604462012-05-14T01:58:00.000+09:002012-05-14T11:05:00.693+09:00テトレーション・フラクタル<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhg6dieaZBdkwvOn5sPkgEG75cZRzsgipet65bblgHfV8EseWza6uJSGha0aD-l3aeUOJCQRwh41VaDJhF8vruDGMTncRV3DJNJ09DV15SS21e_P5W15yhMTZtXyxpRnIYaaKNcnHoKpSfa/s1600/tetration.png" imageanchor="1" style="clear:left; float:left;margin-right:1em; margin-bottom:1em"><img border="0" height="320" width="320" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhg6dieaZBdkwvOn5sPkgEG75cZRzsgipet65bblgHfV8EseWza6uJSGha0aD-l3aeUOJCQRwh41VaDJhF8vruDGMTncRV3DJNJ09DV15SS21e_P5W15yhMTZtXyxpRnIYaaKNcnHoKpSfa/s320/tetration.png" /></a></div>
5年ほど前、このブログで<a href="http://handasse.blogspot.com/2007/07/blog-post_8234.html">虚数のテトレーション</a>という記事を書いたことがある。<a href="http://en.wikipedia.org/wiki/Tetration">テトレーション(tetration)</a>とは、自らのべき乗を指定された回数反復する演算のことで、<sup><i>n</i></sup><i>a</i> と表現する。<sup>3</sup>5 の場合、5<sup>5<sup>5</sup></sup> = 1.911×10<sup>2185</sup> となる。Pythonの関数で表現すれば以下のようになる。<br /><br />
<script src="https://gist.github.com/2689028.js"> </script><br />
以前の記事では、<sup>∞</sup><i>i</i> が 0.43828+0.36059<i>i</i> に収束することを見つけたのだが、今回はそれに関連したフラクタルについて紹介したい。<br /><br />
テトレーション <sup><i>n</i></sup><i>a</i> の <i>a</i> には複素数を指定することができる。このとき、<i>a = x + yi</i> として、それを複素平面に置く。ここで正の整数<i>n</i>を大きくしていき、発散と判定された<i>n</i>に対応した色で平面を色分けする。発散しない場合、予め決めておいた<i>n</i>の上限値を使う。これで得られる図を<a href="http://www.tetration.org/">テトレーション・フラクタル(tetration fractal)</a>と呼ぶ。<br /><br />
<sup><i>n</i></sup><i>(x + yi)</i>=<i>a + bi</i>としたとき、<sup><i>n</i>+1</sup><i>(x + yi)</i>=<i>a' + b'i</i> は以下のように計算できる。<br /><br />
<div align="center">
<a href="http://formula.s21g.com/?a%5E%5Cprime%20%3D%20e%5E%7B-%5Cfrac%7B1%7D%7B2%7D%5Cpi%20b%7D%5Ccos%5Cfrac%7B%5Cpi%20a%7D%7B2%7D"><img src="http://formula.s21g.com/?a%5E%5Cprime%20%3D%20e%5E%7B-%5Cfrac%7B1%7D%7B2%7D%5Cpi%20b%7D%5Ccos%5Cfrac%7B%5Cpi%20a%7D%7B2%7D.png" /></a><br />
<a href="http://formula.s21g.com/?b%5E%5Cprime%20%3D%20e%5E%7B-%5Cfrac%7B1%7D%7B2%7D%5Cpi%20b%7D%5Csin%5Cfrac%7B%5Cpi%20a%7D%7B2%7D"><img src="http://formula.s21g.com/?b%5E%5Cprime%20%3D%20e%5E%7B-%5Cfrac%7B1%7D%7B2%7D%5Cpi%20b%7D%5Csin%5Cfrac%7B%5Cpi%20a%7D%7B2%7D.png" /></a><br />
</div>
<span id="fullpost"><br />
Pythonでの実装は<a href="http://code.activestate.com/">ActiveSate Code</a>に<a href="http://code.activestate.com/recipes/577917-tetration-fractal/">Tetration Fractal</a>として載っている。今回はこれを少し修正したコード、C++で書いたコード、更にTBBで並列化したコードを用意して、それらの実行速度のベンチマークを取ってみた。<br /><br />
<table border="1" cellspacing="0" cellpadding="2" align="center">
<tr align="center">
<td width="100">Python</td>
<td width="100">C++</td>
<td width="100">C++ with TBB</td>
</tr>
<tr align="center">
<td>258.732秒</td>
<td>16.999秒</td>
<td>0.976秒</td>
</tr>
</table><br />
ここで、複素平面の領域は(-1.5, 0)-(-0.75, 0.75)であり、最大繰り返し数は256としている。画像の大きさは1,024×1,024とした。実行環境はIntel Xeon X5650の2CPUで、12コア・24スレッドとなっている。<br /><br />
以下に、ソースコードを示す。詳細についてはコードを読んでいただきたい。<br /><br />
Pythonによる実装:<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">% ./tetration.py tetration.png -1.5 0.0 0.75 0.75 1024 1024</span><br /><br />
<script src="https://gist.github.com/2689040.js?file=tetration.py"></script><br /><br />
C++による実装:<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">% g++ tetration.cpp -o tetration -O3</span><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">% ./tetration tetration.out -1.5 0.0 0.75 0.75 1024 1024</span><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">% ./image_tetration.py tetration.out tetration.png</span><br /><br />
<script src="https://gist.github.com/2689074.js?file=tetration.cpp"></script><br /><br />
TBBを利用したC++による並列化実装:<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">% g++ tetration_tbb.cpp -o tetration_tbb -ltbb -O3</span><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">% ./tetration_tbb tetration.out -1.5 0.0 0.75 0.75 1024 1024</span><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">% ./image_tetration.py tetration.out tetration.png</span><br /><br />
<script src="https://gist.github.com/2689086.js?file=tetration_tbb.cpp"></script><br /><br />
C++の実行ファイルで生成された出力ファイルを画像ファイルに変換するスクリプトを以下に示す。<br /><br />
<script src="https://gist.github.com/2689099.js?file=image_tetration.py"></script><br />
</span>noxhttp://www.blogger.com/profile/01142433117618538300noreply@blogger.com2tag:blogger.com,1999:blog-6153301078335336992.post-33997218546802297352012-03-31T01:20:00.004+09:002012-03-31T13:57:28.507+09:00IPythonの埋め込みプロットが素晴らしい<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi3F7W2dxEiRgZHJSQMeU_MOaiVRWYCsG3q2FgTXqYfrUSuf_ZPUlK8VwNnHujVjaGYZ6wTkDfBWgo4_yI162uUI_MMExR56ZgTUTXiZ30MsQySNE4qjW-h7I83f2ZMXAyWrCHMkHOddzIk/s1600/ipython_mandelbrot.png"><img style="float:left; margin:0 10px 10px 0;cursor:pointer; cursor:hand;width: 400px; height: 354px;" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi3F7W2dxEiRgZHJSQMeU_MOaiVRWYCsG3q2FgTXqYfrUSuf_ZPUlK8VwNnHujVjaGYZ6wTkDfBWgo4_yI162uUI_MMExR56ZgTUTXiZ30MsQySNE4qjW-h7I83f2ZMXAyWrCHMkHOddzIk/s400/ipython_mandelbrot.png" border="0" alt=""id="BLOGGER_PHOTO_ID_5725920357692947746" /></a>先日、<a href="http://partake.in/events/ac0fcc7d-a289-4e2a-bb8e-1965aab8b17b">Tokyo.SciPy #3</a>に参加して、SciPyの生みの親である<a href="http://codeconf.s3.amazonaws.com/2011/pycodeconf/talks/PyCodeConf2011%20-%20Travis%20Oliphant.pdf">Travis Olipant氏のセッション</a>を聞いた。その時に最近の<a href="http://ipython.org/">IPython</a>では埋め込みプロットができることを知ったので早速入れてみたところ、これがとてもクールでカッコよかったのでここで紹介したい。<br /><br />
埋め込みプロットはIPythonのバージョン0.11からできるようになっている。ただ、自分が使っているLinux環境では普通に持ってくるとバージョン0.10が入ってしまうので使えない。そこで、IPythonの公式サイトから最新版のバージョン0.12を持ってきて入れてみた。<br /><br />
ただ、埋め込みプロットができるのはターミナル版ではなくQt版のシェルなので、それに関連するライブラリなども一緒に入れなくてはならない。IPythonの起動時に何々のライブラリが足りないとかいろいろと文句を言われるが、足りないものをyumなりapt-getなりソースコードをコンパイルなりして順次入れていけば使えるようになると思う。<br /><br />
あと、IPythonの設定ファイルの構成が0.10から大きく変わった。以前は .ipython/ 内の pythonrc を利用していたのだが、それがなくなって、代わりに .ipython/profile_default/ 内の ipython_config.py と ipython_qtconsole_config.py が使われるようになった。デフォルトの設定ファイルは以下のコマンドで作ることができる。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">$ ipython profile create</span><br /><br />
そして、埋め込みプロットを使うためのQt版IPythonの起動コマンドは以下の通り。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">$ ipython qtconsole --pylab=inline</span><br /><br />
これでQt版シェルが立ち上がって埋め込みプロットを表示できるようになる。<br />
<span id="fullpost"><br />
せっかくなので実際に埋め込みプロットを使ってみる。ここでは自分のお気に入りのマンデルブロ集合をプロットしてみた。因みに、シェル上での複数行のソースコードの編集も改良されたので入力がかなり楽になった。冒頭にスクリーンショットを貼っておく。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">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))</span><br /><br />
</span>noxhttp://www.blogger.com/profile/01142433117618538300noreply@blogger.com0tag:blogger.com,1999:blog-6153301078335336992.post-45920860563187163842012-02-19T21:43:00.003+09:002012-02-19T21:50:07.464+09:00Android端末をサーバにしてHTML5を使ったお絵かきBBSを作成する<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjvLnKJIgL5tnuQoDLhTp_9ipuORT0_Rm_wGJNMGZawhVCzzFXoyZcQbuAV836slgvMbr50ub1ggEF-khsrI24Ep-ieDhuFoVEFpiKC11gG3lBwLi0psyx7zyCw1HnDF7KNEX47mfedC4nZ/s1600/image_bbs01.png"><img style="float:left; margin:0 10px 10px 0;cursor:pointer; cursor:hand;width: 189px; height: 320px;" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjvLnKJIgL5tnuQoDLhTp_9ipuORT0_Rm_wGJNMGZawhVCzzFXoyZcQbuAV836slgvMbr50ub1ggEF-khsrI24Ep-ieDhuFoVEFpiKC11gG3lBwLi0psyx7zyCw1HnDF7KNEX47mfedC4nZ/s320/image_bbs01.png" border="0" alt=""id="BLOGGER_PHOTO_ID_5710827226913498626" /></a>最近、HTML5に触れる機会があり、その良さが何となく伝わってきたので、何かしら簡単なコードを書いてみたくなった。そこで、ちっちゃいけどリッチ、というギャップを楽しむためにAndroid端末を使うことにした。具体的には、Android端末を<a href="http://code.google.com/p/android-scripting/">SL4A</a>で<a href="http://handasse.blogspot.com/2010/09/pythonandroid5.html">ウェブサーバにして</a>、HTML5をインターフェイスにしたお絵かきBBSをPythonで書いてみた。<br /><br />
BBSでは絵とテキストを扱うことができて、それらはAndroid端末上でSQLiteのデータベースで管理される。利用者の利便を考えて、名前などはクッキーで保存する。3G回線や無線LANなどで接続することを考慮して、IPアドレスも取得できるようにした。また、書き込み時にサーバのAndroid端末が振動して書き込みがあったことを知らせてくれる。因みに、NTTドコモの3G回線で使うためには、グローバルIPが振られるmopera Uなどのサービスが必要で、spモードでは利用できない。<br /><br />
実際に作ってみたプログラムのスクリーンショットを冒頭に入れてみた。必要最低限の機能しかないが、それでも文章と画像を保存できるちゃんとした掲示板だ。草の根BBSのころを考えると隔世の感がある。<br /><br />
以下に今回作成したソースコードを載せておく。こんな短いコードでもちゃんと機能するのが不思議な感じだ。<br />
<span id="fullpost"><br />
<strong>image_bbs.py</strong><br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;"># -*- 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()</span><br />
</span>noxhttp://www.blogger.com/profile/01142433117618538300noreply@blogger.com2tag:blogger.com,1999:blog-6153301078335336992.post-78999089368653784212012-01-01T06:23:00.001+09:002012-01-01T06:34:25.801+09:00恭賀新年今年は2012年です。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">bool leap_year(int y) { return !(y%4)^!(y%100)^!(y%400); }</span><br /><br />
そして、<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">leap_year(2012)</span> が真となるので閏年です。平年と比べて一日増えるわけですが、その一日たりとも無駄にしない充実した一年にしたいと思います。<br /><br />
本年もどうぞ宜しくお願い申し上げます。<br />noxhttp://www.blogger.com/profile/01142433117618538300noreply@blogger.com0tag:blogger.com,1999:blog-6153301078335336992.post-87130265852196506912011-12-23T17:02:00.004+09:002011-12-23T17:26:15.032+09:00円周率を1万桁まで求める<a href="http://ja.wikipedia.org/wiki/%E3%83%9E%E3%83%81%E3%83%B3%E3%81%AE%E5%85%AC%E5%BC%8F">マチンの公式</a>で円周率を1万桁まで求めるプログラムを、Python, Erlang, Haskell, C++で書いてみた。出力結果の最後の数桁はずれているかも。また、実行時間を測ったりもしているが、1万桁程度ならどのコードでも瞬時に求まる。<br /><br />
まずは素直にPythonで。実行時オプションで桁数指定、実行時間測定付き。オプションなしで1万桁まで求める。<br /><br />
<strong>pi.py</strong><br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">#!/usr/bin/env python
import sys, time
N = 10**10000
def arctan(m):
global N
c = N
a = b = c / m
m2 = m * m
s = k = 1
while c:
b /= m2
k += 2
c, s = b / k, -s
a += c * s
return a
def main(args):
global N
if len(args) > 1: N = 10**int(args[1])
t1 = time.time()
pi = str((arctan(5) * 4 - arctan(239)) * 4)
t2 = time.time()
print pi[0] + '.' + pi[1:]
print "Time: %f" % (t2 - t1)
if __name__ == "__main__": main(sys.argv)</span><br /><br />
<strong>実行:</strong><br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">$ ./pi.py</span><br /><br />
次にErlangで求めてみる。<br /><br />
<strong>pi.erl</strong><br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">-module(pi).
-export([pi/0]).
pi()->N=e(10,10000),(a(5,N)*4-a(239,N))*4.
e(B,N)->e(B,N,1).
e(_,0,R)->R;
e(B,N,R)->e(B,N-1,R*B).
a(X,N)->a(X,N div X,N div X,N,1,1).
a(_,A,_,0,_,_)->A;
a(M,A,B,_,S,K)->B_=B div(M*M),C=B_ div(K+2),S_=-S,A_=A+C*S_,a(M,A_,B_,C,S_,K+2).</span><br /><br />
<strong>実行:</strong><br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">$ erl
1> c(pi).
2> pi:pi().</span><br /><br />
そしてHaskellを使う。ソースコードは一行だ。<br /><br />
<strong>pi.hs</strong><br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">n=10^10^4;p=a 5 n*4-a 239 n;a m n=t m(n`div`m)(n`div`m)n 1 3;t _ a _ 0 _ _=a;t m a b c s k=t m(a-y*s)x y(-s)(k+2)where x=b`div`m^2;y=x`div`k</span><br /><br />
<strong>実行:</strong><br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">$ ghci
Prelude> :l pi
Prelude> p*4</span><br /><br />
最後にC++のコード。多倍長演算でちょっと長い。実行時間測定付き。<br />
<span id="fullpost"><br />
<strong>pi.cpp</strong><br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">#include <iostream>
#include <iomanip>
#include <algorithm>
#include <cmath>
#ifdef _MSC_VER
#include <ctime>
inline double get_time()
{
return static_cast<double>(std::clock()) / CLOCKS_PER_SEC;
}
#else
#include <sys/time.h>
inline double get_time()
{
timeval tv;
gettimeofday(&tv, 0);
return tv.tv_sec + 1e-6 * tv.tv_usec;
}
#endif
typedef unsigned long long number;
const int FIG = 8;
const number BASE = 0xffffffffULL;
const number NN = 100000000ULL;
// C = 1.0 / (4 * log10(2)) / FIG;
#define C (0.830482023722/FIG)
#define M (X/FIG+2)
#define N (static_cast<int>(X*C+2))
template<int X> class PI {
private:
number ans[N];
public:
void add(number* a, number* b);
void sub(number* a, number* b);
void div(number* a, number d);
number intdiv(number x, number y);
void mul(number* a, number d);
void atan(number* a, number m);
number or_numbers(number* a);
bool has_bit(number* a);
void decimal(number* a, number* w);
void print();
void calc();
};
// a <- arctan(1 / m)
template<int X>
void PI<X>::atan(number* a, number m)
{
number b[N] = { 1 };
number c[N] = { 1 };
div(b, m);
std::copy(b, b + N, a);
const number m2 = m * m;
for (int i = 1; has_bit(c); i++) {
div(b, m2);
std::copy(b, b + N, c);
div(c, i * 2 + 1);
if (i & 1) sub(a, c);
else add(a, c);
}
}
// a <- a + b
template<int X>
void PI<X>::add(number* a, number* b)
{
for (int i = 0; i < N; i++) {
number x = a[i] + b[i];
if (x & ~BASE) {
a[i] = x & BASE;
int j;
for (j = i - 1; a[j] == BASE; j--) a[j] = 0;
a[j]++;
} else a[i] = x;
}
}
// a <- a - b (a > b)
template<int X>
void PI<X>::sub(number* a, number* b)
{
for (int i = 0; i < N; i++) {
if (a[i] < b[i]) {
a[i] = (BASE + 1) + a[i] - b[i];
int j;
for (j = i - 1; a[j] == 0; j--) a[j] = BASE;
a[j]--;
} else a[i] -= b[i];
}
}
// a <- a / d (d <= BASE)
template<int X>
void PI<X>::div(number* a, number d)
{
number res = 0;
for (int i = 0; i < N; i++) {
res <<= (FIG * 4);
number x = a[i] + res;
number q = x / d;
a[i] = q;
res = x - q * d;
}
}
template<int X>
number PI<X>::intdiv(number x, number y)
{
const number maxbit = ~((number)-1>>1);
number a = 0;
number b = 0;
int cnt = sizeof(number) * 8;
while (cnt--) {
a <<= 1ULL;
if (x & maxbit) a |= 1ULL;
x <<= 1ULL;
b <<= 1ULL;
if (a >= y) { ++b; a -= y; }
}
return b;
}
// a <- a * d (d <= BASE)
template<int X>
void PI<X>::mul(number* a, number d)
{
number q = 0;
for (int i = N - 1; i >= 0; i--) {
number x = a[i] * d + q;
a[i] = x & BASE;
q = x >> (FIG * 4);
}
}
// OR numbers
template<int X>
number PI<X>::or_numbers(number* a)
{
number ret = 0;
for (int i = 0; i < N; i++) ret |= a[i];
return ret;
}
template<int X>
bool PI<X>::has_bit(number* a)
{
for (int i = 0; i < N; i++) if (a[i]) return true;
return false;
}
// hex to decimal number
template<int X>
void PI<X>::decimal(number* a, number* w)
{
number b[N];
std::copy(a, a + N, b);
w[0] = b[0];
b[0] = 0;
for (int i = 1; i < M; i++) {
mul(b, NN);
w[i] = b[0];
b[0] = 0;
}
}
// print answer
template<int X>
void PI<X>::print()
{
number w[M];
decimal(ans, w);
std::cout << std::setw(4) << std::right << w[0] << ".";
for (int i = 1; i < M; i++) {
std::cout << std::setw(FIG) << std::setfill('0') << w[i] << " ";
if(i % 6 == 0) std::cout << std::endl << " ";
}
std::cout << std::endl;
}
// calculate PI
// Machin's formula: 4 * arctan(1/5) - arctan(1/239) = PI / 4
template<int X>
void PI<X>::calc()
{
number x[N];
atan(ans, 5);
mul(ans, 4);
atan(x, 239);
sub(ans, x);
mul(ans, 4);
}
int main()
{
std::ios_base::sync_with_stdio(false);
PI<10000> pi;
double start_time = get_time();
pi.calc();
double end_time = get_time();
pi.print();
std::cerr << "Time: " << end_time - start_time << std::endl;
return 0;
}</span><br /><br />
<strong>実行:</strong><br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">$ c++ pi.cpp -O3 -o pi
$ ./pi</span><br />
</span>noxhttp://www.blogger.com/profile/01142433117618538300noreply@blogger.com0tag:blogger.com,1999:blog-6153301078335336992.post-47954965281243707632011-11-30T00:09:00.008+09:002011-12-01T00:28:40.239+09:00好きな・嫌いなコンピュータ言語のアンケート社内で好きなコンピュータ言語と嫌いなコンピュータ言語のアンケートを取ってみた。ことの発端は<a href="http://twitter.com/#!/foota/status/134991192315277314">自分のツイート</a>なのだが、思いがけず賛同を得ることができたので、実際にアンケートを社内で行うことにしたのだ。各自で社内Wikiを編集してもらってもよかったのだけど、プログラムを作る会社でもあるし、投票システムに興味もあったので、Google App Engine (GAE)で作ってみることにした。以前にGAEで掲示板などを作成していたので、それらを流用して時間をかけずにすぐに仕上げることができた。以下のリンクから実際に投票できる。<br /><br />
<div style="text-align: center;"><strong><a href="http://opinion-poll.appspot.com/vote?id=lang">好きな・嫌いなコンピュータ言語のアンケート</a></strong></div><br /><br />
表示されているコンピュータ言語一覧から、最も好きな言語と最も嫌いな言語をそれぞれ一つずつ、それ以外にも好きな言語や嫌いな言語があれば、それらも選択(複数可)する。投票したい言語が一覧にない場合は「言語の追加」で自由に追加できるようになっている。言語を選択して投票ボタンを押すと、投票数とそのグラフが更新される。表示されるグラフには、「最も好きな・最も嫌いな言語」と「好きな・嫌いな言語」の2つがあり、青のグラフが「好きな言語」、赤のグラフが「嫌いな言語」を表している。「好きな・嫌いな言語」のグラフについては、「最も好きな・最も嫌いな言語」と「好きな・嫌いな言語」を足し合わせた票数となっている。<br /><br />
また、コメントも投稿できるようになっていて、自分の好きな言語を啓蒙するのも良し、嫌いな言語について文句を言うのもアリだ。自由に書き込める。<br />
<span id="fullpost"><br />
現時点での社内での票は以下のグラフに示す通りで、幅広く票が入っているように思う。傾向としては、C, C++, C#, ASM, Python, Rubyあたりが人気で、Visual Basic, Java, Perlが不人気みたいだ。<br /><br />
<a onblur="try {parent.deselectBloggerImageGracefully();} catch(e) {}" href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgdgfsXNrNsrQtjxi5QASZXDatMx9GldQSx7-1M6t66ASFBsoJ-7YgCyqWDh6eg0_2RLfUlsgtkquanElT7kt_7kOIl7dzMq2FIJXNVtEZt49dsj81xAZwAcG_m_RmxQN4k18VkJN04jwtO/s1600/lang01.png"><img style="display:block; margin:0px auto 10px; text-align:center;cursor:pointer; cursor:hand;width: 400px; height: 92px;" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgdgfsXNrNsrQtjxi5QASZXDatMx9GldQSx7-1M6t66ASFBsoJ-7YgCyqWDh6eg0_2RLfUlsgtkquanElT7kt_7kOIl7dzMq2FIJXNVtEZt49dsj81xAZwAcG_m_RmxQN4k18VkJN04jwtO/s400/lang01.png" border="0" alt=""id="BLOGGER_PHOTO_ID_5680438727750169282" /></a><br />
<a onblur="try {parent.deselectBloggerImageGracefully();} catch(e) {}" href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi162z2s2t7KayrAbazw_GKrgIDGvu1lHUxqbQSaLdX5h7RYmF2DDOWK4KYsRIUERNPVZeGyV67JKrZQsW7DI6DDosvm5uDW8fdH_TUIZa2BMCoa-G4pG4t_HItJKglQxpg9R2eNtINZAFy/s1600/lang02.png"><img style="display:block; margin:0px auto 10px; text-align:center;cursor:pointer; cursor:hand;width: 400px; height: 99px;" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi162z2s2t7KayrAbazw_GKrgIDGvu1lHUxqbQSaLdX5h7RYmF2DDOWK4KYsRIUERNPVZeGyV67JKrZQsW7DI6DDosvm5uDW8fdH_TUIZa2BMCoa-G4pG4t_HItJKglQxpg9R2eNtINZAFy/s400/lang02.png" border="0" alt=""id="BLOGGER_PHOTO_ID_5680441594636784130" /></a><br />
今回作成したコードの簡単な説明を書いておく。上述の「好きな・嫌いなコンピュータ言語」だけでなく、初期変数を設定することで、「好きな・嫌いな動物」「好きな・嫌いな食べ物」などのように、簡単に任意のアンケートを作成できるようになっている。グラフについては<a href="http://code.google.com/apis/chart/">Google Chart API</a>を利用しているが、1,000ピクセルを超えることができないので、表示数を1,000ピクセル以内に抑える処理を入れている。また、頻繁にデータベースにアクセスしないようにmemcacheも利用している。詳細についてはコードを読んで欲しい。<br /><br />
以下に今回作成したGAEのコードを示しておく。<br /><br />
<strong>main.py</strong><br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""main.py Copyright (C) 2011 nox"""
import os, cgi, re, math, datetime
import urllib, Cookie
import wsgiref.handlers
from google.appengine.api import memcache
from google.appengine.ext import webapp
from google.appengine.ext import db
from google.appengine.ext.webapp import template
# 許可する投票タイトル.
TITLES = { "test": u"投票テスト",
"lang": u"好きな・嫌いなコンピュータ言語のアンケート" }
# 表示メッセージ.
MESSAGES = { "test": u"投票をお願いします。",
"lang": u"最も好きな言語と最も嫌いな言語をそれぞれ1つ、他に好きな言語や嫌いな言語があればそれも選択(複数可)して投票をお願いします。" }
# 投票する属性. (タイトル / love / like / hate / dislike)
PROPS = { "test": (u"テスト", u"最高", u"良い", u"最低", u"悪い"),
"lang": (u"言語", u"最も好き", u"好き", u"最も嫌い", u"嫌い") }
# チャートのタイトル. (love and hate / like and dislike [including love and hate])
CHARTS = { "test": (u"最高・最低のテスト", u"良い・悪いテスト"),
"lang": (u"最も好きな・最も嫌いな言語", u"好きな・嫌いな言語") }
# 最大表示件数
TARGET_LIMIT = 300
COMMENT_LIMIT = 100
class MainPage(webapp.RequestHandler):
def get(self):
self.response.out.write(u"""
<html>
<head>
<title>Opinion Poll</title>
</head>
<body>
<div><a href="vote">投票システム</a></div>
</body>
</html>""")
# 投票データ.
class VoteData(db.Model):
id = db.StringProperty(multiline=False)
content = db.StringProperty(multiline=False)
love = db.IntegerProperty()
like = db.IntegerProperty()
hate = db.IntegerProperty()
dislike = db.IntegerProperty()
# コメントデータ.
class CommentData(db.Model):
id = db.StringProperty(multiline=False)
username = db.StringProperty(multiline=False)
content = db.TextProperty()
date = db.DateTimeProperty(auto_now_add=True)
# 投票システム・表示.
class Vote(webapp.RequestHandler):
def get(self):
global TITLES, TARGET_LIMIT
id = self.request.get("id")
if id not in TITLES.keys(): return
# データベースまたはMemcacheからデータの読み込み.
cache_votes_key = "votes_%s" % id
try:
votes = memcache.get(cache_votes_key)
except:
try: memcache.delete(cache_votes_key)
except: pass
votes = None
if votes is None:
votes = VoteData.all().order("content").filter("id =", id).fetch(TARGET_LIMIT)
try: memcache.add(cache_votes_key, votes)
except: pass
cache_comments_key = "comments_%s" % id
try:
comments = memcache.get(cache_comments_key)
except:
try: memcache.delete(cache_comments_key)
except: pass
comments = None
if comments is None:
comments = CommentData.all().order("-date").filter("id =", id).fetch(COMMENT_LIMIT)
try: memcache.add(cache_comments_key, comments)
except: pass
# チャート表示 (1,000ピクセル以下にデータ数を収める).
target2, target4 = [], []
target2_max, target4_max = 0, 0
for v in votes:
v.love = int(v.love) if v.love else 0
v.hate = int(v.hate) if v.hate else 0
v.like = int(v.like) if v.like else 0
v.dislike = int(v.dislike) if v.dislike else 0
if v.love or v.hate or v.like or v.dislike:
target4.append((v.content, v.love + v.like, v.hate + v.dislike))
target4_max = max(target4_max, v.love + v.like, v.hate + v.dislike)
if v.love or v.hate:
target2.append((v.content, v.love, v.hate))
target2_max = max(target2_max, v.love, v.hate)
target2_max = (target2_max - 1) / 10 if target2_max > 0 else 0
target2_max = (target2_max + 1) * 10
target4_max = (target4_max - 1) / 10 if target4_max > 0 else 0
target4_max = (target4_max + 1) * 10
n = 0
while True:
target2 = filter(lambda t: max(t[1:]) > n, target2)
if len(target2) <= 16: break
n += 1
if target2:
text = ""
if n > 0: text = u" (%d票以上)" % (n + 1)
fig = 10**(int(math.log10(target2_max)))
if target2_max / fig < 5: fig /= 2
chart_scale = "|".join([str(i) for i in range(0, target2_max + 1, fig)])
chart_target2_url = u"http://chart.apis.google.com/chart?chs=%dx200&chd=t:%s|%s&chds=0,%d&cht=bvg&chco=4d89f9,f9894d&chxt=y,x&chxl=0:|%s|1:|%s&chtt=%s%s" % (len(target2) * 60 + 30, ",".join([str(t[1]) for t in target2]), ",".join([str(t[2]) for t in target2]), target2_max, chart_scale, "|".join([urllib.quote(t[0].encode('utf-8')) for t in target2]), CHARTS[id][0], text)
else: chart_target2_url = ""
n = 0
while True:
target4 = filter(lambda t: max(t[1:]) > n, target4)
if len(target4) <= 16: break
n += 1
if target4:
text = ""
if n > 0: text = u" (%d票以上)" % (n + 1)
fig = 10**(int(math.log10(target4_max)))
if target4_max / fig < 5: fig /= 2
chart_scale = "|".join([str(i) for i in range(0, target4_max + 1, fig)])
chart_target4_url = u"http://chart.apis.google.com/chart?chs=%dx200&chd=t:%s|%s&chds=0,%d&cht=bvg&chco=4d89f9,f9894d&chxt=y,x&chxl=0:|%s|1:|%s&chtt=%s%s" % (len(target4) * 60 + 30, ",".join([str(t[1]) for t in target4]), ",".join([str(t[2]) for t in target4]), target4_max, chart_scale, "|".join([urllib.quote(t[0].encode('utf-8')) for t in target4]), CHARTS[id][1], text)
else: chart_target4_url = ""
# 挨拶.
H = (datetime.datetime.utcnow() + datetime.timedelta(hours=9)).hour
if H >= 5 and H < 11: greeting = u"おはようございます。"
elif H >= 11 and H < 17: greeting = u"こんにちは。"
elif H >= 17 or H < 5: greeting = u"こんばんは。"
else: greeting = u"こんにちは。"
greeting += MESSAGES[id]
# ウェブ上で表示させるデータ.
template_values = {
"id": id,
"title": TITLES[id],
"target": PROPS[id][0],
"love": PROPS[id][1],
"like": PROPS[id][2],
"hate": PROPS[id][3],
"dislike": PROPS[id][4],
"votes": votes,
"comments": comments,
"chart_target2_url": chart_target2_url,
"chart_target4_url": chart_target4_url,
"greeting": greeting,
}
path = os.path.join(os.path.dirname(__file__), "vote.html")
self.response.out.write(template.render(path, template_values))
# 投票ターゲットの追加.
class AddVote(webapp.RequestHandler):
def post(self):
vote = VoteData()
vote.id = self.request.get("id")
try:
try: memcache.delete("votes_%s" % vote.id)
except: pass
except:
pass
if not self.request.get("content").strip():
self.redirect("/vote?id=%s" % vote.id)
return
content = cgi.escape(self.request.get("content").strip()[:100])
v = VoteData.all().filter("id =", vote.id).filter("content =", content).get()
if not v:
vote.content = content
vote.put()
self.redirect("/vote?id=%s" % vote.id)
# 投票.
class PostVote(webapp.RequestHandler):
def post(self):
vote = VoteData()
vote.id = self.request.get("id")
love = self.request.get_all("love")
like = self.request.get_all("like")
hate = self.request.get_all("hate")
dislike = self.request.get_all("dislike")
try:
try: memcache.delete("votes_%s" % vote.id)
except: pass
except:
pass
if not (love or like or hate or dislike):
self.redirect("/vote?id=%s" % vote.id)
return
for d in love:
v = VoteData.all().filter("id =", vote.id).filter("content =", d).get()
if not v.love: v.love = 1
else: v.love = int(v.love) + 1
v.put()
for d in like:
v = VoteData.all().filter("id =", vote.id).filter("content =", d).get()
if not v.like: v.like = 1
else: v.like = int(v.like) + 1
v.put()
for d in hate:
v = VoteData.all().filter("id =", vote.id).filter("content =", d).get()
if not v.hate: v.hate = 1
else: v.hate = int(v.hate) + 1
v.put()
for d in dislike:
v = VoteData.all().filter("id =", vote.id).filter("content =", d).get()
if not v.dislike: v.dislike = 1
else: v.dislike = int(v.dislike) + 1
v.put()
self.redirect("/vote?id=%s" % vote.id)
# コメントの投稿.
class PostComment(webapp.RequestHandler):
def insert_whitespace(self, matchobj):
return re.sub("(?<!<br) ", "&nbsp;", matchobj.group(0).replace("\t", " "))
def post(self):
global TARGET_LIMIT
comment = CommentData()
comment.id = self.request.get("id")
try:
try: memcache.delete("comments_%s" % comment.id)
except: pass
except:
pass
if not self.request.get("content").strip():
self.redirect("/vote?id=%s" % comment.id)
return
if self.request.get("username").strip():
comment.username = cgi.escape(self.request.get("username").strip()[:255])
else:
comment.username = "Anonymous"
comment.content = cgi.escape(self.request.get("content").strip()[:2000])
comment.content = re.sub("\r\n|\r|\n", "<br />", comment.content)
comment.content = re.sub("\[link\](.*?)\[/link\]", "<a href=\"\g<1>\">\g<1></a>", comment.content)
comment.content = re.sub("\[code\](.*?)\[/code\]", self.insert_whitespace, comment.content)
comment.content = re.sub("\[code\](.*?)\[/code\]", "<span style=\"font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;\">\g<1></span>", comment.content)
comment.date = comment.date + datetime.timedelta(hours=9)
comment.put()
self.redirect("/vote?id=%s" % comment.id)
def main():
application = webapp.WSGIApplication([("/", MainPage),
("/vote", Vote),
("/vote/add", AddVote),
("/vote/post", PostVote),
("/vote/comment", PostComment),
],
debug=False)
wsgiref.handlers.CGIHandler().run(application)
if __name__ == "__main__": main()</span><br /><br />
<strong>vote.html</strong><br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;"><?xml version="1.0" encoding="UTF-8"?>
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.1//EN" "http://www.w3.org/TR/xhtml11/DTD/xhtml11.dtd">
<html xmlns="http://www.w3.org/1999/xhtml" xml:lang="ja">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<link type="text/css" rel="stylesheet" href="/stylesheets/main.css" />
<title>{{ title }}</title>
</head>
<body>
<script type="text/javascript">
function resize_textarea(ev){
var textarea = ev.target || ev.srcElement;
var value = textarea.value;
var lines = 1;
for (var i = 0, l = value.length; i < l; i++)
if (value.charAt(i) == '\n') lines++;
if (lines > 5) {
if (lines > 20) textarea.setAttribute("rows", 20);
else textarea.setAttribute("rows", lines);
} else textarea.setAttribute("rows", 5);
}
</script>
<div style="font-weight: bold; font-size: large; padding-bottom: 20px;">{{ title }}</div>
<div>
<p>{{ greeting }}</p>
</div>
<form action="/vote/add" method="post">
<input type="text" name="content" size="20" maxlength="100"></input>
<input type="submit" value="{{ target }}の追加" />
<input type="hidden" name="id" value="{{ id }}">
</form>
<form action="/vote/post" method="post">
<div>
<table border="0" width="400">
<th>{{ target }}</th><th>{{ love }}</th><th>{{ like }}</th><th>{{ hate }}</th><th>{{ dislike }}</th>
{% for vote in votes %}
<tr>
<td><a href="http://ja.wikipedia.org/wiki/{{ vote.content }}">{{ vote.content }}</a></td>
<td align="center"><input type="radio" name="love" value="{{ vote.content }}" /> {% if vote.love %} {{ vote.love }} {% else %} 0 {% endif %}</td>
<td align="center"><input type="checkbox" name="like" value="{{ vote.content }}" /> {% if vote.like %} {{ vote.like }} {% else %} 0 {% endif %}</td>
<td align="center"><input type="radio" name="hate" value="{{ vote.content }}" /> {% if vote.hate %} {{ vote.hate }} {% else %} 0 {% endif %}</td>
<td align="center"><input type="checkbox" name="dislike" value="{{ vote.content }}" /> {% if vote.dislike %} {{ vote.dislike }} {% else %} 0 {% endif %}</td>
</tr>
{% endfor %}
<tr align=center>
<td colspan=5><input type="submit" value="投票" style="width: 100px; height: 25px" /></td>
</tr>
</table>
<input type="hidden" name="id" value="{{ id }}">
</div>
</form><br />
<div>
{% if chart_target2_url %}
<p><img src="{{ chart_target2_url }}" /></p>
{% endif %}
{% if chart_target4_url %}
<p><img src="{{ chart_target4_url }}" /></p>
{% endif %}
</div>
コメント記入: <span style="font-size:70%;color:#666666;">リンク: [link]~[/link], ソースコード: [code]~[/code], 1,000文字まで</span><br />
<form action="/vote/comment" method="post">
<div><textarea name="content" rows="5" cols="60" onkeyup="resize_textarea(event)"></textarea></div>
名前: <input type="text" name="username" size="31" maxlength="255" value="{{ username }}"></input> <input type="submit" value="コメント" /> <span style="font-size:70%;color:#666666;">無記名で匿名投稿</span>
<input type="hidden" name="id" value="{{ id }}">
</form>
<div>
{% for comment in comments %}
{{ comment.username }} [{{ comment.date.year }}年{{ comment.date.month }}月{{ comment.date.day }}日{{ comment.date.hour }}時{{ comment.date.minute }}分{{ comment.date.second }}秒]:
<blockquote>{{ comment.content }}</blockquote>
{% endfor %}
</div>
<!-- ここより下は変更しないでください。不具合やご意見などがありましたら、下記のnoxのリンク先で報告していただけると嬉しいです。 -->
<div><img src="http://code.google.com/appengine/images/appengine-silver-120x30.gif" alt="Powered by Google App Engine" /> <span style="font-size:70%;color:#666666;"><a href="http://handasse.blogspot.com/">ソースコードおよびその解説</a></span><br /><span style="font-size:70%;color:#666666;">Copyright &copy; 2011 <a href="http://handasse.blogspot.com/">nox</a>. All rights reserved.</span></div>
</body>
</html></span><br />
</span>noxhttp://www.blogger.com/profile/01142433117618538300noreply@blogger.com0tag:blogger.com,1999:blog-6153301078335336992.post-52645640193183844842011-10-19T00:07:00.006+09:002011-10-19T21:22:40.982+09:00コンピュータにラストワンを解かせてみた先日、地元のまつりがあった。そこの会場となっている区民館では子供用ボーリングやプラ板遊びなどの催しが開かれていたのだが、その中の一つに「ひとりぼっち」というゲームを行なっているのを見つけたので、一緒にいた小4の自分の娘に「このゲーム面白そうだね」と言ったところ、「ラストワンでしょ」と返され、よく知っているゲームとのことだった。遊び方は次の通り。全部で33のマスがあり、真ん中以外の32個のマスにコマを置く。コマは縦・横方向に隣のコマを一つ飛び越えて進むことができる。飛び越えられたコマは取り除かれる。これを繰り返して、最後のコマを真ん中のマスに置けばクリアとなる。32コマを使うステージ以外にもT字型やピラミッド型にコマを配置するなど、さまざまなバリエーションがある。<br /><br />
<a onblur="try {parent.deselectBloggerImageGracefully();} catch(e) {}" href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiijFlu-keJv8dxPlQmyGJfrLnhz4oBvU_X1xphqpRoGx9MD8b2z-Bmmvu0aD9K1WF8Jl4-0Qeohp25DukTHE81iWs5F6TsL1deK8V7EHKBG8HuzDsxsUnJDXvIEaOYmekm1Vcr-LuAbH4Z/s1600/lastone_ex.png"><img style="display:block; margin:0px auto 10px; text-align:center;cursor:pointer; cursor:hand;width: 400px; height: 176px;" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiijFlu-keJv8dxPlQmyGJfrLnhz4oBvU_X1xphqpRoGx9MD8b2z-Bmmvu0aD9K1WF8Jl4-0Qeohp25DukTHE81iWs5F6TsL1deK8V7EHKBG8HuzDsxsUnJDXvIEaOYmekm1Vcr-LuAbH4Z/s400/lastone_ex.png" border="0" alt=""id="BLOGGER_PHOTO_ID_5665081540739080946" /></a><br />
ネットで調べたところ、<a href="http://blogs.yahoo.co.jp/showagamer/64820991.html">1985年にバンダイから同名の電子ゲームが出ていた</a>ようだ。リンク先では動画もあるのでそれを見てもらえば遊び方は一目瞭然だろう。普通に解いても面白そうだが、コンピュータを使って解かせたくなったので、早速プログラムを作ってみた。<br /><br />
余談だが、子供のころ、コンピュータは何でも解いてくれる魔法の箱だと思っていた。円周率を何桁まで解いたと聞いて、ワクワクしたものだ。ラストワンのようなゲームをコンピュータに解かせたいという気持ちは、子供のころのコンピュータに対する憧れから来るのかもしれない。<br />
<span id="fullpost"><br />
閑話休題。ラストワンは33マスあり、コマが「ある」「ない」の2状態のみを持つので、33ビットを使って盤面の状態を表すことができる。そして、ある状態から取り得る次の状態をすべて求めて、それを優先順位付きキューに入れて探索した。キューが空になるか、ゴールの状態が見つかるまで探索を行う。優先順位付きキューにしたのは探索を短縮するためにヒューリスティックな関数でスコアを決めようと思ったからなんだけど、32コマあるステージも含め、どれも一瞬で解けてしまうので、どうやら高度なヒューリスティック関数を用意しなくても問題なかったらしい。<br /><br />
<a onblur="try {parent.deselectBloggerImageGracefully();} catch(e) {}" href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiWNYbF1V_7oEUy9cWCz1zM81R1jeB3MIb2YgM6OtdGqzW_PzabMheavUNqbk6oWgWdOogbaa1LVh0zUH9h-8QrtIraZW48LnIDp8E9KbvDWu2vGuVRCCx8ZDft0KiCyRGjjwLO0iLj5BPI/s1600/ex_t.png"><img style="display:block; margin:0px auto 10px; text-align:center;cursor:pointer; cursor:hand;width: 180px; height: 180px;" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiWNYbF1V_7oEUy9cWCz1zM81R1jeB3MIb2YgM6OtdGqzW_PzabMheavUNqbk6oWgWdOogbaa1LVh0zUH9h-8QrtIraZW48LnIDp8E9KbvDWu2vGuVRCCx8ZDft0KiCyRGjjwLO0iLj5BPI/s400/ex_t.png" border="0" alt=""id="BLOGGER_PHOTO_ID_5664849589974839794" /></a><br />
上図のように初期配置としてT字型にコマを置く場合、以下に示すような入力データを用意して、これをプログラムに渡す。<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">o</span> (オー)でコマのあるマスを表し、<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">.</span> (ドット)で空のマスを表している。<br /><br />
<strong>lastone.dat</strong><br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;"> ...
...
.......
..ooo..
...o...
.o.
...</span><br /><br />
入力データのファイル名が lastone.dat であれば、以下のように実行すればよい。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">% ./lastone lastone.dat</span><br /><br />
実行結果は以下のようになる。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;"> ...
...
.......
..ooo..
...o...
.o.
...
...
...
.......
.o..o..
...o...
.o.
...
...
...
.......
.o.oo..
.......
...
...
...
...
.......
.oo....
.......
...
...
...
...
.......
...o...
.......
...
...</span><br /><br />
最後に作成したコードを示す。<br /><br />
<strong>lastone.cpp</strong><br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">// Last One solver by nox, 14 Oct. 2011
#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <queue>
#include <algorithm>
#include <utility>
using namespace std;
typedef unsigned long long state_t;
typedef int score_t;
const int N = 33; // number of cells
class LastOne {
private:
vector<vector<pair<int, int> > > D; // table of directions
state_t start, goal;
vector<state_t> ans; // solution
public:
LastOne(const string datafile="");
~LastOne() { }
void read(const string datafile);
int movable(state_t state);
void next(state_t state, vector<state_t>& neighbors);
bool search();
void show(state_t state);
void answer();
};
// initialize and make table
LastOne::LastOne(const string datafile)
{
if (!datafile.empty()) read(datafile);
vector<pair<int, int> > v;
v.push_back(make_pair(1, 2));
v.push_back(make_pair(3, 8));
D.push_back(v);
v.clear();
v.push_back(make_pair(4, 9));
D.push_back(v);
v.clear();
v.push_back(make_pair(1, 0));
v.push_back(make_pair(5, 10));
D.push_back(v);
v.clear();
v.push_back(make_pair(4, 5));
v.push_back(make_pair(8, 15));
D.push_back(v);
v.clear();
v.push_back(make_pair(9, 16));
D.push_back(v);
v.clear();
v.push_back(make_pair(4, 3));
v.push_back(make_pair(10, 17));
D.push_back(v);
v.clear();
v.push_back(make_pair(7, 8));
v.push_back(make_pair(13, 20));
D.push_back(v);
v.clear();
v.push_back(make_pair(8, 9));
v.push_back(make_pair(14, 21));
D.push_back(v);
v.clear();
v.push_back(make_pair(3, 0));
v.push_back(make_pair(7, 6));
v.push_back(make_pair(9, 10));
v.push_back(make_pair(15, 22));
D.push_back(v);
v.clear();
v.push_back(make_pair(4, 1));
v.push_back(make_pair(8, 7));
v.push_back(make_pair(10, 11));
v.push_back(make_pair(16, 23));
D.push_back(v);
v.clear();
v.push_back(make_pair(5, 2));
v.push_back(make_pair(9, 8));
v.push_back(make_pair(11, 12));
v.push_back(make_pair(17, 24));
D.push_back(v);
v.clear();
v.push_back(make_pair(10, 9));
v.push_back(make_pair(18, 25));
D.push_back(v);
v.clear();
v.push_back(make_pair(11, 10));
v.push_back(make_pair(19, 26));
D.push_back(v);
v.clear();
v.push_back(make_pair(14, 15));
D.push_back(v);
v.clear();
v.push_back(make_pair(15, 16));
D.push_back(v);
v.clear();
v.push_back(make_pair(8, 3));
v.push_back(make_pair(14, 13));
v.push_back(make_pair(16, 17));
v.push_back(make_pair(22, 27));
D.push_back(v);
v.clear();
v.push_back(make_pair(9, 4));
v.push_back(make_pair(15, 14));
v.push_back(make_pair(17, 18));
v.push_back(make_pair(23, 28));
D.push_back(v);
v.clear();
v.push_back(make_pair(10, 5));
v.push_back(make_pair(16, 15));
v.push_back(make_pair(18, 19));
v.push_back(make_pair(24, 29));
D.push_back(v);
v.clear();
v.push_back(make_pair(17, 16));
D.push_back(v);
v.clear();
v.push_back(make_pair(18, 17));
D.push_back(v);
v.clear();
v.push_back(make_pair(13, 6));
v.push_back(make_pair(21, 22));
D.push_back(v);
v.clear();
v.push_back(make_pair(14, 7));
v.push_back(make_pair(22, 23));
D.push_back(v);
v.clear();
v.push_back(make_pair(15, 8));
v.push_back(make_pair(21, 20));
v.push_back(make_pair(23, 24));
v.push_back(make_pair(27, 30));
D.push_back(v);
v.clear();
v.push_back(make_pair(16, 9));
v.push_back(make_pair(22, 21));
v.push_back(make_pair(24, 25));
v.push_back(make_pair(28, 31));
D.push_back(v);
v.clear();
v.push_back(make_pair(17, 10));
v.push_back(make_pair(23, 22));
v.push_back(make_pair(25, 26));
v.push_back(make_pair(29, 32));
D.push_back(v);
v.clear();
v.push_back(make_pair(18, 11));
v.push_back(make_pair(24, 23));
D.push_back(v);
v.clear();
v.push_back(make_pair(19, 12));
v.push_back(make_pair(25, 24));
D.push_back(v);
v.clear();
v.push_back(make_pair(22, 15));
v.push_back(make_pair(28, 29));
D.push_back(v);
v.clear();
v.push_back(make_pair(23, 16));
D.push_back(v);
v.clear();
v.push_back(make_pair(24, 17));
v.push_back(make_pair(28, 27));
D.push_back(v);
v.clear();
v.push_back(make_pair(27, 22));
v.push_back(make_pair(31, 32));
D.push_back(v);
v.clear();
v.push_back(make_pair(28, 23));
D.push_back(v);
v.clear();
v.push_back(make_pair(29, 24));
v.push_back(make_pair(31, 20));
D.push_back(v);
}
// read LastOne data
void LastOne::read(const string datafile)
{
fstream fs(datafile.c_str(), ios_base::in);
string line;
goal = 0x10000;
start = 0ULL;
while (getline(fs, line)) {
for (string::iterator p = line.begin(); p != line.end(); ++p) {
if (*p == 'o') { start <<= 1; start |= 1ULL; }
else if (*p == '.') start <<= 1;
}
}
fs.close();
}
// number of movable pieces
int LastOne::movable(state_t state)
{
int s = 0;
for (int i = 0; i < N; i++)
if ((state >> i) & 1)
for (vector<pair<int, int> >::iterator p = D[i].begin(); p != D[i].end(); ++p)
s += ((state >> p->first) & 1) & !((state >> p->second) & 1);
return s;
}
// next states
void LastOne::next(state_t state, vector<state_t>& neighbors)
{
neighbors.clear();
for (int i = 0; i < N; i++)
if ((state >> i) & 1)
for (vector<pair<int, int> >::iterator p = D[i].begin(); p != D[i].end(); ++p)
if (((state >> p->first) & 1) & !((state >> p->second) & 1))
neighbors.push_back(state ^ (1ULL << p->first) ^ (1ULL << p->second) ^ (1ULL << i));
}
// search solution
bool LastOne::search()
{
vector<state_t> checked(1, start);
priority_queue<pair<score_t, vector<state_t> > > queue;
queue.push(make_pair(1, checked));
vector<state_t> neighbors;
while (!queue.empty()) {
pair<score_t, vector<state_t> > q(queue.top());
queue.pop();
state_t last = q.second.back();
if (last == goal) {
ans = q.second;
return true;
}
next(last, neighbors);
for (vector<state_t>::iterator p = neighbors.begin(); p != neighbors.end(); ++p) {
if (binary_search(checked.begin(), checked.end(), *p)) continue;
if (movable(*p) == 0 && *p != goal) continue;
checked.insert(upper_bound(checked.begin(), checked.end(), *p), *p);
vector<state_t> path(q.second);
path.push_back(*p);
queue.push(make_pair(static_cast<score_t>(path.size()), path));
}
}
return false;
}
// show state
void LastOne::show(state_t state)
{
string line;
while (state) {
if (state & 1ULL) line.push_back('o');
else line.push_back('.');
state >>= 1;
}
int n = N - line.size();
for (int i = 0; i < n; i++) line.push_back('.');
reverse(line.begin(), line.end());
cout << " " << line.substr(0, 3) << endl;
cout << " " << line.substr(3, 3) << endl;
cout << line.substr(6, 7) << endl;
cout << line.substr(13, 7) << endl;
cout << line.substr(20, 7) << endl;
cout << " " << line.substr(27, 3) << endl;
cout << " " << line.substr(30, 3) << endl;
}
// show solution
void LastOne::answer()
{
for (vector<state_t>::iterator p = ans.begin(); p != ans.end(); ++p) {
show(*p);
cout << endl;
}
}
int main(int argc, char* argv[])
{
if (argc < 2) {
cerr << "Usage: " << argv[0] << " datafile" << endl;
exit(1);
}
LastOne lastone(argv[1]);
if (lastone.search()) lastone.answer();
else cout << "No solution." << endl;
return 0;
}</span><br />
</span>noxhttp://www.blogger.com/profile/01142433117618538300noreply@blogger.com0tag:blogger.com,1999:blog-6153301078335336992.post-15436165630585950242011-10-14T01:03:00.006+09:002011-10-17T13:15:20.133+09:00オフィスビル内の最短経路<a onblur="try {parent.deselectBloggerImageGracefully();} catch(e) {}" href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi5U8CeBv_VAertcuCO2o3GJL5Rv_vEQPR7YrsK_kLPK2xB0nnl0SMasxuiTsHs3WUwB8Hv8KnqgPZE0JDVTNaof2fVpq64Qr24v2fU-LJGB_xwFb34lkmsUiJFx0XEsSD-SlOIQ8dx2p8a/s1600/office_3f_which_t.png"><img style="float:left; margin:0 10px 10px 0;cursor:pointer; cursor:hand;width: 227px; height: 294px;" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi5U8CeBv_VAertcuCO2o3GJL5Rv_vEQPR7YrsK_kLPK2xB0nnl0SMasxuiTsHs3WUwB8Hv8KnqgPZE0JDVTNaof2fVpq64Qr24v2fU-LJGB_xwFb34lkmsUiJFx0XEsSD-SlOIQ8dx2p8a/s400/office_3f_which_t.png" border="0" alt=""id="BLOGGER_PHOTO_ID_5663008666911204866" /></a>現在のオフィスでは、エントランスからどの通路を通るのが最も近いのか、未だに話題になることがある。エントランスからオフィスにつながるエレベータまでの道が二通りあり、どちらも似たような距離なので同じオフィスに行く人達なのにそこで別々になってしまうこともしばしば。左の図が問題のフロアなんだけど、エントランスから青丸で示したエレベータまで行く必要がある。だったらコンピュータに解かせてしまえばいいじゃない、ということでプログラムを作って最短経路と正確な距離を調べてみた。定規で測ったほうが早いという意見は聞こえません。<br /><br />
以前、<a href="http://handasse.blogspot.com/2009/10/python-2.html">Python: 画像で与えられた迷路に対し2点間の最短経路を求める</a>というブログ記事で、画像から直接最短距離を求めるプログラムを書いたことがあるので、まずはそれを使って調べてみた。図面の邪魔な文字だけ消して、あとはスタート地点とゴール地点の座標を引数で指定するだけだから簡単なんだけど、これだと問題があることに気がついた。というのはCCL (connected component labeling; 連結成分ラベリング)で移動できる領域を調べて、それをA*で解いているのだけど、隣接するノード(ピクセル)を上下左右斜めの8方向のみで調べているから正しい距離を出すことができないのだ。実際に作成した図は以下のようにギザギザになってしまい最短ルートを通らない。なので通り道の分岐や方向が変化する箇所をノードにしてそれをA*で求めることにした。<br />
<span id="fullpost"><br />
<a onblur="try {parent.deselectBloggerImageGracefully();} catch(e) {}" href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgS_ngtO2-Ax55iSzvEkX8vk41kgKeVOFkpcuskYbZx5nDqQ9yBOmo7jBg52FtNc8i9nqu-LfABhJrjhtREBmFKYNjbsGzUJENc8Ulf-jRkPevWQx62zvIZqwid6_W_Yhs397KyPGFgTNKQ/s1600/result_01a_t.png"><img style="display:block; margin:0px auto 10px; text-align:center;cursor:pointer; cursor:hand;width: 228px; height: 295px;" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgS_ngtO2-Ax55iSzvEkX8vk41kgKeVOFkpcuskYbZx5nDqQ9yBOmo7jBg52FtNc8i9nqu-LfABhJrjhtREBmFKYNjbsGzUJENc8Ulf-jRkPevWQx62zvIZqwid6_W_Yhs397KyPGFgTNKQ/s400/result_01a_t.png" border="0" alt=""id="BLOGGER_PHOTO_ID_5663009140727149282" /></a><br />
全部で10個のノード(画像のピクセル座標で指定)を作って調べたところ、以下のような結果になった。左側を通る通路の距離が324.53、右側を通る通路の距離が322.51となり僅かに右側の通路のほうが近かった。ただし、右側の方は途中でフロア案内と少し重なるし、6つあるエレベータのどれを使うかによっても結果が変わってしまうぐらいなので、実質的には同じ距離といってしまっても差し支えないと思う。<br /><br />
<a onblur="try {parent.deselectBloggerImageGracefully();} catch(e) {}" href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjJIh3PgF2XpzHVvJUjVS-B6af9AuKhnVbIohQlt3_tkBB0xBgZ9e1P2meXxM0uLqpPYye4Q2768BGkYBCrHAFRpxl2nGCKjPJAnjUE8h1GN5NofrBgn6VOQDqhdxJyq0QhNKi4_cAdTkaZ/s1600/route_01a_dist_t.png"><img style="display:block; margin:0px auto 10px; text-align:center;cursor:pointer; cursor:hand;width: 228px; height: 295px;" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjJIh3PgF2XpzHVvJUjVS-B6af9AuKhnVbIohQlt3_tkBB0xBgZ9e1P2meXxM0uLqpPYye4Q2768BGkYBCrHAFRpxl2nGCKjPJAnjUE8h1GN5NofrBgn6VOQDqhdxJyq0QhNKi4_cAdTkaZ/s400/route_01a_dist_t.png" border="0" alt=""id="BLOGGER_PHOTO_ID_5663010017921635138" /></a><br />
<a onblur="try {parent.deselectBloggerImageGracefully();} catch(e) {}" href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhk1wYVKdXgiIffbelKkTcw-H0rjI-EuuF9i9jMx2EzgS-ODnSTzvwVKYYLCQ0snVMJbiidClC4hkYL0-EWUPjNM_Yo0PRe_PEislr0p80DJpwvGOgr7e5rY0IyRus7yjnOmtBY9BIEY-HF/s1600/route_01_dist_t.png"><img style="display:block; margin:0px auto 10px; text-align:center;cursor:pointer; cursor:hand;width: 228px; height: 295px;" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhk1wYVKdXgiIffbelKkTcw-H0rjI-EuuF9i9jMx2EzgS-ODnSTzvwVKYYLCQ0snVMJbiidClC4hkYL0-EWUPjNM_Yo0PRe_PEislr0p80DJpwvGOgr7e5rY0IyRus7yjnOmtBY9BIEY-HF/s400/route_01_dist_t.png" border="0" alt=""id="BLOGGER_PHOTO_ID_5663010020264014130" /></a><br />
結局、それぞれの距離にはほとんど差がないので、結論としては「好きな方を使え」ということになるだろう。自分は広い通路の方が好みなので右側の通路を利用している。最後にソースコードとその使い方を示しておく。<br /><br />
nodes.datに、一行を"ピクセルX座標 ピクセルY座標 隣接ノード番号リスト"としてノード数の分だけ記述する。ノード番号は0から始まり、先頭行と最終行のノードはそれぞれスタート地点、ゴール地点を表している。また、フロアの図面画像をinput.pngとして渡すと、最短距離のルートが書きこまれた画像をoutput.pngとして出力する。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">% ./route_solver.py nodes.dat input.png output.png</span><br /><br />
<strong>nodes.dat</strong><br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;"> 27 273 1 3 6
21 128 0 2
125 23 1 9
117 221 0 4 5 6
141 178 3 7
167 178 3 7
195 178 0 3 7
192 110 4 5 6 8
193 91 7 9
151 41 2 8</span><br /><br />
<strong>route_solver.py</strong><br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">#!/usr/bin/env python
# -*- coding: utf-8 -*-
import sys, os, math, heapq, Image, ImageDraw
class Node:
def __init__(self):
self.pos = (0, 0)
self.neighbors = []
self.dist = float("inf")
self.path = []
# A* Search Algorithm
class AStar:
def __init__(self, start, goal):
self.start, self.goal = start, goal
self.start.dist = 0.0
self.start.path = [self.start]
def distance(self, node1, node2):
return math.sqrt((node1.pos[0] - node2.pos[0])**2 + (node1.pos[1] - node2.pos[1])**2)
def heuristic(self, node):
return self.distance(node, self.goal)
def search(self):
queue = []
heapq.heappush(queue, (self.heuristic(self.start), self.start))
while queue:
score, last = heapq.heappop(queue)
if last == self.goal: return last.path, last.dist
for nb in last.neighbors:
newdist = last.dist + self.distance(last, nb)
if nb.dist < newdist: continue
nb.dist = newdist
nb.path = last.path + [nb]
heapq.heappush(queue, (nb.dist + self.heuristic(nb), nb))
return [], 0.0
def draw_path(in_image, out_image, path):
im = Image.open(in_image)
im = im.convert("RGB")
draw = ImageDraw.Draw(im)
draw.line([(int(p.pos[0]), int(p.pos[1])) for p in path], fill=255, width=2)
im.save(out_image)
def main(args):
if len(args) < 4:
print >>sys.stderr, "Usage: %s datafile in_image out_image" % os.path.basename(args[0])
sys.exit(1)
datafile, in_image, out_image = args[1:4];
print "Reading nodes data..."
data = [map(int, l.strip().split()) for l in file(datafile).readlines()]
nodes = [Node() for i in range(len(data))]
for i, node in enumerate(nodes):
node.pos = tuple(data[i][:2])
for j in data[i][2:]:
node.neighbors.append(nodes[j])
print " Number of nodes: %d" % len(nodes)
print "A* searching phase..."
start, goal = nodes[0], nodes[-1]
print " Start: (%d, %d), Goal: (%d, %d)" % (start.pos + goal.pos)
astar = AStar(start, goal)
path, dist = astar.search()
if path:
print " Found route:"
for p in path: print " (%d, %d)" % p.pos
print " Distance: %f" % dist
else:
print " No route."
print "Drawing path on the imagefile..."
draw_path(in_image, out_image, path)
if __name__ == "__main__": main(sys.argv)</span><br /><br />
追記 (2011/10/17): ソースコードが間違っていたので修正。今回の結果には影響なし。<br />
</span>noxhttp://www.blogger.com/profile/01142433117618538300noreply@blogger.com0tag:blogger.com,1999:blog-6153301078335336992.post-49186699121782675412011-09-27T04:06:00.002+09:002015-04-06T11:15:20.696+09:00GoogleのURL短縮サービスgoo.glのAPIをPythonで利用するしばらく前から<a href="http://goo.gl/">goo.gl</a>のAPI (Google URL Shortener API)が使えるようになっていたことは知っていたけど、実際に使ったことなかったので試しに使ってみた。import文を含めなければ3行で短縮URLの作成・復元の自動判別を行うことができた。条件演算子などを使えば1行にまとめられる。簡単。<br /><br />
ただし、正規表現は適当だし、エラー処理はしてないので、ちゃんと使う場合はその辺の処理を行う必要あり。あと、本格的に使うのであれば<a href="http://code.google.com/intl/ja/apis/urlshortener/v1/getting_started.html#APIKey">APIキーを取得</a>すべき。また、<a href="http://code.google.com/p/google-api-python-client/">Google APIs Client Library for Python</a>を<a href="http://code.google.com/p/google-api-python-client/source/browse/samples/urlshortener/urlshortener.py">使う方法</a>もある。<br /><br />
http://handasse.blogspot.com/ を短縮URLに変換:<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">% googl.py http://handasse.blogspot.com/
http://goo.gl/RiIiD</span><br /><br />
http://goo.gl/RiIiD を元のURLに戻す:<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">% googl.py http://goo.gl/RiIiD
http://handasse.blogspot.com/</span><br /><br />
<strong>googl.py</strong><br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">import sys,re
from urllib2 import urlopen as U, Request as R
from json import loads as J
API,URL="https://www.googleapis.com/urlshortener/v1/url",sys.argv[1]
if re.match('http://goo\.gl/.+',URL):print J(U(API+'?shortUrl=%s'%URL).read())['longUrl']
else:print J(U(R(API,'{"longUrl":"%s"}'%URL,{'Content-Type':'application/json'})).read())['id']</span><br /><br />
<strong>追記</strong> (2015/3/26):<br /><br />
現在、API Keyを利用しないとエラーが返ってくるようなので、以下のように修正した。<br /><br />
<script src="https://gist.github.com/foota/45a41e635caeaf0b4aea.js"></script><br />
コード中のYOUR_API_KEYの部分を実際のAPI Keyと差し替えれば利用できる。API Keyの取得するためには、まず、<a href="https://console.developers.google.com/">Google Developers Console</a>にアクセスし、適当なプロジェクトに入り、APIと認証のAPIを選択し、Other popular APIsの中のURL Shortener APIを選ぶ。そこのページでEnable APIのボタンをクリックし、「APIコンソールでレポートを表示する」のリンクが表示されるので、そのページにアクセスする。メニューのAPI Accessを選び、そこのCreate new Browser key...ボタンをクリックする。これでAPI Keyが作成されるので、これを上述のYOUR_API_KEYと差し替えれば利用できるようになる。<br />
noxhttp://www.blogger.com/profile/01142433117618538300noreply@blogger.com2tag:blogger.com,1999:blog-6153301078335336992.post-74260564494328749472011-08-08T00:48:00.003+09:002011-08-08T01:04:39.890+09:00暇潰しで書いたプログラム<a onblur="try {parent.deselectBloggerImageGracefully();} catch(e) {}" href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjU-6ErxPtG7ROPx8lao8o7gY6sa4kQthCIMvzGTd6-1ZdU2hHNBvtYLQOQL1paVv0XsJU0ojcuoiIN252PPCvJNY2D_xhpH-vqm0O09YapVa5H6KQpmao-lJzTWXLBzmAIH12oV7devk7J/s1600/code.png"><img style="float:left; margin:0 10px 10px 0;cursor:pointer; cursor:hand;width: 200px; height: 200px;" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjU-6ErxPtG7ROPx8lao8o7gY6sa4kQthCIMvzGTd6-1ZdU2hHNBvtYLQOQL1paVv0XsJU0ojcuoiIN252PPCvJNY2D_xhpH-vqm0O09YapVa5H6KQpmao-lJzTWXLBzmAIH12oV7devk7J/s320/code.png" border="0" alt=""id="BLOGGER_PHOTO_ID_5638142066670877042" /></a>
ちょっと出先で手持ち無沙汰になったとき、なにかコードが書きたくなることってあると思う。その時に書くプログラムは何でもよくて、言語も選ばない。でも、駅で電車を待っている間とか、注文したメニューが来るまでとか、ノートPCを引っ張り出してまでやりたくはない。そんなとき、Android端末を片手に<a href="http://code.google.com/p/android-scripting/">SL4A</a>でコードを書くっていうのは丁度よい。<br /><br />
というわけで、ちょっと外出したときの暇な時間にSL4AのPythonで書いたコードが以下のコード。本当にくだらないプログラムなわけだけど、暇潰しにはなった。何をするプログラムか興味のある方は読み解いて欲しい。暇潰しぐらいにしかならないけど。まあ、勘のよい人ならファイル名だけでわかってしまうかも。 <br /><br />
Android端末を持っている方なら、SL4Aのメニューを開いてAddを選択し、Scan Barcodeで冒頭のQRコードを読み込むことでソースコードを取得できる。SL4Aで書いたプログラムだけど、Android端末以外でもPythonの動作する環境なら大抵は実行できると思う。<br /><br />
<strong>Zm9yIGkgaW4gcmFuZ2UoMSwxMDEpOnByaW50J0ZpenonKigxLWklMykrJ0J1enonKigxLWklNSlvciBp.py</strong><br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">exec(r"""%s'''gzga*%%%%,hmkl*ocr*nco`fc"a8ajp*mpf*a+\\0+.p%%%s'kormpv"mq.`cqg469gzga*`cqg46,fgamfgqvpkle*mq,rcvj,`cqglcog*]]dkng]]+Y8/1_++')))%%+++''')))"""%(lambda _:(_,_))("exec(''.join(map(lambda _:chr(ord(_)^2),"))</span><br />noxhttp://www.blogger.com/profile/01142433117618538300noreply@blogger.com2tag:blogger.com,1999:blog-6153301078335336992.post-81817016349168492712011-07-11T03:24:00.012+09:002011-07-11T16:32:39.143+09:00Scalaのアクターモデルでマンデルブロ集合を並列計算<a onblur="try {parent.deselectBloggerImageGracefully();} catch(e) {}" href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg9o6DfsDz3yufw-CBfkFPxEuwf-sm-DCGqnClicoYgiXQZcYCsTSBznlZ6vz-3VBtsyu0Oj1N4slvBWkN0QAzjiS8olV3AFYY4bAytAF8NztgknCQyMS0gPKSFileUFvmiCPy8V8itdswv/s1600/mandelbrot.png"><img style="float:left; margin:0 10px 10px 0;cursor:pointer; cursor:hand;width: 200px; height: 200px;" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg9o6DfsDz3yufw-CBfkFPxEuwf-sm-DCGqnClicoYgiXQZcYCsTSBznlZ6vz-3VBtsyu0Oj1N4slvBWkN0QAzjiS8olV3AFYY4bAytAF8NztgknCQyMS0gPKSFileUFvmiCPy8V8itdswv/s200/mandelbrot.png" border="0" alt=""id="BLOGGER_PHOTO_ID_5627801316169846786" /></a>
最近、本格的にScalaを使用するつもりで環境を整えた。ビルドツールとして<a href="https://github.com/harrah/xsbt/wiki">Simple Build Tool (sbt)</a>を利用し、エディタは<a href="https://github.com/aemoncannon/ensime">ensime</a>を入れたEmacsにしている。利用するに当たって<a href="http://aemon.com/file_dump/ensime_manual.html#tth_sEc2.2">.emacsに必要な記述を加えておく</a>こと。初めてのプロジェクトの場合、適当に作成したディレクトリ内でsbtを起動してプロジェクトの雛形を作った後、Emacsを立ち上げて M-x ensime-config-gen で下準備をする。ここまでは初回のみの作業となる。その後は M-x ensime を実行して、src/main/scala/内でコードを書くだけ。これで、コーディング中にタブで補完してくれるし、文法などが間違っている場合に赤でハイライトもしてくれる。コンパイルや実行などは、C-c C-v s とすればsbtがEmacs上で立ち上がるので compile や run するだけ。簡単。<br /><br />
Scalaと言えばアクターでしょ、やっぱり。ということで今回はアクターモデルを使用したコードを書いてみた。自分はScalaをいじるよりも前にErlangをいじっていたので余計にそう考えてしまうのかもしれない。それでもマルチコアやらメニーコアが溢れているこの時代、並行処理の重要性は誰もが認識しているはず。<br /><br />
書いたコードはマンデルブロ集合の計算だ。毎回飽きもせずにまたそれかと思うかもしれないが、並列処理させるのが簡単だし、いつも同じ内容にしておけば比較も容易。なによりもマンデルブロ集合の画像はともて魅力的だと自分では思っているので、今回もマンデルブロ集合だ。冒頭に作成した画像を載せたけど、やっぱり美しいと思う。<br />
<span id="fullpost"><br />
今回のアクターモデルでは、scala.actors.Actorから派生させずに、scala.actors.Actor.actorを使用した。クラスを作るほど大行ではないし、その方がスッキリするからだ。そして、計算のループ部分は末尾再帰とパターンマッチを使っている。関数型言語として使うならやっぱりその方がしっくりくる。ただ、ちゃんと末尾再帰させるように、final修飾子をつけるなど、気をつけること。それと@tailrecアノテーション(scala.annotation.tailrec)は常に付けたほうが良い。@tailrecアノテーションを末尾再帰させたい関数につけておけば、末尾再帰としてコンパイルされなかった時にエラーを返してくれる。アクターを使用しないコードもついでに作成した。画像の色付け部分は<a href="http://d.hatena.ne.jp/a-san/20100213#p1">ScalaでMandelbrot</a>、時間測定は<a href="http://diaspar-journal.blogspot.com/2009/04/blog-post.html#lang-scala">関数の実行時間を計測するには?</a>を参考にさせてもらった。<br /><br />
実行環境としてXeon X5650@2.67GHz * 2 (12コア/24スレッド; メモリ48GB)を使用した。マンデルブロ集合の条件としては、画像の大きさが4,800 x 4,800ピクセルで、(-0.005, -0.005)から(0.005, 0.005)を計算の範囲とし、繰り返し数10,000で10.0を閾値としている。因みにこの範囲は集合に含まれるので計算量としては最大となる(画像としては真っ黒なので面白くはないが)。実行すると最後に画像ファイル(mandelbrot.png)を出力して終了する。<br /><br />
実行した結果、アクターを使用しなかったコードの実行時間が<strong>1007.18秒</strong>だったのに対し、アクターを使用したコードでは<strong>51.636秒</strong>であった。<strong>約20倍の高速化</strong>である。コア数12・スレッド数24であることを考えても、この高速化は十分なのではないだろうか。ところで、アクターを使用しないシングルスレッドでのコードでも以前作成したC++のコード(1067.39秒)より速かった。それに、TBBを使用した並列化したコードでも56.8956秒の計算時間がかかっており、アクターを使用したScalaのコードはそれよりも速い。色の作成の部分など少し異なるところもあるのでそのまま比較することはできないかもしれないが、少なくとも今回の場合ではScalaで十分な速度を得られることは確認できた。余談だが、Telsaのボードを2枚挿ししたマシン上で実行した<a href="http://handasse.blogspot.com/2011/06/cuda-40gpu.html">CUDAのコードでは3.73秒</a>なので次元が違ったりする。<br /><br />
しかし、実際にコードを書くときは注意深くなる必要があるかもしれない。例えば以下のコードがある。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;"> i match {
case MAX_LOOP => -1
case _ if zr2 + zi2 >= TH => i
case _ => {
val zrzi = zr * zi
calcLoop(i + 1, cr, ci, zr2 - zi2 + cr, zrzi + zrzi + ci)
}
}</span><br /><br />
これを以下のように変更する。<br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;"> (i, zr2 + zi2) match {
case (MAX_LOOP, _) => -1
case (_, z) if z >= TH => i
case _ => {
val zrzi = zr * zi
calcLoop(i + 1, cr, ci, zr2 - zi2 + cr, zrzi + zrzi + ci)
}
}</span><br /><br />
これだけのことで実行速度が6倍程遅くなる。これに気付かずにいると本来のパフォーマンスを出すことができない。予想よりもプログラムが遅いと感じた場合は再度コードを確認したほうがいいだろう。<br /><br />
<a href="http://www.amazon.co.jp/gp/product/4822284239/ref=as_li_qf_sp_asin_il?ie=UTF8&tag=footawwwservi-22&linkCode=as2&camp=247&creative=1211&creativeASIN=4822284239"><img style="float:left;margin:0 10px 10px 0;" border="0" src="http://ws.assoc-amazon.jp/widgets/q?_encoding=UTF8&Format=_SL110_&ASIN=4822284239&MarketPlace=JP&ID=AsinImage&WS=1&tag=footawwwservi-22&ServiceVersion=20070822" ></a><img src="http://www.assoc-amazon.jp/e/ir?t=footawwwservi-22&l=as2&o=9&a=4822284239" width="1" height="1" border="0" alt="" style="border:none !important; margin:0px !important;" /><a href="http://www.amazon.co.jp/gp/product/4844327453/ref=as_li_tf_il?ie=UTF8&tag=footawwwservi-22&linkCode=as2&camp=247&creative=1211&creativeASIN=4844327453"><img style="float:right;margin:0 10px 10px 0;" border="0" src="http://ws.assoc-amazon.jp/widgets/q?_encoding=UTF8&Format=_SL110_&ASIN=4844327453&MarketPlace=JP&ID=AsinImage&WS=1&tag=footawwwservi-22&ServiceVersion=20070822" ></a><img src="http://www.assoc-amazon.jp/e/ir?t=footawwwservi-22&l=as2&o=9&a=4844327453" width="1" height="1" border="0" alt="" style="border:none !important; margin:0px !important;float:left;margin:0 10px 10px 0;" />自分は普段は主にC++とPythonを使っていて、目的によって使い分けている。ある程度しっかりした作りで高速に動作させたい場合はC++だし、プロトタイプとして作ったり、結果の解析用コードをサクっと作りたい場合はPythonを使う。Scalaはどのような用途に使えるのか。まず、Javaの資産を使えるのは大きい。ネットワーク処理や画像処理などをある程度本格的につくろうとしたときに便利そうだ。あとは言語の仕様として処理を簡潔に書き表せるという点が良い。ライブラリが強力で簡潔に書ける。この点だけで十分に利用する価値があるだろう。<br /><br />
最後に自分の持っているScalaの書籍を2つ挙げておく。<a href="http://www.amazon.co.jp/gp/product/4822284239/ref=as_li_qf_sp_asin_tl?ie=UTF8&tag=footawwwservi-22&linkCode=as2&camp=247&creative=1211&creativeASIN=4822284239">Scalaプログラミング入門</a><img src="http://www.assoc-amazon.jp/e/ir?t=footawwwservi-22&l=as2&o=9&a=4822284239" width="1" height="1" border="0" alt="" style="border:none !important; margin:0px !important;" />は実践的だし、<a href="http://www.amazon.co.jp/gp/product/4844327453/ref=as_li_qf_sp_asin_tl?ie=UTF8&tag=footawwwservi-22&linkCode=as2&camp=247&creative=1211&creativeASIN=4844327453">Scalaスケーラブルプログラミング</a><img src="http://www.assoc-amazon.jp/e/ir?t=footawwwservi-22&l=as2&o=9&a=4844327453" width="1" height="1" border="0" alt="" style="border:none !important; margin:0px !important;" />は網羅的なので、どちらも持っていて損はないと思う。<br /><br />
作成したコードを以下に示しておく。<br /><br />
<strong>mandelbrot.scala (アクターあり; 並行処理)</strong><br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">// mandelbrot.scala with actor by nox, 2011.07.09
import scala.actors._, Actor._
import scala.annotation.tailrec
import java.io.FileOutputStream
import java.awt.image.BufferedImage
import java.awt.Color
import javax.imageio.ImageIO
import scala.testing.Benchmark
/**
* マンデルブロ集合.
* SW : 画像の幅.
* SH : 画像の高さ.
* T : 集合の上の位置.
* L : 集合の左の位置.
* W : 集合の幅.
* H : 集合の高さ.
* MAX_LOOP : イテレーションの最大回数.
* TH : 閾値.
*/
final class Mandelbrot(SW: Int, SH: Int, T: Double, L: Double, W: Double, H: Double, MAX_LOOP: Int, TH: Double) {
/**
* イテレーションの回数から色を決定して任意のピクセルにセットする.
*/
def setPixel(x: Int, y: Int, a: Any, image: BufferedImage) =
a match {
case -1 => image setRGB(x, y, 0)
case i: Int => image setRGB(x, y, Color HSBtoRGB(i / 100.0f, 1.0f, 1.0f))
}
/**
* マンデルブロ集合の計算部分.
*/
@tailrec
def calcLoop(i: Int, cr: Double, ci: Double, zr: Double, zi: Double): Int = {
val zr2 = zr * zr
val zi2 = zi * zi
i match {
case MAX_LOOP => -1
case _ if zr2 + zi2 >= TH => i
case _ => {
val zrzi = zr * zi
calcLoop(i + 1, cr, ci, zr2 - zi2 + cr, zrzi + zrzi + ci)
}
}
}
/**
* アクター.
* ある幅に存在するピクセルの計算を受け持つ.
*/
def calc = actor {
react {
case (ix: Int, image: BufferedImage) => {
for (iy <- 0 until SH)
setPixel(ix, iy, calcLoop(0, L + (ix.toDouble / SW) * W, T + (iy.toDouble / SH) * H, 0.0, 0.0), image)
reply(ix)
}
}
}
/**
* 幅のピクセル数分だけアクターにメッセージを送る.
*/
@tailrec
def runLoop(ix: Int, results: Array[Future[Any]], image: BufferedImage): Unit =
ix match {
case SW => results foreach(_ apply)
case _ => {
results(ix) = calc !! ((ix, image))
runLoop(ix + 1, results, image)
}
}
/**
* 実行.
* 計算部分の処理時間計測および画像ファイルの書き出し.
*/
def run = {
val results = new Array[Future[Any]](SW)
val image = new BufferedImage(SW, SH, BufferedImage TYPE_3BYTE_BGR)
val t = (new Benchmark { def run = runLoop(0, results, image) } runBenchmark 1)(0) / 1000.0
println("Time: " + t + " s")
val out = new FileOutputStream("mandelbrot.png")
ImageIO write(image, "png", out)
out close
}
}
object MandelbrotMain {
def main(args: Array[String]) =
new Mandelbrot(4800, 4800, -0.005, -0.005, 0.01, 0.01, 10000, 10.0) run
//new Mandelbrot(640, 480, -1.0, -2.0, 2.6666667, 2.0, 1000, 10.0) run
//new Mandelbrot(800, 800, 0.025185, -1.401565, 0.0005, 0.0005, 10000, 10.0) run
}</span><br /><br />
<strong>mandelbrot.scala (アクターなし; 逐次処理)</strong><br /><br />
<span style="font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;">// mandelbrot.scala without actor by nox, 2011.07.09
import scala.annotation.tailrec
import java.io.FileOutputStream
import java.awt.image.BufferedImage
import java.awt.Color
import javax.imageio.ImageIO
import scala.testing.Benchmark
/**
* マンデルブロ集合.
* SW : 画像の幅.
* SH : 画像の高さ.
* T : 集合の上の位置.
* L : 集合の左の位置.
* W : 集合の幅.
* H : 集合の高さ.
* MAX_LOOP : イテレーションの最大回数.
* TH : 閾値.
*/
final class Mandelbrot(SW: Int, SH: Int, T: Double, L: Double, W: Double, H: Double, MAX_LOOP: Int, TH: Double) {
/**
* イテレーションの回数から色を決定して任意のピクセルにセットする.
*/
def setPixel(x: Int, y: Int, a: Any, image: BufferedImage) =
a match {
case -1 => image setRGB(x, y, 0)
case i: Int => image setRGB(x, y, Color HSBtoRGB(i / 100.0f, 1.0f, 1.0f))
}
/**
* マンデルブロ集合の計算部分.
*/
@tailrec
def calcLoop(i: Int, cr: Double, ci: Double, zr: Double, zi: Double): Int = {
val zr2 = zr * zr
val zi2 = zi * zi
i match {
case MAX_LOOP => i
case _ if zr2 + zi2 >= TH => i
case _ => {
val zrzi = zr * zi
calcLoop(i + 1, cr, ci, zr2 - zi2 + cr, zrzi + zrzi + ci)
}
}
}
def calc(cr: Double, ci: Double): Int =
calcLoop(0, cr, ci, 0.0, 0.0)
/**
* 画像のピクセル数分だけループする.
*/
@tailrec
def runLoop(ix: Int, iy: Int, image: BufferedImage): Unit =
(ix, iy) match {
case (SW, _) => None
case (_, SH) => runLoop(ix + 1, 0, image)
case (_, _) => {
setPixel(ix, iy, calc(L + (ix.toDouble / SW) * W, T + (iy.toDouble / SH) * H), image)
runLoop(ix, iy + 1, image)
}
}
/**
* 実行.
* 計算部分の処理時間計測および画像ファイルの書き出し.
*/
def run = {
val image = new BufferedImage(SW, SH, BufferedImage TYPE_3BYTE_BGR)
val t = (new Benchmark { def run = runLoop(0, 0, image) } runBenchmark 1)(0) / 1000.0
println("Time: " + t + " s")
val out = new FileOutputStream("mandelbrot.png")
ImageIO write(image, "png", out)
out close
}
}
object MandelbrotMain {
def main(args: Array[String]) =
new Mandelbrot(4800, 4800, -0.005, -0.005, 0.01, 0.01, 10000, 10.0) run
//new Mandelbrot(640, 480, -1.0, -2.0, 2.6666667, 2.0, 1000, 10.0) run
//new Mandelbrot(800, 800, 0.025185, -1.401565, 0.0005, 0.0005, 10000, 10.0) run
}</span><br />
</span>noxhttp://www.blogger.com/profile/01142433117618538300noreply@blogger.com0