2009年12月25日金曜日

日本全国コンビニ店舗分布地図: 高解像度インタラクティブ版

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

前回に引き続き今回もProcessingによる日本全国コンビニ店舗分布地図について書いてみる。前回は大きさが固定された地図で拡大・縮小や移動ができなかったので、高解像度インタラクティブ版として、それをできるようにしてみた。

まず、日本地図がPNG画像なのでこれをベクタ画像に変更する。フォーマットはSVGだ。前回と同様に「カビパン男と私」で提供されている日本地図のSVGファイルを若干加工して利用させてもらった。

次にこのSVGデータをProcessingで利用する方法について述べる。「ビジュアライジング・データ」にもSVGの利用方法(processing.candy.*)が書いてあるのだが、実はここに書いてある情報は古くて現在では使用できない。現在では、PShapeを使ってSVGを利用する。

PShape mapShape = loadShape("japan.svg");

こんな感じだ。画面の(x, y)座標に表示するにはdraw関数内で以下のように記述すればいい。

shape(mapShape, x, y);

詳しくはProcessingのリファレンスに書かれているので、下記のソースコードなどを参考にして調べてみて欲しい。

さて、これでSVG画像が利用できるようになったので、今度はこれを使って、マウスの左ボタンのドラッグで平行移動を行えるようにしてみる。やり方としてはmouseDragged関数を使う。移動したピクセル分だけ画像も移動する。mouseDragged関数内で移動分のオフセットを変数に入れるようにし、上述のshape関数でオフセット分だけ移動した位置に描画すればいい。

次に画像の拡大・縮小を行う。マウスの右ボタンのドラッグで操作する。マウスを上に移動すれば拡大、下に移動すれば縮小だ。SVG画像の縮小・拡大はscaleメソッドを使う。例えば、上記の例で元の画像を1.5倍にしたければ以下のようにする。

mapShape.scale(1.5);

scaleメソッドによる拡大・縮小は画像の位置やそのオフセットなどもちゃんと正しい位置になるようにしなければならないのだが、それらについてはソースコードを参考にして欲しい。

上記はSVG画像による日本地図の描画について説明したが、コンビニデータの描画についても位置のオフセットと拡大率・縮小率で行うことができ、Placeクラスのdrawメソッドをちょっといじればいい。

以下が完成版。最初に全部のコンビニを表示してしまうと分布が解りづらいので、今回はセブンイレブンとローソンのみを最初に表示するようにした。コンビニの表示・非表示はコンビニ名をクリックすることで切り替えることができる。

This browser does not have a Java Plug-in.
Get the latest Java Plug-in here.


以下にソースコードを示しておく。

convenience.pde

