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; }

続きを読む...