String convData = "data.tsv.gz"; String japanMap = "japan.svg"; PShape mapShape; int totalCount; Place[] places; int placeCount = 0; final float minX = -0.19834; final float maxX = 0.2425; final float minY = -0.1875; final float maxY = 0.1845; PImage mapImage; PFont font; float offset_x = 0.0; float offset_y = 0.0; int center_x = 0; int center_y = 0; float zoom = 1.0; int mx_start, my_start, my_pos; String[] convName = { "Seven-Eleven", "LAWSON", "FamilyMart", "Circle K", "Sunkus", "Daily YAMAZAKI", "MINISTOP", "am/pm", "POPLAR", "Three F", "COMMUNITY STORE", "Cocostore", "SHOP99", "Seicomart", "SAVE ON" }; color[] convColor = { #ff0000, #007fff, #00ff00, #ffff00, #ff00ff, #00ffff, #7f0000, #0000ff, #007f00, #7f7f00, #7f007f, #007f7f, #7f7f7f, #ff7f00, #7f00ff }; boolean[] convIds = { true, true, false, false, false, false, false, false, false, false, false, false, false, false, false }; int numConv = 15; public void setup() { mapShape = loadShape(japanMap); smooth(); loop(); size(int(mapShape.width), int(mapShape.height), JAVA2D); shapeMode(CORNER); mapShape.disableStyle(); font = loadFont("CourierNewPSMT-14.vlw"); textMode(SCREEN); textFont(font); readData(); } public void draw() { background(0); noStroke(); fill(32); shape(mapShape, int(offset_x * zoom) + center_x, int(offset_y * zoom) + center_y); for (int i = 0; i < placeCount; i++) places[i].draw(); for (int i = 0; i < numConv; i++) { fill(convColor[i]); text(convName[i], 15, 15 * (i + 1)); if (convIds[i]) text("*", 5, 15 * (i + 1)); } } float TX(float x) { float offset = width * (1.0 - zoom) * 0.5; return map(x, minX, maxX, offset, width - offset); } float TY(float y) { float offset = height * (1.0 - zoom) * 0.5; return map(y, minY, maxY, height - offset, offset); } void mousePressed() { mx_start = int(mouseX - offset_x * zoom); my_start = int(mouseY - offset_y * zoom); loop(); for (int i = 0; i < numConv; i++) if (mouseX < convName[i].length() * 8 + 15 && mouseY > 15 * i && mouseY < 15 * (i + 1)) convIds[i] = !convIds[i]; my_pos = mouseY; } void mouseReleased() { noLoop(); } void mouseDragged() { if (mouseButton == LEFT) { offset_x = (mouseX - mx_start) / zoom; offset_y = (mouseY - my_start) / zoom; } else if (mouseButton == RIGHT) { if (mouseY - my_pos > 4) { float zoom_out = 9.8 / 10.0; mapShape.scale(zoom_out); zoom *= zoom_out; center_x = int(width * (1.0 - zoom) * 0.5); center_y = int(height * (1.0 - zoom) * 0.5); my_pos = mouseY; } else if (mouseY - my_pos < -4) { float zoom_in = 10.0 / 9.8; mapShape.scale(zoom_in); zoom *= zoom_in; center_x = int(width * (1.0 - zoom) * 0.5); center_y = int(height * (1.0 - zoom) * 0.5); my_pos = mouseY; } } } void readData() { new Slurper(); } Place parsePlace(String line) { String pieces[] = split(line, '\t'); int id = int(pieces[0]); String name = pieces[1]; float x = float(pieces[2]); float y = float(pieces[3]); return new Place(id, name, x, y); }

Place.pde

class Place { int id; String name; float x, y; public Place(int id, String name, float x, float y) { this.id = id; this.name = name; this.x = x; this.y = y; } public void draw() { if (convIds[this.id]) { int xx = (int)TX(x) + int(offset_x * zoom); int yy = (int)TY(y) + int(offset_y * zoom); set(xx, yy, convColor[this.id]); } } }

Slurper.pde

class Slurper implements Runnable { Slurper() { Thread thread = new Thread(this); thread.start(); } public void run() { try { BufferedReader reader = createReader(convData); String line = reader.readLine(); totalCount = int(line); places = new Place[totalCount]; while ((line = reader.readLine()) != null) { places[placeCount] = parsePlace(line); placeCount++; } } catch (IOException e) { e.printStackTrace(); } noLoop(); } }

続きを読む...

2009年12月14日月曜日

Processing: 日本全国コンビニ店舗分布地図

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

「ビジュアライジング・データ Processingによる情報視覚化手法」を読んでみたのだが、Processingの情報処理とその視覚化がとても興味深いものだった。そこで、日本地図を使って何らかの情報を視覚化したいと考え、それならば全国に存在するコンビニ店舗の分布を調べてみようと思い、Processingで作ってみることにした。

因みに以前にも「Processingで分子動力学計算」や「Processingを使ってWebカメラを監視カメラにする」などのブログ記事を書いているように、プログラミング言語Processingはビジュアル関連で広範に応用でき、そのポテンシャルはとても大きい。

まずは、コンビニの店舗情報が必要なので、gooのコンビニ店舗検索を利用し、ここから住所情報を得ることにした。ただ、一度に検索できる件数は5,000件までなので5,000件以上存在するコンビニに対してはPythonで都道府県ごとに検索してマージした。データはHTMLで書かれているのでPythonのre.findallのパターンマッチを使って一気に必要なデータに変換した。

次に、Google MapsのHTTPリクエスト経由のジオコーディングを使って住所情報を経度・緯度に変換した。この辺もPythonを使えばちょちょいのちょいだ。以下のAPIを利用する。

http://maps.google.com/maps/geo?q=住所情報&output=json&sensor=false&key=APIキー

ただし、変換は1日につき15,000件の制限が掛かっているので4万件以上ある住所情報を取得するには3日かかることになる。もっともIPアドレスによる制限なので別のIPを使えば問題ないようだ。

日本地図上にデータを表示するためには、経度・緯度情報からXY座標に変換する必要がある。日本地図は国土地理院によるとユニバーサル横メルカトル(UTM)図法が使われているらしい。簡単に言うと世界地図でよく使われているメルカトル図法の赤道に合わせている中心線を6°ごとの経度線に合わせて作成する方法らしい。そのための変換関数をPythonで以下のように作成した。本当は縮尺の微調整が入るようなのだが、今回はそこまで厳密にする必要はないのでその辺は省いた。

def UTM(lat, lng): lng0 = 35 lat0_list = ((abs(lat - 129.0), 129.0, -3.0), (abs(lat - 135.0), 135.0, 0.0), (abs(lat - 141.0), 141.0, 3.0)) lat0 = sorted(lat0_list)[0][1] lat0_diff = sorted(lat0_list)[0][2] y = (lng - lng0) / 180.0 * math.pi x = math.log(math.tan(math.pi / 4.0 + (lat - lat0) / 180.0 * math.pi / 2.0)) + math.log(math.tan(math.pi / 4.0 + lat0_diff / 180.0 * math.pi / 2.0)) * 2.0 return x, y

以下に、作成したデータの一部を示しておく。データはタブ区切りで、先頭からコンビニID、店舗名、X座標、Y座標になっている。また、ネットワーク経由で経時的にデータを取得する場合、最初にデータ数を知る必要があるため先頭行に全データ数を記述している。

0 旭川豊岡4条5丁目 0.129276 0.152887 0 旭川旭神2条 0.128996 0.152486 0 愛別町 0.132287 0.155449 0 美深西1条 0.128267 0.165430 0 厚真町 0.119953 0.134899

作成した座標データを表示させるための日本地図だが、「カビパン男と私」というサイトの地図データを利用させてもらった。ここの静止画データをペイントソフトで少しだけ加工してできあがりだ。

さて、必要なデータは全て揃ったので、ここからProcessingを使っての視覚化作業になる。とは言っても、前述の「ビジュアライジング・データ」の6章「散布図」を参考に作成したのでそれほど苦労するところはなかった。詳しいことは書籍を参考にして欲しい。日本地図の画像と座標データのスケーリングについては、日本の東西南北の大ざっぱな位置から大体のスケールを求めてあとは手作業で修正した。

以下が作成したプログラムになる。ProcessingはJavaアプレットにエクスポートする機能があるので簡単にウェブに貼り付けることができる。プログラムの使い方だが、最初は全てのコンビニの店舗が表示されている。左上のコンビニ名をクリックすることで表示・非表示を設定できる。データ数が多いのでどうしても店舗が重なってしまうが、そのあたりは表示・非表示を切り替えて対応して欲しい。

This browser does not have a Java Plug-in.
Get the latest Java Plug-in here.



コンビニデータについて。使用したコンビニを店舗数を含めて以下に挙げておく。全ての店舗を合わせると44,447店舗となる。ただし、今回作成したデータで正しく位置情報を取得できなかったものについては省いて使用したので、実際に地図上に表示されているデータは若干少なくなっている。

セブンイレブン 12,056店舗 ローソン 8,580店舗 ファミリーマート 8,463店舗 サークルK 3,342店舗 サンクス 3,449店舗 デイリーヤマザキ 1,185店舗 ミニストップ 1,980店舗 am/pm 989店舗 ポプラ 706店舗 スリーエフ 680店舗 コミュニティ・ストア 153店舗 ココストア 257店舗 SHOP99 936店舗 セイコーマート 967店舗 セーブオン 540店舗

今回のプログラミングはなかなか楽しいものだった。コンビニ店舗情報を視覚化することで、セブンイレブンは店舗数が一番多いが四国や東北北部には店舗がなく都市型だとか、ローソンは地方にも比較的広がっているのだとか、店舗数が少ないコンビニは地域密着型なのだろうとか、いろいろ見えてくるのが面白い。Processingはもっと認知され使われてもいいのではないだろうか。

追記(2009/12/14):

今回はgooのコンビニ店舗検索に載っているコンビニのみを扱ったのだが、コメントで他のコンビニについても言及されたので、セイコーマートとセーブオンの2つのコンビニについて追加してみた。所在地情報はコンビニまっぷから取得した。

追記(2009/12/25):

地図の平行移動や拡大・縮小ができる高解像度インタラクティブ版を作成してみた。

最後に作成したソースコードを示しておく。

convenience.pde

String convData = "data.tsv.gz"; String japanMap = "japan.png"; int totalCount; Place[] places; int placeCount = 0; float minX = -0.130; float maxX = 0.243; float minY = -0.157; float maxY = 0.186; PImage mapImage; PFont font; String[] convName = { "Seven-Eleven", "Lawson", "Family Mart", "Circle K", "Sunkus", "Daily YAMAZAKI", "MINISTOP", "am/pm", "POPLAR", "Three F", "COMMUNITY STORE", "Cocostore", "SHOP99" }; color[] convColor = { #ff0000, #007fff, #00ff00, #ffff00, #ff00ff, #00ffff, #7f0000, #0000ff, #007f00, #7f7f00, #7f007f, #007f7f, #7f7f7f }; boolean[] convIds = { true, true, true, true, true, true, true, true, true, true, true, true, true }; int numConv = 13; public void setup() { mapImage = loadImage(japanMap); size(mapImage.width, mapImage.height, P3D); font = loadFont("CourierNewPSMT-14.vlw"); textMode(SCREEN); textFont(font); readData(); } public void draw() { image(mapImage, 0, 0); for (int i = 0; i < placeCount; i++) places[i].draw(); for (int i = 0; i < numConv; i++) { fill(convColor[i]); text(convName[i], 15, 16 * (i + 1)); if (convIds[i]) text("*", 5, 16 * (i + 1)); } } float TX(float x) { return map(x, minX, maxX, 0, width); } float TY(float y) { return map(y, minY, maxY, height, 0); } void mousePressed() { for (int i = 0; i < numConv; i++) if (mouseX < convName[i].length() * 8 + 15 && mouseY > 16 * i && mouseY < 16 * (i + 1)) convIds[i] = !convIds[i]; } void readData() { new Slurper(); } Place parsePlace(String line) { String pieces[] = split(line, '\t'); int id = int(pieces[0]); String name = pieces[1]; float x = float(pieces[2]); float y = float(pieces[3]); return new Place(id, name, x, y); }

Place.pde

class Place { int id; String name; float x, y; public Place(int id, String name, float x, float y) { this.id = id; this.name = name; this.x = x; this.y = y; } void draw() { if (convIds[this.id]) { int xx = (int)TX(x); int yy = (int)TY(y); set(xx, yy, convColor[this.id]); } } }

Slurper.pde

class Slurper implements Runnable { Slurper() { Thread thread = new Thread(this); thread.start(); } public void run() { try { BufferedReader reader = createReader(convData); String line = reader.readLine(); totalCount = int(line); places = new Place[totalCount]; while ((line = reader.readLine()) != null) { places[placeCount] = parsePlace(line); placeCount++; } } catch (IOException e) { e.printStackTrace(); } } }

続きを読む...

2009年12月12日土曜日

CUDAで作成した分子動力学計算プログラムを書き直してみた

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

以前に、はじめてのCUDAプログラミングで分子動力学計算というブログ記事を書いたことがある。最近このなんちゃって分子動力学計算(MD)プログラムのソースコードを読み直してみたのだが、かなりひどい。しかも、「CUDAプログラミング」でググってみると、この記事が2番目にくる。こんないい加減なコードを参考にされたら読んだ人にも迷惑が掛かるので修正することにした。それほどCUDAに慣れているわけではないが、前回のコードよりはましだと思う。

それにしてもコードの内容がひどい。意味もなく__syncthreads()が入っているし、複数のスレッドから同じグローバルメモリに書き込みしてるし、レジスタを活用してないし、計算の順序は非効率だし、ダメダメだ。

そこでまず、CUDA Visual Profilerで関数のパフォーマンス測定を行った。因みにこのプロファイラはCUDAプログラミングツールキットに含まれているもので、CUDAでプログラミングを行うには必須だと思う。使い方は簡単で、Windowsであればcudaprof.exeを実行して、FileメニューからNew...を選び、プロジェクトの名前と場所を設定して、実行するCUDAプログラムを指定するだけだ。引数が必要ならそれも指定しておく。デフォルトでは4回実行され解析される。それぞれの関数がどれだけ時間がかかったのかプロットされるので、どこがボトルネックになっているのか一目瞭然だ。

まあ、予想していた通り、Calc関数が全体の99.8%を占めていたので、ここから修正を行った。まず、iとjを使っていたループをjだけにした。iはスレッド数で分割されていたが、ブロックと合わせてそれを消した。次に、jのループ内で不必要にグローバルメモリにアクセスしないようにした。例えばchg[i]というグローバル変数はループ外でローカル変数に入れてそれを使うようにする。あとは、できるだけ除算などの演算を減らすようにした。

次に、系全体を中心に戻す関数があるのだが、中心座標を求める部分はホスト側の関数で実装し、系全体を戻す関数のみをCUDAの関数とした。もともとこの処理は毎回やらなくてもよいものなので、指定したステップ毎に1回行うようにした。今回は100ステップに1回としている。なので、全体の処理時間から見れば無視しても良く、頑張ってすべてGPUで処理する必要はないと考えた。

さらに、全体のエネルギーを求める計算を見直してみた。ここは酷くて、前回のコードではCalc関数ごとに全エネルギーの和を求めるという訳が分からないことをしている。そこで、全エネルギーを求める部分の関数を別に書き、さらに並列リダクションを実装した。並列リダクションについてはOptimizing Parallel Reduction in CUDAが詳しい。今回の実装もこれを参考にさせてもらった。ただ、数千個程度の粒子数であれば、デバイスからのメモリコピーを考えてもホスト側で計算させてしまう方が速いかもしれない。一応1024個でスイッチするようにしてあるが、もっと大きな値で良いと思う。

以下に実際の処理速度の結果を示す。テストデータとして4,139個の粒子数を持つデータを用意して、それを100ステップ計算させた。動作環境としてIntel Core2 6600@2.40GHz / NVIDIA GeForce 8800 GTSを使用した。



以前書いたC++のプログラムは307.1秒かかっている。前回のCUDAプログラムで35.2秒、今回のコードが6.6秒だった。C++と比べると46.2倍、前回のCUDAプログラムと比べても5.3倍高速に動作した。

最後に、今回作成したソースコードを示しておく。コンパイル方法や実行方法は前回の記事を参考にして欲しい。

simple_md_gpu.cu

//////////////////////////////////////////////////////////////// // Simple Molecular Dynamics using GPU by nox, 2009.12.10. // // Perform molecular dynamics for pseudo-atoms without internal // interactions, such as bonds, angles, dihedrals. Only vdW and // Coulomb interactions. #include <iostream> #include <iomanip> #include <fstream> #include <sstream> #include <algorithm> #include <numeric> #include <cmath> #include <cutil_inline.h> using namespace std; const int BLOCK = 256; __device__ float Distance(float* crd, int i, int j) { float dx = crd[i*3] - crd[j*3]; float dy = crd[i*3+1] - crd[j*3+1]; float dz = crd[i*3+2] - crd[j*3+2]; return sqrtf(dx * dx + dy * dy + dz * dz); } __global__ void Center(float* crd, float cx, float cy, float cz, int num_atoms) { int i = blockIdx.x * blockDim.x + threadIdx.x; if (i >= num_atoms) return; crd[i*3] -= cx; crd[i*3+1] -= cy; crd[i*3+2] -= cz; } // Calculate energies and forces // Total energy = vdW + Coulomb // vdW // U = eps * [(Rij / r[i])^12 - 2 * (Rij / r[i])^6] // F = -12 * eps / Rij * [(Rij / r[i])^13 - (Rij / r[i])^7] * r_xyz / r[i] // Coulomb // U = SUM_i>j qiqj / r[i] // F = SUM_j qiqj / r[i]^3 * r_xyz __global__ void Calc(float* crd, float* f, float* ene, float* chg, int num_atoms) { int i = blockIdx.x * blockDim.x + threadIdx.x; if (i >= num_atoms) return; const float eps = 0.2f; const float Rij = 2.5f; float r, rij, r12, r6, r3, r_, q, x; float e0, f0, f1, f2; e0 = f0 = f1 = f2 = 0.0f; float c0 = crd[i*3]; float c1 = crd[i*3+1]; float c2 = crd[i*3+2]; float q0 = chg[i]; for (int j = 0; j < num_atoms; j++) { if (i == j) continue; r = Distance(crd, i, j); r_ = 1.0f / r; q = q0 * chg[j]; rij = Rij * r_; r3 = rij * rij * rij; r6 = r3 * r3; r12 = r6 * r6; x = ((12 * eps) * (r12 - r6) + q * r_) * r_ * r_; e0 += eps * (r12 - 2.0f * r6) + q * r_; f0 += x * (c0 - crd[j*3]); f1 += x * (c1 - crd[j*3+1]); f2 += x * (c2 - crd[j*3+2]); } ene[i] = e0; f[i*3] = f0; f[i*3+1] = f1; f[i*3+2] = f2; } template<unsigned int blockSize> __global__ void Energy(float* ene, float* oene, unsigned int n) { extern __shared__ float sdata[]; unsigned int tid = threadIdx.x; unsigned int i = blockIdx.x * (blockSize * 2) + tid; unsigned int gridSize = blockSize * 2 * gridDim.x; sdata[tid] = 0; while (i < n) { sdata[tid] += ene[i] + ene[i+blockSize]; i += gridSize; } __syncthreads(); if (blockSize >= 512) { if (tid < 256) { sdata[tid] += sdata[tid+256]; } __syncthreads(); } if (blockSize >= 256) { if (tid < 128) { sdata[tid] += sdata[tid+128]; } __syncthreads(); } if (blockSize >= 128) { if (tid < 64) { sdata[tid] += sdata[tid+ 64]; } __syncthreads(); } if (tid < 32) { if (blockSize >= 64) sdata[tid] += sdata[tid+32]; if (blockSize >= 32) sdata[tid] += sdata[tid+16]; if (blockSize >= 16) sdata[tid] += sdata[tid+ 8]; if (blockSize >= 8) sdata[tid] += sdata[tid+ 4]; if (blockSize >= 4) sdata[tid] += sdata[tid+ 2]; if (blockSize >= 2) sdata[tid] += sdata[tid+ 1]; } if (tid == 0) oene[blockIdx.x] = sdata[0]; } __global__ void Move(float* crd, float* f, int num_atoms, float cap_range) { int i = blockIdx.x * blockDim.x + threadIdx.x; if (i >= num_atoms) return; crd[i*3] += f[i*3] * 0.01f; crd[i*3+1] += f[i*3+1] * 0.01f; crd[i*3+2] += f[i*3+2] * 0.01f; float r = crd[i*3] * crd[i*3] + crd[i*3+1] * crd[i*3+1] + crd[i*3+2] * crd[i*3+2]; float dr = r - cap_range * cap_range; if (dr > 0.0f) { f[i*3] = -crd[i*3] / cap_range * dr * 0.01f; f[i*3+1] = -crd[i*3+1] / cap_range * dr * 0.01f; f[i*3+2] = -crd[i*3+2] / cap_range * dr * 0.01f; } } class SimpleMD { private: dim3 grid; dim3 threads; int width; int num_steps; int num_atoms; int num_center; float cap_range; float *h_crd, *h_f, *h_ene, *h_chg; float *d_crd, *d_f, *d_ene, *d_chg, *d_oene; float tene; public: SimpleMD(int n, int nc, char*); ~SimpleMD(); void ReadCrd(char*); void PrintCrd(); void CenterPosition(); unsigned int NextPow2(unsigned int x); void GetBlocksAndThreads(int n, int& blocks, int& threads); void TotalEnergyReduce(int threads, int blocks); void TotalEnergy(); void Run(); }; SimpleMD::SimpleMD(int n, int nc, char* fname) : num_steps(n), num_center(nc) { ReadCrd(fname); h_f = new float[num_atoms * 3]; fill(h_f, h_f + num_atoms * 3, 0.0f); h_ene = new float[NextPow2(num_atoms)]; fill(h_ene, h_ene + NextPow2(num_atoms), 0.0f); width = (num_atoms - 1) / BLOCK + 1; grid.x = width; grid.y = 1; grid.z = 1; threads.x = BLOCK; threads.y = 1; threads.z = 1; cudaMalloc((void**)&d_crd, sizeof(float) * num_atoms * 3); cudaMalloc((void**)&d_f, sizeof(float) * num_atoms * 3); cudaMalloc((void**)&d_oene, sizeof(float) * num_atoms); cudaMalloc((void**)&d_ene, sizeof(float) * NextPow2(num_atoms)); cudaMalloc((void**)&d_chg, sizeof(float) * num_atoms); } SimpleMD::~SimpleMD() { cudaFree(d_chg); cudaFree(d_ene); cudaFree(d_f); cudaFree(d_crd); delete[] h_ene; delete[] h_f; delete[] h_chg; delete[] h_crd; } void SimpleMD::ReadCrd(char* fname) { fstream fs(fname, ios_base::in); string line; stringstream ss; if (!fs.is_open()) { cerr << "File open error: " << fname << endl; exit(1); } getline(fs, line); cout << line << endl; getline(fs, line); ss.str(line); ss >> num_atoms; ss.clear(); getline(fs, line); ss.str(line); ss >> cap_range; ss.clear(); h_crd = new float[num_atoms * 3]; h_chg = new float[num_atoms]; for (int i = 0; i < num_atoms; i++) { getline(fs, line); ss.str(line); ss >> h_crd[i*3] >> h_crd[i*3+1] >> h_crd[i*3+2] >> h_chg[i]; ss.clear(); ss.str(""); } fs.close(); } void SimpleMD::PrintCrd() { cout << endl; cout << num_atoms << endl; cout << cap_range << endl; for (int i = 0; i < num_atoms; i++) { for (int j = 0; j < 3; j++) cout << " " << fixed << setw(10) << setprecision(6) << h_crd[i*3+j]; cout << " " << fixed << setw(8) << setprecision(4) << h_chg[i]; cout << endl; } } void SimpleMD::CenterPosition() { float d[3]; cudaMemcpy(h_crd, d_crd, sizeof(float) * num_atoms * 3, cudaMemcpyDeviceToHost); for (int i = 0; i < num_atoms; i++) for (int j = 0; j < 3; j++) d[j] += h_crd[i*3+j]; for (int i = 0; i < 3; i++) d[i] /= num_atoms; cudaMemcpy(d_crd, h_crd, sizeof(float) * num_atoms * 3, cudaMemcpyHostToDevice); Center<<<grid, threads>>>(d_crd, d[0], d[1], d[2], num_atoms); cudaThreadSynchronize(); } unsigned int SimpleMD::NextPow2(unsigned int x) { --x; x |= x >> 1; x |= x >> 2; x |= x >> 4; x |= x >> 8; x |= x >> 16; return ++x; } void SimpleMD::GetBlocksAndThreads(int n, int& blocks, int& threads) { threads = (n < BLOCK * 2) ? NextPow2((n + 1) / 2) : BLOCK; blocks = (n + (threads * 2 - 1)) / (threads * 2); blocks = min(width, blocks); } void SimpleMD::TotalEnergyReduce(int threads, int blocks) { switch (threads) { case 512: Energy<512><<<blocks, threads, sizeof(float) * threads>>>(d_ene, d_oene, num_atoms); break; case 256: Energy<256><<<blocks, threads, sizeof(float) * threads>>>(d_ene, d_oene, num_atoms); break; case 128: Energy<128><<<blocks, threads, sizeof(float) * threads>>>(d_ene, d_oene, num_atoms); break; case 64: Energy< 64><<<blocks, threads, sizeof(float) * threads>>>(d_ene, d_oene, num_atoms); break; case 32: Energy< 32><<<blocks, threads, sizeof(float) * threads>>>(d_ene, d_oene, num_atoms); break; case 16: Energy< 16><<<blocks, threads, sizeof(float) * threads>>>(d_ene, d_oene, num_atoms); break; case 8: Energy< 8><<<blocks, threads, sizeof(float) * threads>>>(d_ene, d_oene, num_atoms); break; case 4: Energy< 4><<<blocks, threads, sizeof(float) * threads>>>(d_ene, d_oene, num_atoms); break; case 2: Energy< 2><<<blocks, threads, sizeof(float) * threads>>>(d_ene, d_oene, num_atoms); break; case 1: Energy< 1><<<blocks, threads, sizeof(float) * threads>>>(d_ene, d_oene, num_atoms); break; } } void SimpleMD::TotalEnergy() { if (num_atoms > 1024) { int th, bl; GetBlocksAndThreads(num_atoms, bl, th); int s = bl; while (s > 1) { GetBlocksAndThreads(s, bl, th); TotalEnergyReduce(th, bl); s = (s + (th * 2 - 1)) / (th * 2); } cudaThreadSynchronize(); cudaMemcpy(&tene, d_oene, sizeof(float), cudaMemcpyDeviceToHost); tene /= 2.0f; } else { cudaMemcpy(h_ene, d_ene, sizeof(float) * num_atoms, cudaMemcpyDeviceToHost); tene = accumulate(h_ene, h_ene + num_atoms, 0.0f) / 2.0f; } } void SimpleMD::Run() { cudaMemcpy(d_crd, h_crd, sizeof(float) * num_atoms * 3, cudaMemcpyHostToDevice); cudaMemcpy(d_f, h_f, sizeof(float) * num_atoms * 3, cudaMemcpyHostToDevice); cudaMemcpy(d_ene, h_ene, sizeof(float) * NextPow2(num_atoms), cudaMemcpyHostToDevice); cudaMemcpy(d_chg, h_chg, sizeof(float) * num_atoms, cudaMemcpyHostToDevice); for (int i = 0; i < num_steps; i++) { if (i % num_center == 0) CenterPosition(); Calc<<<grid, threads>>>(d_crd, d_f, d_ene, d_chg, num_atoms); cudaThreadSynchronize(); TotalEnergy(); Move<<<grid, threads>>>(d_crd, d_f, num_atoms, cap_range); cudaThreadSynchronize(); cout << "Energy (" << setw(7) << i + 1 << "): "; cout << fixed << setw(15) << setprecision(5) << tene << endl; } cudaMemcpy(h_crd, d_crd, sizeof(float) * num_atoms * 3, cudaMemcpyDeviceToHost); cudaMemcpy(h_f, d_f, sizeof(float) * num_atoms * 3, cudaMemcpyDeviceToHost); cudaMemcpy(h_ene, d_ene, sizeof(float) * num_atoms, cudaMemcpyDeviceToHost); } int main(int argc, char** argv) { if (argc != 3) { cerr << "Usage: " << argv[0] << " input_file number_of_steps" << endl; cerr << "Input: line 1 : title" << endl; cerr << " line 2 : number of atoms" << endl; cerr << " line 3 : radius of droplet" << endl; cerr << " line 4-: x-crd y-crd z-crd charge" << endl; exit(1); } stringstream ss; int nstep; cutilDeviceInit(1, argv); ss.str(argv[2]); ss >> nstep; int ncenter = 100; SimpleMD md(nstep, ncenter, argv[1]); md.PrintCrd(); md.Run(); md.PrintCrd(); return 0; }

続きを読む...

2009年12月9日水曜日

C++: 構造体を格納したSTLコンテナに対してソート・探索・削除などのアルゴリズムを適用する

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

C++に慣れている人にとっては当たり前のことかもしれないけど、あまりC++に親しんでいない場合、構造体を格納したSTLコンテナに対してアルゴリズム<algorithm>を有効に活用していないかもしれない。そこで、構造体を格納したvectorなどのSTLコンテナでソートや探索、削除などのアルゴリズムの利用方法を書いておく。

struct A { int n; int* p; };

上記のような構造体はよく見かける形だと思う。構造体Aに整数型変数のnとポインタ型変数のpがあり、例えばnに配列の要素数、pにその配列を確保したりする。こういった構造体を以下のようにvectorなどのSTLコンテナを使って格納することは多々ある。

vector<A> A_list;

これで構造体Aをコンテナに格納できるわけだ。ところで、STLコンテナを使用する一つの理由として便利なアルゴリズムが利用できることが挙げられるだろう。sortやfindなどの定番のアルゴリズムは、それほどC++を利用していない人でもSTLを使ったことがあるのであれば知っていると思う。これらはとても便利だが、上記の構造体を入れたコンテナに対して以下のようには利用できない。

sort(A_list.begin(), A_list.end());

当たり前のことだが、どのような基準でソートするのかコンパイラには分からないからだ。なので、どのような基準でソートさせるのかを関数オブジェクトを使って教えればいい。因みに関数オブジェクトとは、operator()が定義されているクラスのことだ。

例えば構造体Aのメンバ変数であるnを基準に昇順にソートさせたいとする。その場合は以下のような関数オブジェクトを作って、sortアルゴリズムの第3引数に渡せばいい。

class LessAn { public: bool operator()(const A& a, const A& b) { return a.n < b.n; } };

sort(A_list.begin(), A_list.end(), LessAn());

たったこれだけだ。

他のアルゴリズムに対しても基本的には同じように関数オブジェクトを使えばいい。探索(find_if)や削除(remove_if)、繰り返し処理(for_each)などを使ったアルゴリズムの利用例を以下に挙げておく。

1から30までのランダムな要素数nと0から29までのランダムな要素の配列pを持つの構造体Aを10個作成してvectorコンテナA_listに格納した。それに対し複数のアルゴリズムを適用した出力結果を以下に示している。

初期化した構造体Aのコンテナ(A_list): n : p[] 12 : 17 4 10 29 4 18 18 22 14 5 5 1 28 : 1 11 25 2 27 6 21 24 2 3 22 22 21 26 8 5 17 6 11 18 9 22 17 19 25 24 23 21 3 : 3 3 14 22 : 1 23 28 17 14 22 27 27 19 23 21 19 28 16 5 20 12 18 16 10 2 4 29 : 26 15 20 9 10 20 6 21 3 8 9 23 24 4 6 20 16 26 11 28 24 9 26 13 17 28 8 12 9 12 : 3 5 19 18 24 0 27 26 23 6 11 5 15 : 22 0 9 17 3 27 12 16 0 11 6 5 17 15 4 12 : 22 20 10 21 4 16 10 27 11 7 27 7 18 : 3 3 5 29 19 8 11 18 2 16 26 10 3 8 0 11 12 5 11 : 29 24 17 8 3 25 21 2 20 1 26 構造体Aの配列(A.p)をソートして重複データを削除: n : p[] 9 : 1 4 5 10 14 17 18 22 29 18 : 1 2 3 5 6 8 9 11 17 18 19 21 22 23 24 25 26 27 2 : 3 14 17 : 1 2 4 5 10 12 14 16 17 18 19 20 21 22 23 27 28 18 : 3 4 6 8 9 10 11 12 13 15 16 17 20 21 23 24 26 28 11 : 0 3 5 6 11 18 19 23 24 26 27 13 : 0 3 4 5 6 9 11 12 15 16 17 22 27 9 : 4 7 10 11 16 20 21 22 27 13 : 0 2 3 5 8 10 11 12 16 18 19 26 29 11 : 1 2 3 8 17 20 21 24 25 26 29 構造体Aの配列(A.p)の要素数(A.n)でA_listをソート: n : p[] 2 : 3 14 9 : 1 4 5 10 14 17 18 22 29 9 : 4 7 10 11 16 20 21 22 27 11 : 0 3 5 6 11 18 19 23 24 26 27 11 : 1 2 3 8 17 20 21 24 25 26 29 13 : 0 3 4 5 6 9 11 12 15 16 17 22 27 13 : 0 2 3 5 8 10 11 12 16 18 19 26 29 17 : 1 2 4 5 10 12 14 16 17 18 19 20 21 22 23 27 28 18 : 1 2 3 5 6 8 9 11 17 18 19 21 22 23 24 25 26 27 18 : 3 4 6 8 9 10 11 12 13 15 16 17 20 21 23 24 26 28 構造体Aの配列(A.p)を素数のみにする: n : p[] 1 : 3 3 : 5 17 29 2 : 7 11 5 : 3 5 11 19 23 4 : 2 3 17 29 4 : 3 5 11 17 6 : 2 3 5 11 19 29 5 : 2 5 17 19 23 7 : 2 3 5 11 17 19 23 5 : 3 11 13 17 23 構造体Aの配列要素の和が偶数ならA_listから削除: n : p[] 1 : 3 3 : 5 17 29 5 : 3 5 11 19 23 4 : 2 3 17 29 6 : 2 3 5 11 19 29 5 : 3 11 13 17 23 再び構造体Aの配列(A.p)の要素数(A.n)でA_listをソート: n : p[] 1 : 3 3 : 5 17 29 4 : 2 3 17 29 5 : 3 5 11 19 23 5 : 3 11 13 17 23 6 : 2 3 5 11 19 29 構造体Aの配列(A.p)内で指定した数値と一致する要素を持つ構造体Aを探す: 0 : 指定の数値を要素として持つ構造体Aは見つかりませんでした. 1 : 指定の数値を要素として持つ構造体Aは見つかりませんでした. 2 : 4 : 2 3 17 29 3 : 1 : 3 4 : 指定の数値を要素として持つ構造体Aは見つかりませんでした. 5 : 3 : 5 17 29 6 : 指定の数値を要素として持つ構造体Aは見つかりませんでした. 7 : 指定の数値を要素として持つ構造体Aは見つかりませんでした. 8 : 指定の数値を要素として持つ構造体Aは見つかりませんでした. 9 : 指定の数値を要素として持つ構造体Aは見つかりませんでした.

以下、ソースコード。

#include <iostream> #include <iomanip> #include <vector> #include <algorithm> #include <numeric> #include <cstdlib> using namespace std; const int num_A = 10; // A_listコンテナ内の構造体Aの数. const int max_n = 30; // 構造体A内の配列の最大数. // 構造体A : 要素数n, 要素の配列pを持つ. struct A { int n; int* p; // 乱数で配列の要素を作成. A(int n_, int* p_) : n(n_), p(p_) { for (int i = 0; i < n; i++) p[i] = rand() % max_n; } }; // 構造体Aの表示. void print(const A& a) { cout << setw(3) << a.n << " :"; for (int i = 0; i < a.n; i++) cout << " " << a.p[i]; cout << endl; } // A_listコンテナの表示. void print(const vector<A>& a_list) { cout << " n : p[]" << endl; for (vector<A>::const_iterator p = a_list.begin(); p != a_list.end(); ++p) print(*p); cout << endl; } // A_listコンテナを構造体Aの配列要素数でソートするための関数オブジェクト. class LessAn { public: bool operator()(const A& a, const A& b) { return a.n < b.n; // 昇順. } }; // 構造体Aの配列の重複を削除するための関数オブジェクト. class UniqueAp { public: void operator()(A& a) { sort(a.p, a.p + a.n); // 配列内をソート. a.n = unique(a.p, a.p + a.n) - a.p; // 重複を削除. } }; // 任意の数値を素数判定する関数オブジェクト. class IsPrime { private: static vector<int> primes; int i, j; public: bool operator()(int n) { if (n < 2) return false; if (binary_search(primes.begin(), primes.end(), n)) return true; for (vector<int>::iterator p = primes.begin(); p != primes.end(); ++p) { if (n < (*p) * (*p)) return true; if (n % *p == 0) return false; } for (i = *(primes.end() - 1) + 2; i * i <= n; i += 2) { for (j = 0; i > primes[j] * primes[j] && i % primes[j] != 0; j++); if (i < primes[j] * primes[j]) { primes.push_back(i); if (n % i == 0) return false; } } return true; } }; // 素数を格納する静的メンバコンテナ. static int init_primes[] = { 2, 3, 5, 7 }; vector<int> IsPrime::primes(init_primes, init_primes + 4); // 構造体Aの配列の要素を素数のみにするための関数オブジェクト. class PrimeAp { public: void operator()(A& a) { // 素数のみを区分けする. // 順番を保持したいのでpartitionではなくstable_partitionを使用. a.n = stable_partition(a.p, a.p + a.n, IsPrime()) - a.p; } }; // 構造体Aの配列の要素の和が偶数がどうか判定する関数オブジェクト. class SumAp_even { public: bool operator()(const A& a) { // 配列内の要素の和を2で割った余りが0か? return accumulate(a.p, a.p + a.n, 0) % 2 == 0; } }; // 構造体Aの配列に指定の数値を含むか判定する関数オブジェクト. class FindAp { private: int n; public: FindAp(int n_) : n(n_) { } bool operator()(const A& a) { // 配列から指定の数値をfindで探索. // 因みにアルゴリズムはSTLコンテナだけではなくポインタでも利用できる. return find(a.p, a.p + a.n, n) != a.p + a.n; } }; int main() { vector<A> A_list; // 構造体Aのvectorコンテナ. // A_listを初期化. for (int i = 0; i < num_A; i++) { // 構造体Aの配列の要素数を乱数で決定. int n = rand() % max_n + 1; A_list.push_back(A(n, new int[n])); } cout << "初期化した構造体Aのコンテナ(A_list):" << endl; print(A_list); cout << "構造体Aの配列(A.p)をソートして重複データを削除:" << endl; for_each(A_list.begin(), A_list.end(), UniqueAp()); print(A_list); cout << "構造体Aの配列(A.p)の要素数(A.n)でA_listをソート:" << endl; sort(A_list.begin(), A_list.end(), LessAn()); print(A_list); cout << "構造体Aの配列(A.p)を素数のみにする:" << endl; for_each(A_list.begin(), A_list.end(), PrimeAp()); print(A_list); cout << "構造体Aの配列要素の和が偶数ならA_listから削除:" << endl; A_list.erase(remove_if(A_list.begin(), A_list.end(), SumAp_even()), A_list.end()); print(A_list); cout << "再び構造体Aの配列(A.p)の要素数(A.n)でA_listをソート:" << endl; sort(A_list.begin(), A_list.end(), LessAn()); print(A_list); cout << "構造体Aの配列(A.p)内で指定した数値と一致する要素を持つ構造体Aを探す:" << endl; vector<A>::iterator p; for (int i = 0; i < 10; i++) { p = find_if(A_list.begin(), A_list.end(), FindAp(i)); cout << setw(3) << i << " : "; if (p != A_list.end()) print(*p); else cout << "指定の数値を要素として持つ構造体Aは見つかりませんでした." << endl; } return 0; }

続きを読む...

2009年11月16日月曜日

次世代スーパーコンピュータが必要な理由

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

先日の事業仕分けによる次世代スーパーコンピュータ(次世代スパコン)の予算見直しで事実上の「凍結」との結論が出たことで、多くの人々から次世代スパコンについて注目が集まることになったが、情報不足のためか、一部誤解があるようだ。そこで、自分の知っている範囲で次世代スパコンについて記したいと思う。もし自分の知識が至らず間違っている場合は指摘して頂けると有り難い。

次世代スーパーコンピュータ

次世代スパコンは富士通のCPU、"Vinus" SPARC64 VIIIfxか、その後継CPUで構成する公算が高いが、このCPUはスカラ型だ。もともと、NECと日立がベクトル型のCPUを開発する予定であったが、撤退によりベクトル型とスカラ型の混成システムから、スカラ型のみのシステムに変更された。因みに、この富士通のCPU "Vinus"は現時点で世界最高速のCPUであり、国産で高速なCPUを開発できるのかという疑問も払拭できている。

ところで、この撤退により次世代スーパーコンピュータ開発について危惧する声が出たのだが、実のところ次世代スパコン開発に関わるユーザや開発者の一部では、開発がしやすくなったとして喜んでいたりする。もともとの混成システムでは、実用的に使う場合、ベクタ部はベクタ部のみの利用、スカラ部ではスカラ部のみの利用でしかパフォーマンスが出せず、これでベクタ部とスカラ部を一台に入れる意味があるのかというもっともな疑問が出ていたようだ。また、次世代スーパーコンピュータ開発実施本部のプロジェクトリーダーである渡辺氏がNEC出身であることからのあらぬ疑いもされずに済むようになったと思う。

NECと日立が撤退した理由として、百数十億円の負担が重荷になったと答えている。しかし、実際はそんな負担額では済まなかったというのが関係者間での通説だ。つまり、それ以上の巨額の開発資金を企業は自腹で負担しなくてはならなかったのだ。では何故そのような負担をしてまでも次世代スパコンに参加したかったのか。それは「世界一高速なコンピュータを作りました」という事実が企業にとって非常に大きな宣伝になるから。世間一般ではそれで世界一の技術があると見なしてくれるのだ。それに、売り上げももちろん重要だが、技術者の士気も上がることのメリットが大きい。技術者の士気はなかなか目に見えにくいものだが、そういう技術者の士気があってこそ、今日の技術立国日本に成長したのだと思う。そして、国や理研側としては実際にかかる費用よりかなり割安に開発できるのだからある意味Win-Winの関係という訳だ。

因みにスパコン開発のために必要な総費用は1,154億円(現在ではNECと日立の撤退により若干増えて1,230億円)は7年間の総費用であり、それにこの金額はスパコンで利用するアプリケーション開発や研究、センターの建設などすべてを含めた費用であり、ハードウェア作製の資金はこの中の一部が担っている。これを考えると決して高い予算ではないように思う。

次世代スパコンについてもっと詳しく知りたい場合は、理化学研究所 次世代スーパーコンピュータ開発実地本部が参考になる。

研究者からみたスーパーコンピューティング

スパコンを利用する研究者からしてみれば世界トップレベルの速いコンピュータ(CPUだけではなくメモリやネットワークのレイテンシを含む)がちゃんと利用できればそれでいいのだが、それにはまず国がちゃんと音頭を取ってくれないとそんな素敵なコンピュータは使えないわけで、国がやるなら日本のためにということで、前述の次世代スーパーコンピュータということになる。

で、本当に重要なのは次世代スパコンを使った実際の研究であり、そのためにグランドチャレンジ・アプリケーションの研究開発などのプロジェクトがある。これは、ライフサイエンス分野ナノテクノロジー分野が対象で非常に大きな役割を果たすと思われる。これらに加えて産学の研究者・技術者へのインフラとしての意味合いも強い。共同利用を積極的に推進しているのだ。

実はこれには訳があって、前回の地球シミュレータ運用の際に、アプリケーション開発が難しい、特定分野に偏っている、申請が通りにくいなどの問題点が指摘されたので、よりよく研究に利用してもらうにはどうしたら良いかという視点からこのような方針になったようなのだ。そして、すでに複数の大学や研究機関との連携をとりつつある。

という訳で、今後代替措置もなくスパコンが使えなくなったりすると産学の基礎研究、応用研究に大きなダメージが出ることは間違いないと思う。

予算「凍結」について

ここでは政治が云々、政党が云々というつもりはない。どのようにすればスーパーコンピューティングの重要性を理解してもらえるか考えてみる。

今回の事業仕分けによる次世代スパコン「凍結」の際、以下のような意見が出た。

「世界一を目指す理由は何か。2位ではだめなのか」
「一時的にトップを取る意味はどれくらいあるか」
「一番だから良いわけではない」

実際にこう思っている一般の方々もいると思う。それに対して、前述の理研のページで一般の方にも分かる説明が出ている。

トップの座にいられるのは、わずか数年なんですね。

個々の研究者や科学を理解している人たちも早速ブログやTwitterで多くの意見を出している。このような議論は積極的に行うべきだと思うし、多くの人の目に触れれば触れるほど認知度は上がってくるはずだ。また、今回の事業仕分けについて文科省でも意見を募集しているので、ここで意見を述べるのも効果があるかもしれない。

さらに、一部の議員の方にも次世代スパコンについては理解してもらっているようで、民主党の藤末健三議員のブログTwitter木俣佳丈議員のブログでも言及されている。

以上をまとめると、技術力の維持と科学の発展のために次世代スパコンは必要であり、このような理由を積極的にアピールして一般の方々にそれを知ってもらうことが大事なのだと思う。

続きを読む...

2009年11月1日日曜日

Google Wave: 「てめーはおれを怒らせた」ボットを作ってみた

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

「てめーはおれを怒らせた」ボットといきなり言われても何のことだか分からないと思うけど、要はユーザの発言に含まれる特定の文字に反応して動作するGoogle Waveのボット billowlet@appspot.com を作成したのだ。

追記(2010/3/7): Google Wave Robots API v2がリリースされたので、新しいAPIでボットを書き直してみた。古いAPIは今後使用できなくなるようなので、今のうちに新しいAPIに移行した方が良いと思う。


先日、Google Waveの開発者用アカウントをGoogleからもらえたので早速何か作ってみることにした。どうせならGoogle App Engineを利用したかったので簡単に作れそうなボット(bot)を作成してみたわけだ。因みに元ネタのジョジョの奇妙な冒険とは言い回しが似ているというだけで内容とは関係ないので悪しからず。

このボットはユーザの発言内容に含まれる特定の文字に反応して動作し、ネコみたいな喋り方にさせられたり、英語や中国語で喋らされたり、笑いっぱなしにさせられたり、どもらされたりと、ユーザの発言がいろいろ変化する。因みに、「ごめんなさい」もしくは「ゆるして」と言うことで、動作を取り消すことができる。画像は実際の使用例だ。

ボットのアカウントはbillowlet@appspot.comで、使い方はGoogle WaveのContactsパネル内にある+ボタンでこのアカウントを追加するだけだ。そして、このボットアカウントをWaveに参加させればよい。因みに、ボットに使われている仔猫のアイコンについてはフリー素材ROKOの画像を利用させてもらった。

以下に「てめーはおれを怒らせた」ボットのデモのリンクを張っておくが、こちらはWave Sandboxアカウント所持者のみが利用できる。

「てめーはおれを怒らせた」ボット デモ

実際に使ってみて、このGoogle Waveはとても面白いWebアプリケーションだと思った。今回のアプリケーション作成でも将来性を強く感じることができたので、これから多くの人が使うことでどのように発展していくのか楽しみだ。

以下にGoogle App Engine上のソースコードを示しておく。ところで、stringモジュールのtranslate関数がunicodeに対してちゃんと動作しないのは何とかならないものかな。

追記(2009/11/23):

開発者用アカウント(Sandbox)では正常に動作するボットの発言が、Previewアカウントでは表示されないので調べてみたら、どうやらGAEにおけるGoogle Wave APIのバグのようだ。現在、GAEでCreateBlip/CreateChildを使うと正常に動作しない。今のところ対策はないようなので、修正されるまで待つしかないようだ。因みにそれ以外のAPIは動作するようで、ボットの発言はないが自分の発言はちゃんとボットによって変更されて表示される。

追記(2009/11/24):

Google Wave APIのバグが直って、Previewアカウントでも正常に動作するようになったみたい。しかし、今度はしばらくするとボットが反応しなくなるバグがあるようだ。プレビュー版だし不具合が出るのは仕方がないな。未知のバグを見つけたらできるだけ報告するようにしよう。

追記(2009/11/27):

ボットが反応しなくなるバグが直ったようだ。これでボットの開発については問題なくなったかな。

追記(2010/1/31):

最近、ボットのアイコンや名前が表示されないと思っていたら、ボットのプロフィールが表示されないバグがあるみたい。

billowlet.py

#!/usr/bin/env python # -*- coding: utf-8 -*- import re, string, urllib, urllib2, random from google.appengine.ext import db from waveapi import events from waveapi import model from waveapi import robot from waveapi import document from waveapi import simplejson HIRAGANA = u"あいうえおかきくけこさしすせそたちつてとなにぬねのはひふへほまみむめもやゆよわをん" KATAKANA = u"アイウエオカキクケコサシスセソタチツテトナニヌネノハヒフヘホマミムメモヤユヨワヲン" class UserData(db.Model): user = db.StringProperty(multiline=False) status = db.StringProperty(multiline=False) def OnParticipantsChanged(properties, context): added = properties["participantsAdded"] for participant in added: if participant != "billowlet@appspot.com": Notify(context, participant) def OnRobotAdded(properties, context): root_wavelet = context.GetRootWavelet() root_wavelet.CreateBlip().GetDocument().SetText(u"おれはbillowlet。よろしくな。") def Translate(text, from_lang, to_lang): url = "http://ajax.googleapis.com/ajax/services/language/translate?v=1.0&q=%s&langpair=%s%%7C%s" % (urllib.quote(text.encode("utf-8")), from_lang, to_lang) data = simplejson.loads(urllib2.urlopen(url).read()) return data["responseData"]["translatedText"] def BeCat(text): text = re.sub(u"。", u"にゃ。", text) text = re.sub(u"?|\?", u"にゃぁ?", text) text = re.sub(u"!|\!", u"にゃっ!", text) text = re.sub(u"な", u"にゃ", text) if text[-1] not in (u"。", u"?", u"?", u"!", u"!"): text += u"にゃ" return text def DotMarks(text): from_kana = u"かきくけこさしすせそたちつてとはひふへほカキクケコサシスセソタチツテトハヒフヘホ" to_kana = u"がぎぐげござじずぜぞだぢづでどばびぶべぼガギグゲゴザジズゼゾダヂヅデドバビブベボ" for f, t in zip(from_kana, to_kana): text = text.replace(f, t) return text def Laugh(text, ch): for i in range(len(text), 0, -1): text = text[:i] + ch + text[i:] return text.replace(u"。", u"").replace(u"、", u"") + ch * random.randint(3, 10) def Kana(idx): return HIRAGANA[idx] + u"|" + KATAKANA[idx] def OnBlipCreated(properties, context): blip = context.GetBlipById(properties["blipId"]) doc = blip.GetDocument() user = blip.GetCreator() data = UserData.all().filter("user =", user).fetch(1) if data: status = data[0].status else: status = None text = doc.GetText().strip() if status and re.match(u"ごめんなさい|ゆるして", text): blip.CreateChild().GetDocument().SetText(u"ゆるしてやる。やれやれだぜ。") if data: data[0].delete() angry = range(len(HIRAGANA)) random.shuffle(angry) if status: if status == u"ネコになれ。": text = BeCat(text) doc.SetText(text) elif status == u"英国紳士になれ。": text = Translate(text, "ja", "en") doc.SetText(text) elif status == u"台湾に住め。": text = Translate(text, "ja", "zh-TW") doc.SetText(text) elif status == u"インチキ日本人になれ。": text = Translate(Translate(text, "ja", "en"), "en", "ja") doc.SetText(text) elif status == u"大事なことは2回言え。": doc.SetText(text + text + u"\n大事なことなので2回言いました。") elif status == u"笑え。": doc.SetText(Laugh(text, u"w")) elif status == u"濁れ。": doc.SetText(DotMarks(text)) elif status == u"逆さまになれ。": doc.SetText(text[::-1]) elif status == u"どもれ。": doc.SetText((text[0] + u"、") * 3 + text) elif status == u"うざくなれ。": doc.SetText(Laugh(text, u"っ")) else: userdata = UserData() if re.search(Kana(angry[0]), text): status = u"ネコになれ。" elif re.search(Kana(angry[1]), text): status = u"英国紳士になれ。" elif re.search(Kana(angry[2]), text): status = u"台湾に住め。" elif re.search(Kana(angry[3]), text): status = u"インチキ日本人になれ。" elif re.search(Kana(angry[4]), text): status = u"大事なことは2回言え。" elif re.search(Kana(angry[5]), text): status = u"笑え。" elif re.search(Kana(angry[6]), text): status = u"濁れ。" elif re.search(Kana(angry[7]), text): status = u"逆さまになれ。" elif re.search(Kana(angry[8]), text): status = u"どもれ。" elif re.search(Kana(angry[9]), text): status = u"うざくなれ。" if status: blip.CreateChild().GetDocument().SetText(u"てめーはおれを怒らせた。%s" % status) userdata.user = user userdata.status = status userdata.put() def Notify(context, participant): root_wavelet = context.GetRootWavelet() root_wavelet.CreateBlip().GetDocument().SetText(u"%s、おれを怒らせるなよ。" % participant.split("@")[0]) def main(): myRobot = robot.Robot("billowlet", image_url="http://billowlet.appspot.com/assets/icon.gif", version="1", profile_url="http://billowlet.appspot.com/") myRobot.RegisterHandler(events.WAVELET_PARTICIPANTS_CHANGED, OnParticipantsChanged) myRobot.RegisterHandler(events.BLIP_SUBMITTED, OnBlipCreated) myRobot.RegisterHandler(events.WAVELET_SELF_ADDED, OnRobotAdded) myRobot.Run() if __name__ == "__main__": main()

app.yaml

application: billowlet version: 1 runtime: python api_version: 1 handlers: - url: /_wave/.* script: billowlet.py - url: /assets static_dir: assets

続きを読む...

2009年10月26日月曜日

Androidアプリケーションを開発するための準備

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

来月、海外に出張することもあってdocomoのHT-03Aを購入してしまった。本当は海外で使うというのは口実でAndroidアプリケーションを開発したいだけだったりする。ここでは個人的なメモも兼ねて、Windows上でAndroidアプリケーションを開発するための準備について記しておく。

まず、Eclipse Classicをダウンロード(現時点の最新版は3.5.1)して、それを適当なディレクトリに展開する。C:\に展開した場合、C:\eclipse\というディレクトリができる。因みに日本語版でも問題ないと思うが、自分は英語版のEclipseを使っている。

次にAndroid SDKをダウンロード(現時点の最新プラットフォームは2.0 "eclair")して、やはり適当な場所に展開する。C:\android\に展開した場合、C:\android\android-sdk-windows\というディレクトリが作成される。C:\android\android-sdk-windows\toolsを環境変数PATHに追加しておく。

展開したディレクトリにあるSDK Setup.exeを実行すると、Android SDK and AVD Managerが起動するので、そのままプラットフォームをアップデートする。もし、SSLでエラーが出る場合は、Settingsの"Force https://... sources to be fetched using http://..."のチェックボックスをオンにすること。Installed Packagesでインストールされたプラットフォームを確認できる。必要なバージョンのプラットフォームがインストールされていないときは、Available Packagesで必要なバージョンをインストールできる。

Androidのアップデートが終わったら、Eclipseを起動させる。最初に作業ディレクトリ(workspace)を訊かれるが適当なディレクトリを指定する。起動したら[Help]-[Install New Software...]でAddボタンを押し、NameにAndroid Plugin、Locationにhttps://dl-ssl.google.com/android/eclipse/を指定する。Developer Toolsという項目が表示されるので、チェックボックスをチェックしてNextボタンを押し、最後にFinishボタンを押して、Eclipseを再起動する。

次いで、[Window]-[Preferences]でダイアログを開き、ツリーメニューのAndroidを選択する。SDK LocationにC:\android\android-sdk-windowsを指定して、Applyボタンを押す。

Eclipse上からAndroidのエミュレータを利用するために、Android仮想デバイス(Android Virtual Device, AVD)を作成する。[Windows]-[Android SDK and AVD manager]で作成するか、コマンドラインで以下のコマンドを実行する(ここではデバイス名をmy_avdとしているが、どのような名前でも良い)。自分はコマンドラインで作成した。因みにtargetオプションで使用するプラットフォームのバージョンを設定できる。

android create avd --target 3 --name my_avd

あとは、コードを書くだけだ。動作の確認の意味でも、最初はHello, Worldアプリケーションあたりを作成するのが良いと思う。

Androidの実機でデバッグするには、まず、開発しているプロジェクトのManifestファイルのApplicationタグにandroid:debuggable="true"を記述する(または、ダイアログのApplicationタグの項目でDebuggableをtrueにする)。

さらに、実機(自分の場合はHT-03A)のメニューで、[設定]-[アプリケーション]-[開発]-[USBデバッグ]をチェックする。開発中にスリープモードにしたくない場合は、[スリープモードにしない]をチェックしておく。そして、実機と開発するPCをUSBで接続する。そうするとWindowsが新しいハードウェアを見つけるのでWindows Updateはせずに特定の場所からインストールする。場所はC:\android\android-sdk-windows\usb_driverを指定する。最後にWindowsを再起動する。

実機をUSB接続してEclipseを起動し、先ほどandroid:debuggable="true"にしたプロジェクトを開き、[Run]-[Run Configurations]を選択して、ダイアログを開く。ツリーメニューのAndroid Applicationから現在開いているプロジェクトを選択して、右のペインのTargetタグをクリックし、Deployment Target Selection ModeをManualにする。これでアプリケーション実行時にAndroid Device Chooserダイアログが開くので、実機を選択して実行する。実機上でアプリケーションが起動したら成功だ。

因みに作成したアプリケーションを配布するには、Java SEのkeytoolでキーストアと鍵を作成して、EclipseのAndroid Toolsからアプリケーションに署名をする必要がある。Android Marketを利用したい場合は、デベロッパープロフィールを作成して、クレジットカードで$25を支払えば配布できるようになる。

追記(2009/10/28): Android 2.0 (eclair)が正式リリースされたので、それに合わせて内容を修正した。

続きを読む...

2009年10月21日水曜日

無線LANで現在位置を取得してGoogleマップとストリートビューで表示する

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

昨年の8月に書いた記事、「Google App Engine: Gearsで位置情報を取得してストリートビューで表示させる」で、無線LANから取得した現在位置をGoogleマップストリートビューで表示するウェブアプリを作ったが、今回、それをブログなどにもそのまま貼り付けられるようにしてみた。もともと、JavaScriptを使ったアプリなのでそれほど変更するところもなかったのだが、無用な部分を削除し、きちんと画面にエラーメッセージを表示するようにした。動作させるためにはGearsをインストールしている必要がある。

無線LANから現在位置情報を取り出すためにGearsのGeolocation APIを使った。位置情報を取り出したらあとはGoogle Maps APIでそれをGoogleマップとストリートビューで表示させるだけだ。因みに、今年の7月からGoogleマップの画面で左上にある四角いボタンを押すと現在位置を表示できるようになっている。

今回のアプリでは、現在位置を示すマーカーに自作の扇形アイコンを使った。ストリートビューの向きと連動している。ペグマン(人型アイコン)でもよかったのだろうけど、個人的には向いている方向が分かりづらいと思う。



以下にソースファイルを示す。ページを読み込んだ際にGearsにより位置情報を取得している。

<script src="http://maps.google.com/maps?file=api&v=2.x&key=Google Maps API Key" type="text/javascript"></script> <script type="text/javascript" src="http://code.google.com/apis/gears/gears_init.js"></script> <script type="text/javascript"> var map; var pano; var overlayInstance = null; var marker; var info_window = false; var lat, lng; var error_msg = ""; function updatePosition(position) { lat = position.latitude; lng = position.longitude; initialize(); } function initialize_geo() { if (window.google && google.gears) { var geo = google.gears.factory.create("beta.geolocation"); geo.getCurrentPosition(updatePosition, function(positionError) { document.getElementById("error").innerHTML = positionError.message; }); } else { error_msg = "位置情報の取得に<a href=\"http://gears.google.com/\">Gears</a>が必要です。"; document.getElementById("error").innerHTML = error_msg; // 秋葉原. lat = 35.698584; lng = 139.774216; initialize(); } } function initialize() { if (GBrowserIsCompatible()) { var latlng = new GLatLng(lat, lng); var zoom = 15; marker = new GMarker(latlng, {icon: getArrowIcon(0.0), clickable: false}); pano = new GStreetviewPanorama(document.getElementById("pano")); GEvent.addListener(pano, "error", handleNoFlash); GEvent.addListener(pano, "initialized", moveStreet); GEvent.addListener(pano, "yawchanged", yawStreet); map = new GMap2(document.getElementById("map_canvas")); GEvent.addListener(map, "click", moveWalker); map.addControl(new GSmallMapControl()); map.addControl(new GMapTypeControl()); map.setCenter(latlng, zoom); map.addOverlay(marker); map.enableScrollWheelZoom(); pano.setLocationAndPOV(latlng); toggleOverlay(); GEvent.addListener(map, "click", function(overlay, latlng) { pano.setLocationAndPOV(latlng); }); } var client = new GStreetviewClient(); client.getNearestPanorama(latlng, setNewLatLng) } function setNewLatLng(data){ if (data.code != 200) { return; } var new_latlng = data.location.latlng; pano.setLocationAndPOV(new_latlng); map.panTo(new_latlng); marker.setLatLng(new_latlng); } function getArrowIcon(bearing) { var icon = new GIcon(); icon.image = getArrowUrl(bearing); icon.iconSize = new GSize(92, 92); icon.iconAnchor = new GPoint(46, 46); return icon; } function getArrowUrl(bearing) { var id = 3 * Math.round(bearing / 3); return "view_marker_" + id + ".png"; } function moveStreet(location_) { map.setCenter(location_.latlng); marker.setLatLng(location_.latlng); document.getElementById("error").innerHTML = error_msg; } function yawStreet(yaw_) { marker.setImage(getArrowUrl(yaw_)); } function moveWalker(overlay_, latlng_) { var client = new GStreetviewClient(); client.getNearestPanorama(latlng_, setNewLatLng) } function toggleOverlay() { if (!overlayInstance) { overlayInstance = new GStreetviewOverlay(); map.addOverlay(overlayInstance); } else { map.removeOverlay(overlayInstance); overlayInstance = null; } } function toggleMarkerLocation() { if (!info_window) { pov = pano.getPOV(); var displayString = [ marker.getLatLng(), "POV yaw: " + pov.yaw, "POV pitch: " + pov.pitch, "POV zoom: " + pov.zoom ].join("<br />"); map.openInfoWindowHtml(marker.getLatLng(), displayString); info_window = true; } else { map.closeInfoWindow(); info_window = false; } } function handleNoFlash(errorCode) { if (errorCode == GStreetviewPanorama.ErrorValues.FLASH_UNAVAILABLE) { document.getElementById("error").innerHTML = "Flashがサポートされていません。"; } else if (errorCode == GStreetviewPanorama.ErrorValues.NO_NEARBY_PANO) { document.getElementById("error").innerHTML = "ストリートビューの範囲外です。"; } } if (window.attachEvent){ window.attachEvent("onload", initialize_geo); window.attachEvent("onunload", GUnload); } else { window.addEventListener("load", initialize_geo, false); window.addEventListener("unload", GUnload, false); } </script> <table><tr> <td><div id="pano" style="float:left; background-color: black; width: 250px; height: 250px"></div></td> <td><div id="map_canvas" style="width: 250px; height: 250px"></div></td> </tr></table> <div><input type="button" onclick="toggleOverlay()" value="ストリートビューレイヤー" /> <input type="button" onclick="toggleMarkerLocation()" value="現在位置" /> <span id="error"></span></div>

続きを読む...

2009年10月5日月曜日

Python: 画像で与えられた迷路に対し2点間の最短経路を求める

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

迷路の描かれた画像に対して、ピクセルの座標で指定したスタート地点とゴール地点の最短経路を求めるプログラムをPython+PILで書いてみた。使用する画像は、デジカメで撮ったものでも、ウェブから拾ってきたものでも、ペイントソフトで自作したものでも構わない。

まずは使用例を見て欲しい。この画像は携帯カメラで撮った自作の簡単な迷路だ(画像上)。それに対して指定した2点間の最短経路を赤線で示してみた(画像下)。ピクセル単位で計測しているので赤線が若干ガタガタしていて完全な最短経路ではないがほぼ最短と考えていいだろう。迷路画像(画像上)をmaze01.jpgとし、スタート地点の座標が(240, 160)、ゴール地点の座標が(210, 400)の場合、コマンドラインで以下のように実行する。

maze_solver.py maze01.jpg -s 240 160 -g 210 400

これで最短経路を求めることができ、画像ビューアが立ち上がって経路の描かれた画像が表示される(画像下)。画像ビューアではなく画像ファイル(ここではmaze01out.jpgとする)に出力したい場合は、以下のようにする。

maze_solver.py maze01.jpg -s 240 160 -g 210 400 -o maze01out.jpg

このプログラムは2つのパートで成り立っていて、一つは画像認識を使った画像データの領域分けであり、もう一つはその領域内の2点間の経路探索である。それぞれ、連結成分ラベル付け(connected component labeling, CCL)アルゴリズムA*探索アルゴリズムを利用して処理している。詳細についてはそれぞれのリンク先と最後に示したソースコードを参照して欲しい。

次の例は、Wikipediaの迷路の項目で例示されている迷路画像(画像上)を解いてみたものだ(画像中、下)。いずれも迷路の真ん中の小部屋(270, 130)をスタート地点とし、画像中段ではゴールを左上(0, 0)、下段では右下(490, 268)とした。スタートとゴールまでの経路は複数あるが、いずれも最短経路を示している。

画像上段をmaze02.pngとした場合、画像中段、下段は、それぞれ以下のように実行して作成した。

maze_solver.py maze02.png -s 270 130 -g 0 0

maze_solver.py maze02.png -s 270 130 -g 490 268

使い方の詳細は、-hオプションを付けて実行すれば良い。

-aオプションで、探索する画像の大きさを指定している。あまり大きな画像だと時間が掛かるのでデフォルトで60,000ピクセルにしている。ただし、複雑な迷路や写りの悪い写真だと画像がつぶれてしまい、正しい結果にならないかもしれないので、そのようなときは指定するピクセルを大きくする必要がある。元の画像ファイルが指定サイズ以下の場合は実サイズのままで処理される。

-tオプションでは、CCLアルゴリズムで使用する色の閾値を指定している。デジカメなどで撮影した画像は、例えば白い部分でも黄色っぽい白や赤っぽい白など完全に同じ色情報でない場合が多く、大体白に近ければ同じ色と見なすわけだが、その許容範囲がこの閾値となる。因みに閾値を0にすると完全に同じ色情報を持っていないと違う色だと見なされる。

Usage: maze_solver.py [options] imagefile Options: -h, --help show this help message and exit -o OUT_IMAGE, --output-image=OUT_IMAGE output-imagefile -s START, --start=START start position [default: (0, 0)] -g GOAL, --goal=GOAL goal position [default: (0, 0)] -a AREA, --area-size=AREA area size [default: 60000] -t THRESHOLD, --threshold=THRESHOLD threshold of the CCL [default: 30]

以下、ソースコード。

maze_solver.py

#!/usr/bin/env python # -*- coding: utf-8 -*- import sys, os, math, heapq, Image, ImageDraw # Connected Component Labeling - Label Equivalence method class CCL: def __init__(self, imagefile, threshold, area): im = Image.open(imagefile) im = im.convert("RGB") self.orig_size = im.size orig_area = self.orig_size[0] * self.orig_size[1] if orig_area <= area: self.size_rate = 1.0 else: self.size_rate = math.sqrt(float(area) / orig_area) im = im.resize((int(self.orig_size[0] * self.size_rate), int(self.orig_size[1] * self.size_rate))) self.size = im.size self.th = threshold self.W = im.size[0] self.D = im.getdata() self.N = len(self.D) self.L = [] self.R = [] for i in range(len(self.D)): self.L.append(i) self.R.append(i) def diff(self, d1, d2): return abs(d1[0] - d2[0]) + abs(d1[1] - d2[1]) + abs(d1[2] - d2[2]) def scanning(self): has_label = False for i in range(self.N): label = self.N if i - self.W >= 0 and self.diff(self.D[i], self.D[i-self.W]) <= self.th: label = min(label, self.L[i-self.W]) if i + self.W < self.N and self.diff(self.D[i], self.D[i+self.W]) <= self.th: label = min(label, self.L[i+self.W]) if i % self.W != 0 and self.diff(self.D[i], self.D[i-1]) <= self.th: label = min(label, self.L[i-1]) if (i + 1) % self.W != 0 and self.diff(self.D[i], self.D[i+1]) <= self.th: label = min(label, self.L[i+1]) if label < self.L[i]: self.R[self.L[i]] = label has_label = True return has_label def analysis(self): for i in range(self.N): label = self.L[i] if label == i: ref = label label = self.R[ref] while True: ref = label label = self.R[ref] if ref == label: break self.R[i] = label def labeling(self): for i in range(self.N): self.L[i] = self.R[self.R[self.L[i]]] def calculate(self): while True: if self.scanning(): self.analysis() self.labeling() else: return self.L, self.W, self.size_rate, self.orig_size, self.size # A* Search Algorithm class AStar: def __init__(self, data, start, goal): self.data = data self.nrow, self.ncol = len(data), len(data[0]) self.start = start self.goal = goal if data[self.start[0]][self.start[1]] != data[self.goal[0]][self.goal[1]]: print "No route!" sys.exit(0) self.passage = data[self.start[0]][self.start[1]] def heuristic(self, pos): return math.sqrt((pos[0] - self.goal[0])**2 + (pos[1] - self.goal[1])**2) def distance(self, path): if len(path) == 0: return 0 sum = 0 p1 = self.start for p2 in path: sum += math.sqrt((p1[0] - p2[0])**2 + (p1[1] - p2[1])**2) p1 = p2 return int(sum) def neighborhood(self, pos): row, col = pos neighbors = [] if row > 0 and self.data[row-1][col] == self.passage: neighbors.append((row - 1, col)) if row < self.nrow - 1 and self.data[row+1][col] == self.passage: neighbors.append((row + 1, col)) if col > 0 and self.data[row][col-1] == self.passage: neighbors.append((row, col - 1)) if col < self.ncol - 1 and self.data[row][col+1] == self.passage: neighbors.append((row, col + 1)) if row > 0 and col > 0 and self.data[row-1][col-1] == self.passage: neighbors.append((row - 1, col - 1)) if row < self.nrow - 1 and col > 0 and self.data[row+1][col-1] == self.passage: neighbors.append((row + 1, col - 1)) if row > 0 and col < self.ncol - 1 and self.data[row-1][col+1] == self.passage: neighbors.append((row - 1, col + 1)) if row < self.nrow - 1 and col < self.ncol - 1 and self.data[row+1][col+1] == self.passage: neighbors.append((row + 1, col + 1)) return neighbors def search(self): path = [] queue = [] checked = [self.start] heapq.heappush(queue, (self.distance(checked) + self.heuristic(self.start), checked)) while len(queue) > 0: score, path = heapq.heappop(queue) last = path[-1] if last == self.goal: return path neighbors = self.neighborhood(last) for nb in neighbors: if nb in checked: continue checked.append(nb) newpath = path + [nb] heapq.heappush(queue, (self.distance(newpath) + self.heuristic(nb), newpath)) return [] def draw_path(in_image, out_image, path, SR): im = Image.open(in_image) im = im.convert("RGB") draw = ImageDraw.Draw(im) draw.line([(int(p[1] / SR), int(p[0] / SR)) for p in path], fill=255) if out_image: im.save(out_image) else: im.show() def main(args): from optparse import OptionParser usage = "usage: %prog [options] imagefile" parser = OptionParser(usage=usage) parser.add_option("-o", "--output-image", type="string", help="output-imagefile", dest="out_image", default="") parser.add_option("-s", "--start", type="int", nargs=2, help="start position [default: %default]", dest="start", default=(0, 0)) parser.add_option("-g", "--goal", type="int", nargs=2, help="goal position [default: %default]", dest="goal", default=(0, 0)) parser.add_option("-a", "--area-size", type="int", help="area size [default: %default]", dest="area", default=60000) parser.add_option("-t", "--threshold", type="int", help="threshold of the CCL [default: %default]", dest="threshold", default=30) options, args = parser.parse_args() if len(args) == 0: parser.print_help() parser.exit() in_image = args[0] print "Connected component labeling phase..." ccl = CCL(in_image, options.threshold, options.area) D, W, SR, ORIG_SZ, SZ = ccl.calculate() print " Width: %d (%d), Height: %d (%d)" % (ORIG_SZ[0], SZ[0], ORIG_SZ[1], SZ[1]) data = [D[i:i+W] for i in range(0, len(D), W)] print "A* searching phase..." print " Start: (%d, %d), Goal: (%d, %d)" % (options.start + options.goal) start = (int(options.start[1] * SR), int(options.start[0] * SR)) goal = (int(options.goal[1] * SR), int(options.goal[0] * SR)) astar = AStar(data, start, goal) path = astar.search() print "Drawing path on the imagefile..." draw_path(in_image, options.out_image, path, SR) if __name__ == "__main__": main(sys.argv)

続きを読む...

2009年9月15日火曜日

C++: STL algorithmのlower_boundとupper_boundの使い方

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

C++ STLのalgorithmに入っているlower_boundとupper_boundが間違えやすいのでメモ。

vector<int>で定義されたコンテナvの要素のうち、A以上、B以下の要素を求める場合。

vector<int>::iterator a = lower_bound(v.begin(), v.end(), A); vector<int>::iterator b = upper_bound(v.begin(), v.end(), B); if (a != v.end() && b != v.begin() && *a <= *(b - 1)) cout << "[" << *a << ", " << *(b - 1) << "]" << endl; else cout << "no element" << endl;

*aが最小の要素、*(b-1)が最大の要素になる。

A以上、B未満の場合。

vector<int>::iterator a = lower_bound(v.begin(), v.end(), A); vector<int>::iterator b = lower_bound(v.begin(), v.end(), B); if (a != v.end() && b != v.begin() && *a <= *(b - 1)) cout << "[" << *a << ", " << *(b - 1) << "]" << endl; else cout << "no element" << endl;

因みに、lower_boundとupper_boundに渡すコンテナはソートされている必要がある。

#include <iostream> #include <vector> #include <algorithm> using namespace std; int main() { vector<int> v; for (int i = 100; i > 0; i--) v.push_back(i); int A = 30; // 下限. int B = 50; // 上限. vector<int>::iterator a, b; cout << "要素数" << v.size() << "で、"; cout << "最小要素" << *min_element(v.begin(), v.end()) << "、"; cout << "最大要素" << *max_element(v.begin(), v.end()); cout << "のコンテナから、" << endl; // lower_bound, upper_boundに渡すコンテナはソートされている必要がある. sort(v.begin(), v.end()); cout << A << "以上、" << B << "以下の要素を求める: "; a = lower_bound(v.begin(), v.end(), A); b = upper_bound(v.begin(), v.end(), B); if (a != v.end() && b != v.begin() && *a <= *(b - 1)) cout << "最小要素 " << *a << ", 最大要素 " << *(b - 1) << endl; else cout << "要素なし" << endl; cout << A << "以上、" << B << "未満の要素を求める: "; a = lower_bound(v.begin(), v.end(), A); b = lower_bound(v.begin(), v.end(), B); if (a != v.end() && b != v.begin() && *a <= *(b - 1)) cout << "最小要素 " << *a << ", 最大要素 " << *(b - 1) << endl; else cout << "要素なし" << endl; return 0; }

出力結果:

要素数100で、最小要素1、最大要素100のコンテナから、 30以上、50以下の要素を求める: 最小要素 30, 最大要素 50 30以上、50未満の要素を求める: 最小要素 30, 最大要素 49

続きを読む...

2009年9月8日火曜日

JavaScriptで素数を計算して動的に表示する

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

JavaScriptを使って素数を計算して、それを動的に表示してみた。今更という感じもするけど、10行そこそこでこのような動的な表示ができるというのは、とても便利だと思う。とは言え、普段はあまりJavaScriptを使ったHTMLは書いてないな。

までの素数


テキストボックスに数字を入力すると、その数までの素数を計算して表示する。最大入力桁を5桁にしたので、99,999までの素数を求めることができる。括弧内の数値は求めた素数の数を示す。

複数のブラウザ(Google Chrome, Firefox, IE8, Safari4)で確認してみたけど、最近はどのブラウザもJavaScriptが速い。今後、ますますJavaScript、特にAjaxなどは要の技術になっていくのだろうな。

以下にソースを示す。

<script type="text/javascript"><!-- function calcPrimes(data) { var primes = []; if (data.value >= 2) primes.push(2); for (i = 3; i <= data.value; i += 2) { for (j = 0; i > primes[j] * primes[j] && i % primes[j] != 0; j++); if (i < primes[j] * primes[j]) primes.push(i); } document.getElementById("primes").innerHTML = primes.length ? " (" + primes.length + "): " + primes.join(", ") : ""; } // --></script> <form><input type="text" size="5" maxlength="5" onkeyup="calcPrimes(this)">までの素数<span id="primes"></span></form>

続きを読む...

2009年9月1日火曜日

結局Pythonを使ってコマンドラインで動作するTwitterクライアントを作ってしまった

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

Twitterのアカウントを取ってから既に1年半になるが、活発に使っているとは言い難い。その原因の一つとしてTwitterのクライアントがある。どうにも自分が利用するのにピッタリだと思うクライアントが見つからなかったのだ。そこで結局、自分の好みに合わせてコマンドライン上で動作するシンプルなTwitterクライアントをPythonで作ってしまった。しかも、ワンライナー(1行プログラム)。

最初の頃はいくつかのクライアントを使ってみたのだが、PCでの作業はシェルで行うことが多いので別のウィンドウを開きたくなかったり、Windows、Unix、MacなどのOSが変わっても同じクライアントを使いたかったり、GUIじゃなくてCUIで操作したかったり、それほど使い込むつもりがないので極力シンプルでコンパクトになっていて欲しかったり、そもそもクライアントをインストールしたくなかったりと、かなり条件を厳しく求めていたら使えるクライアントがなくなってしまい、結局、公式サイトもしくは自分で作成した掲示板からたまにつぶやくだけになってしまった(ただし、Twitterfeedは利用している)。

しかし、最近になってタイムラインをよく眺めるようになり、これが結構面白いと気がついた。自分と違う考え方に触れるのは楽しい。そして、自分ももう少しつぶやいてみようかと思ったのだが、公式サイトや掲示板からのポストはいちいちブラウザを開かなくてはならず手軽だとは言い難かった。

そこで、先日、PythonのワンライナーでTwitterを使えるようにしたこともあり、もう少し改良してそれをクライアントにすることにした。Pythonスクリプトなら、Pythonが入っていればどのOSでも動作するし、使い方はどこでも一緒で、コマンドラインでそのまま使え、タイムラインをgrepなどで簡単に選択表示できる。それに、たった1行のソースコードなので簡単に中身を確認でき、パスワード漏洩やキーロガーなどを心配しなくてもいい。

tw.py

でタイムラインを取得することができ、

tw.py つぶやき tw.py つぶやき http://handasse.blogspot.com/ ブログのURLです。 tw.py "つぶやき&つぶやき" tw.py "つぶやき つぶやき"

などでTwitterにつぶやくことができる。つぶやきに&などの特殊文字や連続でスペースが入る場合などはダブルクォーテーション(")で括ること。つぶやきに入っているURLは自動的にbit.lyによる短縮URLとなる。

また、Windows PowerShellで使用する場合は、プロファイル($profileで表示されるファイルで、Microsoft.PowerShell_profile.ps1などを指し、PowerShell起動時に読み込まれる)に以下の関数を定義しておくと便利だ。

function tw { python 実行パス\tw.py "$args" }

Unixなどで日本語コードが異なっていても、ソースを変更せずにnkfなどを使ったエイリアスで簡単に対処できる。ただ、特殊文字などを使ったときに問題が出るかもしれない。以下の例は、Unix側の表示がUTF-8でシェルがtcshの場合。

alias tw 'tw.py `echo -n \!* | nkf -s` | nkf -w'

これならば tw とするだけで実行できる。また、PowerShellでは tw | select-string フレンド名、Unixでは tw | grep フレンド名で、指定したフレンドの発言をタイムラインから簡単に抜き出すこともできる。

今回のPythonによるTwitterクライアントのソースを以下に示す。simplejsonを使っている。bit.lyのアカウント名とAPI Keyはbit.lyアカウントのページですぐに取得できる。また、Base64でエンコードされたユーザ名とパスワードの作成は、Pythonのbase64モジュールを使ってもいいし、以前にGoogle App Engineで作ったバイナリ/アスキー変換を利用してもらっても構わない。まあ、エンコードされていてもセキュリティが強固になるわけではないのだが、ぱっと見でばれてしまう平文よりはいいんじゃないかな。ワンライナーにしたのもぱっと見で解りづらくするためだし。

ところで、TwitterのAPIでデータを取得したり送信したりすると、たまに失敗することがあるので、エラーが出たときは再度実行すること。

tw.py

#!/usr/bin/env python import sys,re,urllib,urllib2,xml.sax.saxutils,base64,simplejson;pm=urllib2.HTTPPasswordMgrWithDefaultRealm();pm.add_password(None,'twitter.com',base64.b64decode('エンコードされたユーザ名'),base64.b64decode('エンコードされたパスワード'));urllib2.install_opener(urllib2.build_opener(urllib2.HTTPBasicAuthHandler(pm)));tw=None;len(sys.argv)>1and[globals().__setitem__('tw',' '.join(sys.argv[1:]))];F=lambda x:str(simplejson.loads(urllib2.urlopen('http://api.bit.ly/shorten?version=2.0.1&longUrl=%s&login=bit.lyのアカウント名&apiKey=bit.lyのAPI Key'%urllib.quote(x)).read())['results'][x]['shortUrl']);G=lambda x:re.sub(re.escape(x),F(x),tw);tw and[globals().__setitem__('tw',G(link))for link in sorted(re.findall(r'(http(?:s?)\:\/\/[^\/\ ]+\/.*?)(?:[\ <>\"\{\}\|\\\^\[\]\`]|$)',tw),reverse=True)];tw and[urllib2.urlopen('http://twitter.com/statuses/update.xml',urllib.urlencode({'status':tw.decode('cp932').encode('utf-8')}))]or[sys.stdout.write(''.join(['%s: %s\n'%(d['user']['screen_name'],re.sub(r'\r\n|\n|\r',' ',xml.sax.saxutils.unescape(d['text'])))for d in simplejson.loads(urllib2.urlopen('http://twitter.com/statuses/friends_timeline.json').read())]).encode('cp932','replace'))]

追記(2009/9/9):

今までTwitterのお気に入り機能を使っていなかったけど、ふぁぼったーを見て使いたくなったので、今回のコマンドライン型Twitterクライアントでも使えるように改良した。

まず、tw.pyについてはつぶやきの最後に[つぶやきのID]を表示するようにした。また、fav.pyでつぶやきをお気に入りに加えることができる。使い方は以下の通り。

fav.py [つぶやきのID]

tw.pyで表示されたIDを入力する。コピー&ペーストだと簡単。IDの括弧([])は入れても入れなくても構わない。また、IDを付けずに

fav.py

とすれば、自分のお気に入りを表示することができる。

ただ、お気に入りを使わないのであればtw.pyによるIDの表示は目障りになるかもしれないので、そのときは以前のtw.pyを使った方がいいと思う。

以下、修正したソース。

tw.py

#!/usr/bin/env python import sys,re,urllib,urllib2,xml.sax.saxutils,base64,simplejson;pm=urllib2.HTTPPasswordMgrWithDefaultRealm();pm.add_password(None,'twitter.com',base64.b64decode('エンコードされたユーザ名'),base64.b64decode('エンコードされたパスワード'));urllib2.install_opener(urllib2.build_opener(urllib2.HTTPBasicAuthHandler(pm)));tw=None;len(sys.argv)>1and[globals().__setitem__('tw',' '.join(sys.argv[1:]))];F=lambda x:str(simplejson.loads(urllib2.urlopen('http://api.bit.ly/shorten?version=2.0.1&longUrl=%s&login=bit.lyのアカウント名&apiKey=bit.lyのAPI Key'%urllib.quote(x)).read())['results'][x]['shortUrl']);G=lambda x:re.sub(re.escape(x),F(x),tw);tw and[globals().__setitem__('tw',G(link))for link in sorted(re.findall(r'(http(?:s?)\:\/\/[^\/\ ]+\/.*?)(?:[\ <>\"\{\}\|\\\^\[\]\`]|$)',tw),reverse=True)];tw and[urllib2.urlopen('http://twitter.com/statuses/update.xml',urllib.urlencode({'status':tw.decode('cp932').encode('utf-8')}))]or[sys.stdout.write(''.join(['%s: %s [%d]\n'%(d['user']['screen_name'],re.sub(r'\r\n|\n|\r',' ',xml.sax.saxutils.unescape(d['text'])),d['id'])for d in simplejson.loads(urllib2.urlopen('http://twitter.com/statuses/friends_timeline.json').read())]).encode('cp932','replace'))]

fav.py

#!/usr/bin/env python import sys,re,urllib,urllib2,xml.sax.saxutils,base64,simplejson;pm=urllib2.HTTPPasswordMgrWithDefaultRealm();pm.add_password(None,'twitter.com',base64.b64decode('エンコードされたユーザ名'),base64.b64decode('エンコードされたパスワード'));urllib2.install_opener(urllib2.build_opener(urllib2.HTTPBasicAuthHandler(pm)));len(sys.argv)>1and[urllib2.urlopen('http://twitter.com/favorites/create/%d.xml'%int(sys.argv[1].strip('[]')),{})]or[sys.stdout.write(''.join(['%s: %s [%d]\n'%(d['user']['screen_name'],re.sub(r'\r\n|\n|\r',' ',xml.sax.saxutils.unescape(d['text'])),d['id'])for d in simplejson.loads(urllib2.urlopen('http://twitter.com/favorites.json').read())]).encode('cp932','replace'))]

続きを読む...

2009年8月28日金曜日

PythonのワンライナーでTwitterを使う

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

Twitterでつぶやいたり、タイムラインを取得したりするPythonのワンライナー(1行プログラム)を作ってみた。取り敢えずWindowsで動作は確認した。Pythonさえ入っていればどこでも動くと思う。シェルやcronに組み込んだり、ウェブアプリや自作プログラムで利用したり、Python以外に必要なものがないから手軽に使えるんじゃないかな。ただ、ユーザ名とパスワードは生テキストなのでその辺は気をつけるべきかも。

まず、Twitterでつぶやくワンライナー。

python -c "import urllib,urllib2;pm=urllib2.HTTPPasswordMgrWithDefaultRealm();pm.add_password(None,'twitter.com','username','password');urllib2.install_opener(urllib2.build_opener(urllib2.HTTPBasicAuthHandler(pm)));urllib2.urlopen('http://twitter.com/statuses/update.xml',urllib.urlencode({'status':'つぶやき'.decode('cp932').encode('utf-8')}))"

次に、タイムライン取得。simplejsonを使っている。

python -c "import sys,urllib,urllib2,xml.sax.saxutils,simplejson;pm=urllib2.HTTPPasswordMgrWithDefaultRealm();pm.add_password(None,'twitter.com','username','password');urllib2.install_opener(urllib2.build_opener(urllib2.HTTPBasicAuthHandler(pm)));sys.stdout.write(''.join(['%s: %s\n'%(d['user']['screen_name'],xml.sax.saxutils.unescape(d['text']))for d in simplejson.loads(urllib2.urlopen('http://twitter.com/statuses/friends_timeline.json').read())]).encode('cp932','replace'))"

因みに文字コードはCP932、簡単に言えばShift JISなので(厳密には違うけど)、環境に合わせて修正したり、nkfをかませるなどして欲しい。それにしても、TwitterのAPIは使いやすくていいね。

追記(2009/8/29):

Twitterのつぶやきとタイムライン取得の両方を行うワンライナーも書いてみた。このコードはPython 2.5以上で動作する。

python -c "import sys,urllib,urllib2,xml.sax.saxutils,simplejson;pm=urllib2.HTTPPasswordMgrWithDefaultRealm();pm.add_password(None,'twitter.com','username','password');urllib2.install_opener(urllib2.build_opener(urllib2.HTTPBasicAuthHandler(pm)));sys.stdout.write(''.join(['%s: %s\n'%(d['user']['screen_name'],xml.sax.saxutils.unescape(d['text']))for d in simplejson.loads(urllib2.urlopen('http://twitter.com/statuses/friends_timeline.json').read())]).encode('cp932','replace'))if len(sys.argv)<2 else urllib2.urlopen('http://twitter.com/statuses/update.xml',urllib.urlencode({'status':sys.argv[1].decode('cp932').encode('utf-8')}))" つぶやき

「つぶやき」を書けばそれがTwitterに送信され、書かなければタイムライン取得となる。また、「つぶやき」にスペースなどが入る場合はダブルクォーテーション(")で括ること。

次のコードであればPython 2.4以下でもいけるかな。

python -c "import sys,urllib,urllib2,xml.sax.saxutils,simplejson;pm=urllib2.HTTPPasswordMgrWithDefaultRealm();pm.add_password(None,'twitter.com','username','password');urllib2.install_opener(urllib2.build_opener(urllib2.HTTPBasicAuthHandler(pm)));len(sys.argv)<2and[sys.stdout.write(''.join(['%s: %s\n'%(d['user']['screen_name'],xml.sax.saxutils.unescape(d['text']))for d in simplejson.loads(urllib2.urlopen('http://twitter.com/statuses/friends_timeline.json').read())]).encode('cp932','replace'))]or[urllib2.urlopen('http://twitter.com/statuses/update.xml',urllib.urlencode({'status':sys.argv[1].decode('cp932').encode('utf-8')}))]" つぶやき

続きを読む...

2009年8月27日木曜日

Windows PowerShellを便利に使うための10のミニテクニック

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

Microsoftが提供していてWindows 7では標準搭載になるWindows PowerShellがもっと広まって欲しいこともあって、CodeZineWindowsのコマンドプロンプトを便利に使うための10のミニテクニックのPowerShell版を書いてみた。

まず、Windows PowerShell 2.0 CTP3をダウンロードして、インストールする。さらに便利に使うためにPowerShell Community ExtensionsからPSCX 1.2をダウンロード・インストールする。

さて、これで準備が整った。因みに以下のテクニックは、Windows XP、Vista、Windows 7のどのOSでも使えると思う(ただし、確認したのはXPのみ)。

コマンドプロンプトからエクスプローラに移動する

以下のように起動するだけ。

ii .

エクスプローラからコマンドプロンプトに移動する

これは標準では難しいと思う。自分はエクスプローラの使い勝手には非常に不満を持っているので、ずいぶん前からWindowsではFileVisorを使っている。FileVisorであれば以下のコマンドをホットキーに登録しておくことで、現在のディレクトリをカレントディレクトリとして一発でPowerShellを開くことができる。

C:\WINDOWS\system32\windowspowershell\v1.0\powershell.exe -NoExit -Command Set-Location $P

カレントディレクトリを記憶し、あとで戻ってくる

これはそのまま、pushd .popdが使える。

2つのディレクトリを行ったり来たりする

これもほとんど一緒。

doskey /exename=powershell.exe d1=cd $pwd

一時的にネットワークドライブを割り当てる必要はない

PowerShellではcdでそのままネットワークをまたげるので必要ないと思う。

cd \\computer1\project1\program1

処理結果をクリップボードにコピーする

ocb (Out-Clipboard)を使う。

dir | ocb tree | ocb

因みに、gcb (Get-Clipboard)を使うことで、PowerShell上に貼り付けもできる。

gcb

ずれたインデントを修正する

これもあまり変わらない。

more.com /t4 hello.c

ページごとに止めずに、一気に出力する。

more.com /t4 hello.c | echo

簡単な計算をする

setコマンドは使う必要はなく、そのまま計算させるだけ。

10-20+30 20

(10+20)*30/50 18

0xFF 255

因みに、小数も扱える。また、10進→16進変換は以下のようにすればできる。

"{0:X}" -f 255 FF

ファイルのタイムスタンプを変更する

UNIXのtouchコマンドと全く同じで以下のようにする。

touch file.ext

因みに、以下のようにもできる。

dir file.ext | % { $_.lastwritetime = get-date }

コマンドプロンプトの色やサイズを調整する

背景色を変える。

$host.ui.rawui.backgroundcolor = "blue"

画面の大きさを幅120桁、高さ50行に変える。

$host.ui.rawui.windowsize = new-object system.management.automation.host.size(120, 50)

画面のバッファを幅120桁、高さ3000行に変える。

$host.ui.rawui.buffersize = new-object system.management.automation.host.size(120, 3000)

おわりに

これを機にPowerShellが流行るといいなぁ。Windows 7では標準で付いてくるしね。

続きを読む...