2011年12月23日金曜日

円周率を1万桁まで求める

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

マチンの公式で円周率を1万桁まで求めるプログラムを、Python, Erlang, Haskell, C++で書いてみた。出力結果の最後の数桁はずれているかも。また、実行時間を測ったりもしているが、1万桁程度ならどのコードでも瞬時に求まる。

まずは素直にPythonで。実行時オプションで桁数指定、実行時間測定付き。オプションなしで1万桁まで求める。

pi.py

#!/usr/bin/env python import sys, time N = 10**10000 def arctan(m): global N c = N a = b = c / m m2 = m * m s = k = 1 while c: b /= m2 k += 2 c, s = b / k, -s a += c * s return a def main(args): global N if len(args) > 1: N = 10**int(args[1]) t1 = time.time() pi = str((arctan(5) * 4 - arctan(239)) * 4) t2 = time.time() print pi[0] + '.' + pi[1:] print "Time: %f" % (t2 - t1) if __name__ == "__main__": main(sys.argv)

実行:

$ ./pi.py

次にErlangで求めてみる。

pi.erl

-module(pi). -export([pi/0]). pi()->N=e(10,10000),(a(5,N)*4-a(239,N))*4. e(B,N)->e(B,N,1). e(_,0,R)->R; e(B,N,R)->e(B,N-1,R*B). a(X,N)->a(X,N div X,N div X,N,1,1). a(_,A,_,0,_,_)->A; a(M,A,B,_,S,K)->B_=B div(M*M),C=B_ div(K+2),S_=-S,A_=A+C*S_,a(M,A_,B_,C,S_,K+2).

実行:

$ erl 1> c(pi). 2> pi:pi().

そしてHaskellを使う。ソースコードは一行だ。

pi.hs

n=10^10^4;p=a 5 n*4-a 239 n;a m n=t m(n`div`m)(n`div`m)n 1 3;t _ a _ 0 _ _=a;t m a b c s k=t m(a-y*s)x y(-s)(k+2)where x=b`div`m^2;y=x`div`k

実行:

$ ghci Prelude> :l pi Prelude> p*4

最後にC++のコード。多倍長演算でちょっと長い。実行時間測定付き。

pi.cpp

#include <iostream> #include <iomanip> #include <algorithm> #include <cmath> #ifdef _MSC_VER #include <ctime> inline double get_time() { return static_cast<double>(std::clock()) / CLOCKS_PER_SEC; } #else #include <sys/time.h> inline double get_time() { timeval tv; gettimeofday(&tv, 0); return tv.tv_sec + 1e-6 * tv.tv_usec; } #endif typedef unsigned long long number; const int FIG = 8; const number BASE = 0xffffffffULL; const number NN = 100000000ULL; // C = 1.0 / (4 * log10(2)) / FIG; #define C (0.830482023722/FIG) #define M (X/FIG+2) #define N (static_cast<int>(X*C+2)) template<int X> class PI { private: number ans[N]; public: void add(number* a, number* b); void sub(number* a, number* b); void div(number* a, number d); number intdiv(number x, number y); void mul(number* a, number d); void atan(number* a, number m); number or_numbers(number* a); bool has_bit(number* a); void decimal(number* a, number* w); void print(); void calc(); }; // a <- arctan(1 / m) template<int X> void PI<X>::atan(number* a, number m) { number b[N] = { 1 }; number c[N] = { 1 }; div(b, m); std::copy(b, b + N, a); const number m2 = m * m; for (int i = 1; has_bit(c); i++) { div(b, m2); std::copy(b, b + N, c); div(c, i * 2 + 1); if (i & 1) sub(a, c); else add(a, c); } } // a <- a + b template<int X> void PI<X>::add(number* a, number* b) { for (int i = 0; i < N; i++) { number x = a[i] + b[i]; if (x & ~BASE) { a[i] = x & BASE; int j; for (j = i - 1; a[j] == BASE; j--) a[j] = 0; a[j]++; } else a[i] = x; } } // a <- a - b (a > b) template<int X> void PI<X>::sub(number* a, number* b) { for (int i = 0; i < N; i++) { if (a[i] < b[i]) { a[i] = (BASE + 1) + a[i] - b[i]; int j; for (j = i - 1; a[j] == 0; j--) a[j] = BASE; a[j]--; } else a[i] -= b[i]; } } // a <- a / d (d <= BASE) template<int X> void PI<X>::div(number* a, number d) { number res = 0; for (int i = 0; i < N; i++) { res <<= (FIG * 4); number x = a[i] + res; number q = x / d; a[i] = q; res = x - q * d; } } template<int X> number PI<X>::intdiv(number x, number y) { const number maxbit = ~((number)-1>>1); number a = 0; number b = 0; int cnt = sizeof(number) * 8; while (cnt--) { a <<= 1ULL; if (x & maxbit) a |= 1ULL; x <<= 1ULL; b <<= 1ULL; if (a >= y) { ++b; a -= y; } } return b; } // a <- a * d (d <= BASE) template<int X> void PI<X>::mul(number* a, number d) { number q = 0; for (int i = N - 1; i >= 0; i--) { number x = a[i] * d + q; a[i] = x & BASE; q = x >> (FIG * 4); } } // OR numbers template<int X> number PI<X>::or_numbers(number* a) { number ret = 0; for (int i = 0; i < N; i++) ret |= a[i]; return ret; } template<int X> bool PI<X>::has_bit(number* a) { for (int i = 0; i < N; i++) if (a[i]) return true; return false; } // hex to decimal number template<int X> void PI<X>::decimal(number* a, number* w) { number b[N]; std::copy(a, a + N, b); w[0] = b[0]; b[0] = 0; for (int i = 1; i < M; i++) { mul(b, NN); w[i] = b[0]; b[0] = 0; } } // print answer template<int X> void PI<X>::print() { number w[M]; decimal(ans, w); std::cout << std::setw(4) << std::right << w[0] << "."; for (int i = 1; i < M; i++) { std::cout << std::setw(FIG) << std::setfill('0') << w[i] << " "; if(i % 6 == 0) std::cout << std::endl << " "; } std::cout << std::endl; } // calculate PI // Machin's formula: 4 * arctan(1/5) - arctan(1/239) = PI / 4 template<int X> void PI<X>::calc() { number x[N]; atan(ans, 5); mul(ans, 4); atan(x, 239); sub(ans, x); mul(ans, 4); } int main() { std::ios_base::sync_with_stdio(false); PI<10000> pi; double start_time = get_time(); pi.calc(); double end_time = get_time(); pi.print(); std::cerr << "Time: " << end_time - start_time << std::endl; return 0; }

実行:

$ c++ pi.cpp -O3 -o pi $ ./pi

続きを読む...

2011年11月30日水曜日

好きな・嫌いなコンピュータ言語のアンケート

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

社内で好きなコンピュータ言語と嫌いなコンピュータ言語のアンケートを取ってみた。ことの発端は自分のツイートなのだが、思いがけず賛同を得ることができたので、実際にアンケートを社内で行うことにしたのだ。各自で社内Wikiを編集してもらってもよかったのだけど、プログラムを作る会社でもあるし、投票システムに興味もあったので、Google App Engine (GAE)で作ってみることにした。以前にGAEで掲示板などを作成していたので、それらを流用して時間をかけずにすぐに仕上げることができた。以下のリンクから実際に投票できる。



表示されているコンピュータ言語一覧から、最も好きな言語と最も嫌いな言語をそれぞれ一つずつ、それ以外にも好きな言語や嫌いな言語があれば、それらも選択(複数可)する。投票したい言語が一覧にない場合は「言語の追加」で自由に追加できるようになっている。言語を選択して投票ボタンを押すと、投票数とそのグラフが更新される。表示されるグラフには、「最も好きな・最も嫌いな言語」と「好きな・嫌いな言語」の2つがあり、青のグラフが「好きな言語」、赤のグラフが「嫌いな言語」を表している。「好きな・嫌いな言語」のグラフについては、「最も好きな・最も嫌いな言語」と「好きな・嫌いな言語」を足し合わせた票数となっている。

また、コメントも投稿できるようになっていて、自分の好きな言語を啓蒙するのも良し、嫌いな言語について文句を言うのもアリだ。自由に書き込める。

現時点での社内での票は以下のグラフに示す通りで、幅広く票が入っているように思う。傾向としては、C, C++, C#, ASM, Python, Rubyあたりが人気で、Visual Basic, Java, Perlが不人気みたいだ。



今回作成したコードの簡単な説明を書いておく。上述の「好きな・嫌いなコンピュータ言語」だけでなく、初期変数を設定することで、「好きな・嫌いな動物」「好きな・嫌いな食べ物」などのように、簡単に任意のアンケートを作成できるようになっている。グラフについてはGoogle Chart APIを利用しているが、1,000ピクセルを超えることができないので、表示数を1,000ピクセル以内に抑える処理を入れている。また、頻繁にデータベースにアクセスしないようにmemcacheも利用している。詳細についてはコードを読んで欲しい。

以下に今回作成したGAEのコードを示しておく。

main.py

#!/usr/bin/env python # -*- coding: utf-8 -*- """main.py Copyright (C) 2011 nox""" import os, cgi, re, math, datetime import urllib, Cookie import wsgiref.handlers from google.appengine.api import memcache from google.appengine.ext import webapp from google.appengine.ext import db from google.appengine.ext.webapp import template # 許可する投票タイトル. TITLES = { "test": u"投票テスト", "lang": u"好きな・嫌いなコンピュータ言語のアンケート" } # 表示メッセージ. MESSAGES = { "test": u"投票をお願いします。", "lang": u"最も好きな言語と最も嫌いな言語をそれぞれ1つ、他に好きな言語や嫌いな言語があればそれも選択(複数可)して投票をお願いします。" } # 投票する属性. (タイトル / love / like / hate / dislike) PROPS = { "test": (u"テスト", u"最高", u"良い", u"最低", u"悪い"), "lang": (u"言語", u"最も好き", u"好き", u"最も嫌い", u"嫌い") } # チャートのタイトル. (love and hate / like and dislike [including love and hate]) CHARTS = { "test": (u"最高・最低のテスト", u"良い・悪いテスト"), "lang": (u"最も好きな・最も嫌いな言語", u"好きな・嫌いな言語") } # 最大表示件数 TARGET_LIMIT = 300 COMMENT_LIMIT = 100 class MainPage(webapp.RequestHandler): def get(self): self.response.out.write(u""" <html> <head> <title>Opinion Poll</title> </head> <body> <div><a href="vote">投票システム</a></div> </body> </html>""") # 投票データ. class VoteData(db.Model): id = db.StringProperty(multiline=False) content = db.StringProperty(multiline=False) love = db.IntegerProperty() like = db.IntegerProperty() hate = db.IntegerProperty() dislike = db.IntegerProperty() # コメントデータ. class CommentData(db.Model): id = db.StringProperty(multiline=False) username = db.StringProperty(multiline=False) content = db.TextProperty() date = db.DateTimeProperty(auto_now_add=True) # 投票システム・表示. class Vote(webapp.RequestHandler): def get(self): global TITLES, TARGET_LIMIT id = self.request.get("id") if id not in TITLES.keys(): return # データベースまたはMemcacheからデータの読み込み. cache_votes_key = "votes_%s" % id try: votes = memcache.get(cache_votes_key) except: try: memcache.delete(cache_votes_key) except: pass votes = None if votes is None: votes = VoteData.all().order("content").filter("id =", id).fetch(TARGET_LIMIT) try: memcache.add(cache_votes_key, votes) except: pass cache_comments_key = "comments_%s" % id try: comments = memcache.get(cache_comments_key) except: try: memcache.delete(cache_comments_key) except: pass comments = None if comments is None: comments = CommentData.all().order("-date").filter("id =", id).fetch(COMMENT_LIMIT) try: memcache.add(cache_comments_key, comments) except: pass # チャート表示 (1,000ピクセル以下にデータ数を収める). target2, target4 = [], [] target2_max, target4_max = 0, 0 for v in votes: v.love = int(v.love) if v.love else 0 v.hate = int(v.hate) if v.hate else 0 v.like = int(v.like) if v.like else 0 v.dislike = int(v.dislike) if v.dislike else 0 if v.love or v.hate or v.like or v.dislike: target4.append((v.content, v.love + v.like, v.hate + v.dislike)) target4_max = max(target4_max, v.love + v.like, v.hate + v.dislike) if v.love or v.hate: target2.append((v.content, v.love, v.hate)) target2_max = max(target2_max, v.love, v.hate) target2_max = (target2_max - 1) / 10 if target2_max > 0 else 0 target2_max = (target2_max + 1) * 10 target4_max = (target4_max - 1) / 10 if target4_max > 0 else 0 target4_max = (target4_max + 1) * 10 n = 0 while True: target2 = filter(lambda t: max(t[1:]) > n, target2) if len(target2) <= 16: break n += 1 if target2: text = "" if n > 0: text = u" (%d票以上)" % (n + 1) fig = 10**(int(math.log10(target2_max))) if target2_max / fig < 5: fig /= 2 chart_scale = "|".join([str(i) for i in range(0, target2_max + 1, fig)]) chart_target2_url = u"http://chart.apis.google.com/chart?chs=%dx200&chd=t:%s|%s&chds=0,%d&cht=bvg&chco=4d89f9,f9894d&chxt=y,x&chxl=0:|%s|1:|%s&chtt=%s%s" % (len(target2) * 60 + 30, ",".join([str(t[1]) for t in target2]), ",".join([str(t[2]) for t in target2]), target2_max, chart_scale, "|".join([urllib.quote(t[0].encode('utf-8')) for t in target2]), CHARTS[id][0], text) else: chart_target2_url = "" n = 0 while True: target4 = filter(lambda t: max(t[1:]) > n, target4) if len(target4) <= 16: break n += 1 if target4: text = "" if n > 0: text = u" (%d票以上)" % (n + 1) fig = 10**(int(math.log10(target4_max))) if target4_max / fig < 5: fig /= 2 chart_scale = "|".join([str(i) for i in range(0, target4_max + 1, fig)]) chart_target4_url = u"http://chart.apis.google.com/chart?chs=%dx200&chd=t:%s|%s&chds=0,%d&cht=bvg&chco=4d89f9,f9894d&chxt=y,x&chxl=0:|%s|1:|%s&chtt=%s%s" % (len(target4) * 60 + 30, ",".join([str(t[1]) for t in target4]), ",".join([str(t[2]) for t in target4]), target4_max, chart_scale, "|".join([urllib.quote(t[0].encode('utf-8')) for t in target4]), CHARTS[id][1], text) else: chart_target4_url = "" # 挨拶. H = (datetime.datetime.utcnow() + datetime.timedelta(hours=9)).hour if H >= 5 and H < 11: greeting = u"おはようございます。" elif H >= 11 and H < 17: greeting = u"こんにちは。" elif H >= 17 or H < 5: greeting = u"こんばんは。" else: greeting = u"こんにちは。" greeting += MESSAGES[id] # ウェブ上で表示させるデータ. template_values = { "id": id, "title": TITLES[id], "target": PROPS[id][0], "love": PROPS[id][1], "like": PROPS[id][2], "hate": PROPS[id][3], "dislike": PROPS[id][4], "votes": votes, "comments": comments, "chart_target2_url": chart_target2_url, "chart_target4_url": chart_target4_url, "greeting": greeting, } path = os.path.join(os.path.dirname(__file__), "vote.html") self.response.out.write(template.render(path, template_values)) # 投票ターゲットの追加. class AddVote(webapp.RequestHandler): def post(self): vote = VoteData() vote.id = self.request.get("id") try: try: memcache.delete("votes_%s" % vote.id) except: pass except: pass if not self.request.get("content").strip(): self.redirect("/vote?id=%s" % vote.id) return content = cgi.escape(self.request.get("content").strip()[:100]) v = VoteData.all().filter("id =", vote.id).filter("content =", content).get() if not v: vote.content = content vote.put() self.redirect("/vote?id=%s" % vote.id) # 投票. class PostVote(webapp.RequestHandler): def post(self): vote = VoteData() vote.id = self.request.get("id") love = self.request.get_all("love") like = self.request.get_all("like") hate = self.request.get_all("hate") dislike = self.request.get_all("dislike") try: try: memcache.delete("votes_%s" % vote.id) except: pass except: pass if not (love or like or hate or dislike): self.redirect("/vote?id=%s" % vote.id) return for d in love: v = VoteData.all().filter("id =", vote.id).filter("content =", d).get() if not v.love: v.love = 1 else: v.love = int(v.love) + 1 v.put() for d in like: v = VoteData.all().filter("id =", vote.id).filter("content =", d).get() if not v.like: v.like = 1 else: v.like = int(v.like) + 1 v.put() for d in hate: v = VoteData.all().filter("id =", vote.id).filter("content =", d).get() if not v.hate: v.hate = 1 else: v.hate = int(v.hate) + 1 v.put() for d in dislike: v = VoteData.all().filter("id =", vote.id).filter("content =", d).get() if not v.dislike: v.dislike = 1 else: v.dislike = int(v.dislike) + 1 v.put() self.redirect("/vote?id=%s" % vote.id) # コメントの投稿. class PostComment(webapp.RequestHandler): def insert_whitespace(self, matchobj): return re.sub("(?<!<br) ", "&nbsp;", matchobj.group(0).replace("\t", " ")) def post(self): global TARGET_LIMIT comment = CommentData() comment.id = self.request.get("id") try: try: memcache.delete("comments_%s" % comment.id) except: pass except: pass if not self.request.get("content").strip(): self.redirect("/vote?id=%s" % comment.id) return if self.request.get("username").strip(): comment.username = cgi.escape(self.request.get("username").strip()[:255]) else: comment.username = "Anonymous" comment.content = cgi.escape(self.request.get("content").strip()[:2000]) comment.content = re.sub("\r\n|\r|\n", "<br />", comment.content) comment.content = re.sub("\[link\](.*?)\[/link\]", "<a href=\"\g<1>\">\g<1></a>", comment.content) comment.content = re.sub("\[code\](.*?)\[/code\]", self.insert_whitespace, comment.content) comment.content = re.sub("\[code\](.*?)\[/code\]", "<span style=\"font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;\">\g<1></span>", comment.content) comment.date = comment.date + datetime.timedelta(hours=9) comment.put() self.redirect("/vote?id=%s" % comment.id) def main(): application = webapp.WSGIApplication([("/", MainPage), ("/vote", Vote), ("/vote/add", AddVote), ("/vote/post", PostVote), ("/vote/comment", PostComment), ], debug=False) wsgiref.handlers.CGIHandler().run(application) if __name__ == "__main__": main()

vote.html

<?xml version="1.0" encoding="UTF-8"?> <!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.1//EN" "http://www.w3.org/TR/xhtml11/DTD/xhtml11.dtd"> <html xmlns="http://www.w3.org/1999/xhtml" xml:lang="ja"> <head> <meta http-equiv="Content-Type" content="text/html; charset=UTF-8"> <link type="text/css" rel="stylesheet" href="/stylesheets/main.css" /> <title>{{ title }}</title> </head> <body> <script type="text/javascript"> function resize_textarea(ev){ var textarea = ev.target || ev.srcElement; var value = textarea.value; var lines = 1; for (var i = 0, l = value.length; i < l; i++) if (value.charAt(i) == '\n') lines++; if (lines > 5) { if (lines > 20) textarea.setAttribute("rows", 20); else textarea.setAttribute("rows", lines); } else textarea.setAttribute("rows", 5); } </script> <div style="font-weight: bold; font-size: large; padding-bottom: 20px;">{{ title }}</div> <div> <p>{{ greeting }}</p> </div> <form action="/vote/add" method="post"> <input type="text" name="content" size="20" maxlength="100"></input> <input type="submit" value="{{ target }}の追加" /> <input type="hidden" name="id" value="{{ id }}"> </form> <form action="/vote/post" method="post"> <div> <table border="0" width="400"> <th>{{ target }}</th><th>{{ love }}</th><th>{{ like }}</th><th>{{ hate }}</th><th>{{ dislike }}</th> {% for vote in votes %} <tr> <td><a href="http://ja.wikipedia.org/wiki/{{ vote.content }}">{{ vote.content }}</a></td> <td align="center"><input type="radio" name="love" value="{{ vote.content }}" /> {% if vote.love %} {{ vote.love }} {% else %} 0 {% endif %}</td> <td align="center"><input type="checkbox" name="like" value="{{ vote.content }}" /> {% if vote.like %} {{ vote.like }} {% else %} 0 {% endif %}</td> <td align="center"><input type="radio" name="hate" value="{{ vote.content }}" /> {% if vote.hate %} {{ vote.hate }} {% else %} 0 {% endif %}</td> <td align="center"><input type="checkbox" name="dislike" value="{{ vote.content }}" /> {% if vote.dislike %} {{ vote.dislike }} {% else %} 0 {% endif %}</td> </tr> {% endfor %} <tr align=center> <td colspan=5><input type="submit" value="投票" style="width: 100px; height: 25px" /></td> </tr> </table> <input type="hidden" name="id" value="{{ id }}"> </div> </form><br /> <div> {% if chart_target2_url %} <p><img src="{{ chart_target2_url }}" /></p> {% endif %} {% if chart_target4_url %} <p><img src="{{ chart_target4_url }}" /></p> {% endif %} </div> コメント記入: <span style="font-size:70%;color:#666666;">リンク: [link]~[/link], ソースコード: [code]~[/code], 1,000文字まで</span><br /> <form action="/vote/comment" method="post"> <div><textarea name="content" rows="5" cols="60" onkeyup="resize_textarea(event)"></textarea></div> 名前: <input type="text" name="username" size="31" maxlength="255" value="{{ username }}"></input> <input type="submit" value="コメント" /> <span style="font-size:70%;color:#666666;">無記名で匿名投稿</span> <input type="hidden" name="id" value="{{ id }}"> </form> <div> {% for comment in comments %} {{ comment.username }} [{{ comment.date.year }}年{{ comment.date.month }}月{{ comment.date.day }}日{{ comment.date.hour }}時{{ comment.date.minute }}分{{ comment.date.second }}秒]: <blockquote>{{ comment.content }}</blockquote> {% endfor %} </div> <!-- ここより下は変更しないでください。不具合やご意見などがありましたら、下記のnoxのリンク先で報告していただけると嬉しいです。 --> <div><img src="http://code.google.com/appengine/images/appengine-silver-120x30.gif" alt="Powered by Google App Engine" /> <span style="font-size:70%;color:#666666;"><a href="http://handasse.blogspot.com/">ソースコードおよびその解説</a></span><br /><span style="font-size:70%;color:#666666;">Copyright &copy; 2011 <a href="http://handasse.blogspot.com/">nox</a>. All rights reserved.</span></div> </body> </html>

続きを読む...

2011年10月19日水曜日

コンピュータにラストワンを解かせてみた

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

先日、地元のまつりがあった。そこの会場となっている区民館では子供用ボーリングやプラ板遊びなどの催しが開かれていたのだが、その中の一つに「ひとりぼっち」というゲームを行なっているのを見つけたので、一緒にいた小4の自分の娘に「このゲーム面白そうだね」と言ったところ、「ラストワンでしょ」と返され、よく知っているゲームとのことだった。遊び方は次の通り。全部で33のマスがあり、真ん中以外の32個のマスにコマを置く。コマは縦・横方向に隣のコマを一つ飛び越えて進むことができる。飛び越えられたコマは取り除かれる。これを繰り返して、最後のコマを真ん中のマスに置けばクリアとなる。32コマを使うステージ以外にもT字型やピラミッド型にコマを配置するなど、さまざまなバリエーションがある。


ネットで調べたところ、1985年にバンダイから同名の電子ゲームが出ていたようだ。リンク先では動画もあるのでそれを見てもらえば遊び方は一目瞭然だろう。普通に解いても面白そうだが、コンピュータを使って解かせたくなったので、早速プログラムを作ってみた。

余談だが、子供のころ、コンピュータは何でも解いてくれる魔法の箱だと思っていた。円周率を何桁まで解いたと聞いて、ワクワクしたものだ。ラストワンのようなゲームをコンピュータに解かせたいという気持ちは、子供のころのコンピュータに対する憧れから来るのかもしれない。

閑話休題。ラストワンは33マスあり、コマが「ある」「ない」の2状態のみを持つので、33ビットを使って盤面の状態を表すことができる。そして、ある状態から取り得る次の状態をすべて求めて、それを優先順位付きキューに入れて探索した。キューが空になるか、ゴールの状態が見つかるまで探索を行う。優先順位付きキューにしたのは探索を短縮するためにヒューリスティックな関数でスコアを決めようと思ったからなんだけど、32コマあるステージも含め、どれも一瞬で解けてしまうので、どうやら高度なヒューリスティック関数を用意しなくても問題なかったらしい。


上図のように初期配置としてT字型にコマを置く場合、以下に示すような入力データを用意して、これをプログラムに渡す。o (オー)でコマのあるマスを表し、. (ドット)で空のマスを表している。

lastone.dat

... ... ....... ..ooo.. ...o... .o. ...

入力データのファイル名が lastone.dat であれば、以下のように実行すればよい。

% ./lastone lastone.dat

実行結果は以下のようになる。

... ... ....... ..ooo.. ...o... .o. ... ... ... ....... .o..o.. ...o... .o. ... ... ... ....... .o.oo.. ....... ... ... ... ... ....... .oo.... ....... ... ... ... ... ....... ...o... ....... ... ...

最後に作成したコードを示す。

lastone.cpp

// Last One solver by nox, 14 Oct. 2011 #include <iostream> #include <fstream> #include <string> #include <vector> #include <queue> #include <algorithm> #include <utility> using namespace std; typedef unsigned long long state_t; typedef int score_t; const int N = 33; // number of cells class LastOne { private: vector<vector<pair<int, int> > > D; // table of directions state_t start, goal; vector<state_t> ans; // solution public: LastOne(const string datafile=""); ~LastOne() { } void read(const string datafile); int movable(state_t state); void next(state_t state, vector<state_t>& neighbors); bool search(); void show(state_t state); void answer(); }; // initialize and make table LastOne::LastOne(const string datafile) { if (!datafile.empty()) read(datafile); vector<pair<int, int> > v; v.push_back(make_pair(1, 2)); v.push_back(make_pair(3, 8)); D.push_back(v); v.clear(); v.push_back(make_pair(4, 9)); D.push_back(v); v.clear(); v.push_back(make_pair(1, 0)); v.push_back(make_pair(5, 10)); D.push_back(v); v.clear(); v.push_back(make_pair(4, 5)); v.push_back(make_pair(8, 15)); D.push_back(v); v.clear(); v.push_back(make_pair(9, 16)); D.push_back(v); v.clear(); v.push_back(make_pair(4, 3)); v.push_back(make_pair(10, 17)); D.push_back(v); v.clear(); v.push_back(make_pair(7, 8)); v.push_back(make_pair(13, 20)); D.push_back(v); v.clear(); v.push_back(make_pair(8, 9)); v.push_back(make_pair(14, 21)); D.push_back(v); v.clear(); v.push_back(make_pair(3, 0)); v.push_back(make_pair(7, 6)); v.push_back(make_pair(9, 10)); v.push_back(make_pair(15, 22)); D.push_back(v); v.clear(); v.push_back(make_pair(4, 1)); v.push_back(make_pair(8, 7)); v.push_back(make_pair(10, 11)); v.push_back(make_pair(16, 23)); D.push_back(v); v.clear(); v.push_back(make_pair(5, 2)); v.push_back(make_pair(9, 8)); v.push_back(make_pair(11, 12)); v.push_back(make_pair(17, 24)); D.push_back(v); v.clear(); v.push_back(make_pair(10, 9)); v.push_back(make_pair(18, 25)); D.push_back(v); v.clear(); v.push_back(make_pair(11, 10)); v.push_back(make_pair(19, 26)); D.push_back(v); v.clear(); v.push_back(make_pair(14, 15)); D.push_back(v); v.clear(); v.push_back(make_pair(15, 16)); D.push_back(v); v.clear(); v.push_back(make_pair(8, 3)); v.push_back(make_pair(14, 13)); v.push_back(make_pair(16, 17)); v.push_back(make_pair(22, 27)); D.push_back(v); v.clear(); v.push_back(make_pair(9, 4)); v.push_back(make_pair(15, 14)); v.push_back(make_pair(17, 18)); v.push_back(make_pair(23, 28)); D.push_back(v); v.clear(); v.push_back(make_pair(10, 5)); v.push_back(make_pair(16, 15)); v.push_back(make_pair(18, 19)); v.push_back(make_pair(24, 29)); D.push_back(v); v.clear(); v.push_back(make_pair(17, 16)); D.push_back(v); v.clear(); v.push_back(make_pair(18, 17)); D.push_back(v); v.clear(); v.push_back(make_pair(13, 6)); v.push_back(make_pair(21, 22)); D.push_back(v); v.clear(); v.push_back(make_pair(14, 7)); v.push_back(make_pair(22, 23)); D.push_back(v); v.clear(); v.push_back(make_pair(15, 8)); v.push_back(make_pair(21, 20)); v.push_back(make_pair(23, 24)); v.push_back(make_pair(27, 30)); D.push_back(v); v.clear(); v.push_back(make_pair(16, 9)); v.push_back(make_pair(22, 21)); v.push_back(make_pair(24, 25)); v.push_back(make_pair(28, 31)); D.push_back(v); v.clear(); v.push_back(make_pair(17, 10)); v.push_back(make_pair(23, 22)); v.push_back(make_pair(25, 26)); v.push_back(make_pair(29, 32)); D.push_back(v); v.clear(); v.push_back(make_pair(18, 11)); v.push_back(make_pair(24, 23)); D.push_back(v); v.clear(); v.push_back(make_pair(19, 12)); v.push_back(make_pair(25, 24)); D.push_back(v); v.clear(); v.push_back(make_pair(22, 15)); v.push_back(make_pair(28, 29)); D.push_back(v); v.clear(); v.push_back(make_pair(23, 16)); D.push_back(v); v.clear(); v.push_back(make_pair(24, 17)); v.push_back(make_pair(28, 27)); D.push_back(v); v.clear(); v.push_back(make_pair(27, 22)); v.push_back(make_pair(31, 32)); D.push_back(v); v.clear(); v.push_back(make_pair(28, 23)); D.push_back(v); v.clear(); v.push_back(make_pair(29, 24)); v.push_back(make_pair(31, 20)); D.push_back(v); } // read LastOne data void LastOne::read(const string datafile) { fstream fs(datafile.c_str(), ios_base::in); string line; goal = 0x10000; start = 0ULL; while (getline(fs, line)) { for (string::iterator p = line.begin(); p != line.end(); ++p) { if (*p == 'o') { start <<= 1; start |= 1ULL; } else if (*p == '.') start <<= 1; } } fs.close(); } // number of movable pieces int LastOne::movable(state_t state) { int s = 0; for (int i = 0; i < N; i++) if ((state >> i) & 1) for (vector<pair<int, int> >::iterator p = D[i].begin(); p != D[i].end(); ++p) s += ((state >> p->first) & 1) & !((state >> p->second) & 1); return s; } // next states void LastOne::next(state_t state, vector<state_t>& neighbors) { neighbors.clear(); for (int i = 0; i < N; i++) if ((state >> i) & 1) for (vector<pair<int, int> >::iterator p = D[i].begin(); p != D[i].end(); ++p) if (((state >> p->first) & 1) & !((state >> p->second) & 1)) neighbors.push_back(state ^ (1ULL << p->first) ^ (1ULL << p->second) ^ (1ULL << i)); } // search solution bool LastOne::search() { vector<state_t> checked(1, start); priority_queue<pair<score_t, vector<state_t> > > queue; queue.push(make_pair(1, checked)); vector<state_t> neighbors; while (!queue.empty()) { pair<score_t, vector<state_t> > q(queue.top()); queue.pop(); state_t last = q.second.back(); if (last == goal) { ans = q.second; return true; } next(last, neighbors); for (vector<state_t>::iterator p = neighbors.begin(); p != neighbors.end(); ++p) { if (binary_search(checked.begin(), checked.end(), *p)) continue; if (movable(*p) == 0 && *p != goal) continue; checked.insert(upper_bound(checked.begin(), checked.end(), *p), *p); vector<state_t> path(q.second); path.push_back(*p); queue.push(make_pair(static_cast<score_t>(path.size()), path)); } } return false; } // show state void LastOne::show(state_t state) { string line; while (state) { if (state & 1ULL) line.push_back('o'); else line.push_back('.'); state >>= 1; } int n = N - line.size(); for (int i = 0; i < n; i++) line.push_back('.'); reverse(line.begin(), line.end()); cout << " " << line.substr(0, 3) << endl; cout << " " << line.substr(3, 3) << endl; cout << line.substr(6, 7) << endl; cout << line.substr(13, 7) << endl; cout << line.substr(20, 7) << endl; cout << " " << line.substr(27, 3) << endl; cout << " " << line.substr(30, 3) << endl; } // show solution void LastOne::answer() { for (vector<state_t>::iterator p = ans.begin(); p != ans.end(); ++p) { show(*p); cout << endl; } } int main(int argc, char* argv[]) { if (argc < 2) { cerr << "Usage: " << argv[0] << " datafile" << endl; exit(1); } LastOne lastone(argv[1]); if (lastone.search()) lastone.answer(); else cout << "No solution." << endl; return 0; }

続きを読む...

2011年10月14日金曜日

オフィスビル内の最短経路

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

現在のオフィスでは、エントランスからどの通路を通るのが最も近いのか、未だに話題になることがある。エントランスからオフィスにつながるエレベータまでの道が二通りあり、どちらも似たような距離なので同じオフィスに行く人達なのにそこで別々になってしまうこともしばしば。左の図が問題のフロアなんだけど、エントランスから青丸で示したエレベータまで行く必要がある。だったらコンピュータに解かせてしまえばいいじゃない、ということでプログラムを作って最短経路と正確な距離を調べてみた。定規で測ったほうが早いという意見は聞こえません。

以前、Python: 画像で与えられた迷路に対し2点間の最短経路を求めるというブログ記事で、画像から直接最短距離を求めるプログラムを書いたことがあるので、まずはそれを使って調べてみた。図面の邪魔な文字だけ消して、あとはスタート地点とゴール地点の座標を引数で指定するだけだから簡単なんだけど、これだと問題があることに気がついた。というのはCCL (connected component labeling; 連結成分ラベリング)で移動できる領域を調べて、それをA*で解いているのだけど、隣接するノード(ピクセル)を上下左右斜めの8方向のみで調べているから正しい距離を出すことができないのだ。実際に作成した図は以下のようにギザギザになってしまい最短ルートを通らない。なので通り道の分岐や方向が変化する箇所をノードにしてそれをA*で求めることにした。


全部で10個のノード(画像のピクセル座標で指定)を作って調べたところ、以下のような結果になった。左側を通る通路の距離が324.53、右側を通る通路の距離が322.51となり僅かに右側の通路のほうが近かった。ただし、右側の方は途中でフロア案内と少し重なるし、6つあるエレベータのどれを使うかによっても結果が変わってしまうぐらいなので、実質的には同じ距離といってしまっても差し支えないと思う。



結局、それぞれの距離にはほとんど差がないので、結論としては「好きな方を使え」ということになるだろう。自分は広い通路の方が好みなので右側の通路を利用している。最後にソースコードとその使い方を示しておく。

nodes.datに、一行を"ピクセルX座標 ピクセルY座標 隣接ノード番号リスト"としてノード数の分だけ記述する。ノード番号は0から始まり、先頭行と最終行のノードはそれぞれスタート地点、ゴール地点を表している。また、フロアの図面画像をinput.pngとして渡すと、最短距離のルートが書きこまれた画像をoutput.pngとして出力する。

% ./route_solver.py nodes.dat input.png output.png

nodes.dat

27 273 1 3 6 21 128 0 2 125 23 1 9 117 221 0 4 5 6 141 178 3 7 167 178 3 7 195 178 0 3 7 192 110 4 5 6 8 193 91 7 9 151 41 2 8

route_solver.py

#!/usr/bin/env python # -*- coding: utf-8 -*- import sys, os, math, heapq, Image, ImageDraw class Node: def __init__(self): self.pos = (0, 0) self.neighbors = [] self.dist = float("inf") self.path = [] # A* Search Algorithm class AStar: def __init__(self, start, goal): self.start, self.goal = start, goal self.start.dist = 0.0 self.start.path = [self.start] def distance(self, node1, node2): return math.sqrt((node1.pos[0] - node2.pos[0])**2 + (node1.pos[1] - node2.pos[1])**2) def heuristic(self, node): return self.distance(node, self.goal) def search(self): queue = [] heapq.heappush(queue, (self.heuristic(self.start), self.start)) while queue: score, last = heapq.heappop(queue) if last == self.goal: return last.path, last.dist for nb in last.neighbors: newdist = last.dist + self.distance(last, nb) if nb.dist < newdist: continue nb.dist = newdist nb.path = last.path + [nb] heapq.heappush(queue, (nb.dist + self.heuristic(nb), nb)) return [], 0.0 def draw_path(in_image, out_image, path): im = Image.open(in_image) im = im.convert("RGB") draw = ImageDraw.Draw(im) draw.line([(int(p.pos[0]), int(p.pos[1])) for p in path], fill=255, width=2) im.save(out_image) def main(args): if len(args) < 4: print >>sys.stderr, "Usage: %s datafile in_image out_image" % os.path.basename(args[0]) sys.exit(1) datafile, in_image, out_image = args[1:4]; print "Reading nodes data..." data = [map(int, l.strip().split()) for l in file(datafile).readlines()] nodes = [Node() for i in range(len(data))] for i, node in enumerate(nodes): node.pos = tuple(data[i][:2]) for j in data[i][2:]: node.neighbors.append(nodes[j]) print " Number of nodes: %d" % len(nodes) print "A* searching phase..." start, goal = nodes[0], nodes[-1] print " Start: (%d, %d), Goal: (%d, %d)" % (start.pos + goal.pos) astar = AStar(start, goal) path, dist = astar.search() if path: print " Found route:" for p in path: print " (%d, %d)" % p.pos print " Distance: %f" % dist else: print " No route." print "Drawing path on the imagefile..." draw_path(in_image, out_image, path) if __name__ == "__main__": main(sys.argv)

追記 (2011/10/17): ソースコードが間違っていたので修正。今回の結果には影響なし。

続きを読む...

2011年9月27日火曜日

GoogleのURL短縮サービスgoo.glのAPIをPythonで利用する

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

しばらく前からgoo.glのAPI (Google URL Shortener API)が使えるようになっていたことは知っていたけど、実際に使ったことなかったので試しに使ってみた。import文を含めなければ3行で短縮URLの作成・復元の自動判別を行うことができた。条件演算子などを使えば1行にまとめられる。簡単。

ただし、正規表現は適当だし、エラー処理はしてないので、ちゃんと使う場合はその辺の処理を行う必要あり。あと、本格的に使うのであればAPIキーを取得すべき。また、Google APIs Client Library for Python使う方法もある。

http://handasse.blogspot.com/ を短縮URLに変換:

% googl.py http://handasse.blogspot.com/ http://goo.gl/RiIiD

http://goo.gl/RiIiD を元のURLに戻す:

% googl.py http://goo.gl/RiIiD http://handasse.blogspot.com/

googl.py

import sys,re from urllib2 import urlopen as U, Request as R from json import loads as J API,URL="https://www.googleapis.com/urlshortener/v1/url",sys.argv[1] if re.match('http://goo\.gl/.+',URL):print J(U(API+'?shortUrl=%s'%URL).read())['longUrl'] else:print J(U(R(API,'{"longUrl":"%s"}'%URL,{'Content-Type':'application/json'})).read())['id']

追記 (2015/3/26):

現在、API Keyを利用しないとエラーが返ってくるようなので、以下のように修正した。


コード中のYOUR_API_KEYの部分を実際のAPI Keyと差し替えれば利用できる。API Keyの取得するためには、まず、Google Developers Consoleにアクセスし、適当なプロジェクトに入り、APIと認証のAPIを選択し、Other popular APIsの中のURL Shortener APIを選ぶ。そこのページでEnable APIのボタンをクリックし、「APIコンソールでレポートを表示する」のリンクが表示されるので、そのページにアクセスする。メニューのAPI Accessを選び、そこのCreate new Browser key...ボタンをクリックする。これでAPI Keyが作成されるので、これを上述のYOUR_API_KEYと差し替えれば利用できるようになる。

続きを読む...

2011年8月8日月曜日

暇潰しで書いたプログラム

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

ちょっと出先で手持ち無沙汰になったとき、なにかコードが書きたくなることってあると思う。その時に書くプログラムは何でもよくて、言語も選ばない。でも、駅で電車を待っている間とか、注文したメニューが来るまでとか、ノートPCを引っ張り出してまでやりたくはない。そんなとき、Android端末を片手にSL4Aでコードを書くっていうのは丁度よい。

というわけで、ちょっと外出したときの暇な時間にSL4AのPythonで書いたコードが以下のコード。本当にくだらないプログラムなわけだけど、暇潰しにはなった。何をするプログラムか興味のある方は読み解いて欲しい。暇潰しぐらいにしかならないけど。まあ、勘のよい人ならファイル名だけでわかってしまうかも。

Android端末を持っている方なら、SL4Aのメニューを開いてAddを選択し、Scan Barcodeで冒頭のQRコードを読み込むことでソースコードを取得できる。SL4Aで書いたプログラムだけど、Android端末以外でもPythonの動作する環境なら大抵は実行できると思う。

Zm9yIGkgaW4gcmFuZ2UoMSwxMDEpOnByaW50J0ZpenonKigxLWklMykrJ0J1enonKigxLWklNSlvciBp.py

exec(r"""%s'''gzga*%%%%,hmkl*ocr*nco`fc"a8ajp*mpf*a+\\0+.p%%%s'kormpv"mq.`cqg469gzga*`cqg46,fgamfgqvpkle*mq,rcvj,`cqglcog*]]dkng]]+Y8/1_++')))%%+++''')))"""%(lambda _:(_,_))("exec(''.join(map(lambda _:chr(ord(_)^2),"))

続きを読む...

2011年7月11日月曜日

Scalaのアクターモデルでマンデルブロ集合を並列計算

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

最近、本格的にScalaを使用するつもりで環境を整えた。ビルドツールとしてSimple Build Tool (sbt)を利用し、エディタはensimeを入れたEmacsにしている。利用するに当たって.emacsに必要な記述を加えておくこと。初めてのプロジェクトの場合、適当に作成したディレクトリ内でsbtを起動してプロジェクトの雛形を作った後、Emacsを立ち上げて M-x ensime-config-gen で下準備をする。ここまでは初回のみの作業となる。その後は M-x ensime を実行して、src/main/scala/内でコードを書くだけ。これで、コーディング中にタブで補完してくれるし、文法などが間違っている場合に赤でハイライトもしてくれる。コンパイルや実行などは、C-c C-v s とすればsbtがEmacs上で立ち上がるので compile や run するだけ。簡単。

Scalaと言えばアクターでしょ、やっぱり。ということで今回はアクターモデルを使用したコードを書いてみた。自分はScalaをいじるよりも前にErlangをいじっていたので余計にそう考えてしまうのかもしれない。それでもマルチコアやらメニーコアが溢れているこの時代、並行処理の重要性は誰もが認識しているはず。

書いたコードはマンデルブロ集合の計算だ。毎回飽きもせずにまたそれかと思うかもしれないが、並列処理させるのが簡単だし、いつも同じ内容にしておけば比較も容易。なによりもマンデルブロ集合の画像はともて魅力的だと自分では思っているので、今回もマンデルブロ集合だ。冒頭に作成した画像を載せたけど、やっぱり美しいと思う。

今回のアクターモデルでは、scala.actors.Actorから派生させずに、scala.actors.Actor.actorを使用した。クラスを作るほど大行ではないし、その方がスッキリするからだ。そして、計算のループ部分は末尾再帰とパターンマッチを使っている。関数型言語として使うならやっぱりその方がしっくりくる。ただ、ちゃんと末尾再帰させるように、final修飾子をつけるなど、気をつけること。それと@tailrecアノテーション(scala.annotation.tailrec)は常に付けたほうが良い。@tailrecアノテーションを末尾再帰させたい関数につけておけば、末尾再帰としてコンパイルされなかった時にエラーを返してくれる。アクターを使用しないコードもついでに作成した。画像の色付け部分はScalaでMandelbrot、時間測定は関数の実行時間を計測するには?を参考にさせてもらった。

実行環境としてXeon X5650@2.67GHz * 2 (12コア/24スレッド; メモリ48GB)を使用した。マンデルブロ集合の条件としては、画像の大きさが4,800 x 4,800ピクセルで、(-0.005, -0.005)から(0.005, 0.005)を計算の範囲とし、繰り返し数10,000で10.0を閾値としている。因みにこの範囲は集合に含まれるので計算量としては最大となる(画像としては真っ黒なので面白くはないが)。実行すると最後に画像ファイル(mandelbrot.png)を出力して終了する。

実行した結果、アクターを使用しなかったコードの実行時間が1007.18秒だったのに対し、アクターを使用したコードでは51.636秒であった。約20倍の高速化である。コア数12・スレッド数24であることを考えても、この高速化は十分なのではないだろうか。ところで、アクターを使用しないシングルスレッドでのコードでも以前作成したC++のコード(1067.39秒)より速かった。それに、TBBを使用した並列化したコードでも56.8956秒の計算時間がかかっており、アクターを使用したScalaのコードはそれよりも速い。色の作成の部分など少し異なるところもあるのでそのまま比較することはできないかもしれないが、少なくとも今回の場合ではScalaで十分な速度を得られることは確認できた。余談だが、Telsaのボードを2枚挿ししたマシン上で実行したCUDAのコードでは3.73秒なので次元が違ったりする。

しかし、実際にコードを書くときは注意深くなる必要があるかもしれない。例えば以下のコードがある。

i match { case MAX_LOOP => -1 case _ if zr2 + zi2 >= TH => i case _ => { val zrzi = zr * zi calcLoop(i + 1, cr, ci, zr2 - zi2 + cr, zrzi + zrzi + ci) } }

これを以下のように変更する。

(i, zr2 + zi2) match { case (MAX_LOOP, _) => -1 case (_, z) if z >= TH => i case _ => { val zrzi = zr * zi calcLoop(i + 1, cr, ci, zr2 - zi2 + cr, zrzi + zrzi + ci) } }

これだけのことで実行速度が6倍程遅くなる。これに気付かずにいると本来のパフォーマンスを出すことができない。予想よりもプログラムが遅いと感じた場合は再度コードを確認したほうがいいだろう。

自分は普段は主にC++とPythonを使っていて、目的によって使い分けている。ある程度しっかりした作りで高速に動作させたい場合はC++だし、プロトタイプとして作ったり、結果の解析用コードをサクっと作りたい場合はPythonを使う。Scalaはどのような用途に使えるのか。まず、Javaの資産を使えるのは大きい。ネットワーク処理や画像処理などをある程度本格的につくろうとしたときに便利そうだ。あとは言語の仕様として処理を簡潔に書き表せるという点が良い。ライブラリが強力で簡潔に書ける。この点だけで十分に利用する価値があるだろう。

最後に自分の持っているScalaの書籍を2つ挙げておく。Scalaプログラミング入門は実践的だし、Scalaスケーラブルプログラミングは網羅的なので、どちらも持っていて損はないと思う。

作成したコードを以下に示しておく。

mandelbrot.scala (アクターあり; 並行処理)

// mandelbrot.scala with actor by nox, 2011.07.09 import scala.actors._, Actor._ import scala.annotation.tailrec import java.io.FileOutputStream import java.awt.image.BufferedImage import java.awt.Color import javax.imageio.ImageIO import scala.testing.Benchmark /** * マンデルブロ集合. * SW : 画像の幅. * SH : 画像の高さ. * T : 集合の上の位置. * L : 集合の左の位置. * W : 集合の幅. * H : 集合の高さ. * MAX_LOOP : イテレーションの最大回数. * TH : 閾値. */ final class Mandelbrot(SW: Int, SH: Int, T: Double, L: Double, W: Double, H: Double, MAX_LOOP: Int, TH: Double) { /** * イテレーションの回数から色を決定して任意のピクセルにセットする. */ def setPixel(x: Int, y: Int, a: Any, image: BufferedImage) = a match { case -1 => image setRGB(x, y, 0) case i: Int => image setRGB(x, y, Color HSBtoRGB(i / 100.0f, 1.0f, 1.0f)) } /** * マンデルブロ集合の計算部分. */ @tailrec def calcLoop(i: Int, cr: Double, ci: Double, zr: Double, zi: Double): Int = { val zr2 = zr * zr val zi2 = zi * zi i match { case MAX_LOOP => -1 case _ if zr2 + zi2 >= TH => i case _ => { val zrzi = zr * zi calcLoop(i + 1, cr, ci, zr2 - zi2 + cr, zrzi + zrzi + ci) } } } /** * アクター. * ある幅に存在するピクセルの計算を受け持つ. */ def calc = actor { react { case (ix: Int, image: BufferedImage) => { for (iy <- 0 until SH) setPixel(ix, iy, calcLoop(0, L + (ix.toDouble / SW) * W, T + (iy.toDouble / SH) * H, 0.0, 0.0), image) reply(ix) } } } /** * 幅のピクセル数分だけアクターにメッセージを送る. */ @tailrec def runLoop(ix: Int, results: Array[Future[Any]], image: BufferedImage): Unit = ix match { case SW => results foreach(_ apply) case _ => { results(ix) = calc !! ((ix, image)) runLoop(ix + 1, results, image) } } /** * 実行. * 計算部分の処理時間計測および画像ファイルの書き出し. */ def run = { val results = new Array[Future[Any]](SW) val image = new BufferedImage(SW, SH, BufferedImage TYPE_3BYTE_BGR) val t = (new Benchmark { def run = runLoop(0, results, image) } runBenchmark 1)(0) / 1000.0 println("Time: " + t + " s") val out = new FileOutputStream("mandelbrot.png") ImageIO write(image, "png", out) out close } } object MandelbrotMain { def main(args: Array[String]) = new Mandelbrot(4800, 4800, -0.005, -0.005, 0.01, 0.01, 10000, 10.0) run //new Mandelbrot(640, 480, -1.0, -2.0, 2.6666667, 2.0, 1000, 10.0) run //new Mandelbrot(800, 800, 0.025185, -1.401565, 0.0005, 0.0005, 10000, 10.0) run }

mandelbrot.scala (アクターなし; 逐次処理)

// mandelbrot.scala without actor by nox, 2011.07.09 import scala.annotation.tailrec import java.io.FileOutputStream import java.awt.image.BufferedImage import java.awt.Color import javax.imageio.ImageIO import scala.testing.Benchmark /** * マンデルブロ集合. * SW : 画像の幅. * SH : 画像の高さ. * T : 集合の上の位置. * L : 集合の左の位置. * W : 集合の幅. * H : 集合の高さ. * MAX_LOOP : イテレーションの最大回数. * TH : 閾値. */ final class Mandelbrot(SW: Int, SH: Int, T: Double, L: Double, W: Double, H: Double, MAX_LOOP: Int, TH: Double) { /** * イテレーションの回数から色を決定して任意のピクセルにセットする. */ def setPixel(x: Int, y: Int, a: Any, image: BufferedImage) = a match { case -1 => image setRGB(x, y, 0) case i: Int => image setRGB(x, y, Color HSBtoRGB(i / 100.0f, 1.0f, 1.0f)) } /** * マンデルブロ集合の計算部分. */ @tailrec def calcLoop(i: Int, cr: Double, ci: Double, zr: Double, zi: Double): Int = { val zr2 = zr * zr val zi2 = zi * zi i match { case MAX_LOOP => i case _ if zr2 + zi2 >= TH => i case _ => { val zrzi = zr * zi calcLoop(i + 1, cr, ci, zr2 - zi2 + cr, zrzi + zrzi + ci) } } } def calc(cr: Double, ci: Double): Int = calcLoop(0, cr, ci, 0.0, 0.0) /** * 画像のピクセル数分だけループする. */ @tailrec def runLoop(ix: Int, iy: Int, image: BufferedImage): Unit = (ix, iy) match { case (SW, _) => None case (_, SH) => runLoop(ix + 1, 0, image) case (_, _) => { setPixel(ix, iy, calc(L + (ix.toDouble / SW) * W, T + (iy.toDouble / SH) * H), image) runLoop(ix, iy + 1, image) } } /** * 実行. * 計算部分の処理時間計測および画像ファイルの書き出し. */ def run = { val image = new BufferedImage(SW, SH, BufferedImage TYPE_3BYTE_BGR) val t = (new Benchmark { def run = runLoop(0, 0, image) } runBenchmark 1)(0) / 1000.0 println("Time: " + t + " s") val out = new FileOutputStream("mandelbrot.png") ImageIO write(image, "png", out) out close } } object MandelbrotMain { def main(args: Array[String]) = new Mandelbrot(4800, 4800, -0.005, -0.005, 0.01, 0.01, 10000, 10.0) run //new Mandelbrot(640, 480, -1.0, -2.0, 2.6666667, 2.0, 1000, 10.0) run //new Mandelbrot(800, 800, 0.025185, -1.401565, 0.0005, 0.0005, 10000, 10.0) run }

続きを読む...

2011年6月6日月曜日

CUDA 4.0でマルチGPU化が限りなく簡単になっているという話

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

以前、複数のGPUでマンデルブロ集合を並列計算というブログ記事を書いた。このときはTeslaが1枚で7.45秒、Teslaが2枚で4.26秒で計算できたと報告したが、マルチGPU化は多少面倒だったし、ある程度のオーバーヘッドが掛かっていた。

しかし、CUDA 4.0がリリースされてから状況が一変した。自由にcudaSetDeviceが呼び出せるようになり、またUVA (Unified Virtual Addressing)を利用することにより、それはそれは簡単にマルチGPU化や、ホスト・複数GPU間での一元的なメモリの参照を実現できることになった。

というわけで、CUDA 4.0を使ったマンデルブロ集合のプログラムを以下に示す。シングルGPUからマルチGPUへ変更した部分については赤字で示してある。ただし、簡便のために、GPUデバイスは2枚で決め打ちとし、マンデルブロ集合の幅は2の倍数のピクセルとなるようにしてある。また、UVAは使わなくても良かったのだが、使い方を示すために利用している。

CUDA 4.0によるマルチGPU化によりマンデルブロ集合計算の実行時間(計算条件は前回と同じ)は7.45秒から3.73秒となった。ほぼ2倍の速度である。シングルGPUからの修正箇所はほんの僅か。カーネル関数に至っては一文字たりとも変更していない。良い時代になったものだ。

mandelbrot_multigpu.cu

// madelbrot using multi GPU devices by nox, 2011.06.06 // nvcc -lcutil_x86_64 -arch sm_13 -use_fast_math -prec-sqrt=false -keep -L ~/NVIDIA_GPU_Computing_SDK/C/lib -I ~/NVIDIA_GPU_Computing_SDK/C/common/inc -g mandelbrot_multigpu.cu -o mandelbrot_multigpu #include <iostream> #include <fstream> #include "cutil_inline.h" using namespace std; const int BLOCK_SIZE_X = 16; const int BLOCK_SIZE_Y = 16; __device__ uchar4 coloring(int n, int b) { int d = n % b; d *= 256 / b; int m = static_cast<int>(d / 42.667f); switch (m) { case 0: return make_uchar4(0, 6 * d, 255, 0); case 1: return make_uchar4(0, 255, 255 - 6 * (d - 43), 0); case 2: return make_uchar4(6 * (d - 86), 255, 0, 0); case 3: return make_uchar4(255, 255 - 6 * (d - 129), 0, 0); case 4: return make_uchar4(255, 0, 6 * (d - 171), 0); case 5: return make_uchar4(255 - 6 * (d - 214), 0, 255, 0); default: return make_uchar4(0, 0, 0, 0); } } __global__ void mandelbrot(float t, float l, float w, float h, int sw, int sh, int max_iter, float th, uchar4* d_color) { int ix = blockIdx.x * blockDim.x + threadIdx.x; int iy = blockIdx.y * blockDim.y + threadIdx.y; if (ix >= sw || iy >= sh) return; float ci = t + (static_cast<float>(iy) / sh) * h; float cr = l + (static_cast<float>(ix) / sw) * w; float zi = 0.0f; float zr = 0.0f; float zrzi, zr2, zi2; for (int i = 0; i < max_iter; i++) { zrzi = zr * zi; zr2 = zr * zr; zi2 = zi * zi; zr = zr2 - zi2 + cr; zi = zrzi + zrzi + ci; if (zi2 + zr2 >= th) { d_color[ix*sh+iy] = coloring(i, 256); return; } } d_color[ix*sh+iy] = make_uchar4(0, 0, 0, 0); } void write_mandelbrot(const string outfile, float t, float l, float w, float h, int sw, int sh, int max_iter=256, float th=2.0f) { dim3 num_blocks((sw - 1) / BLOCK_SIZE_X + 1, (sh - 1) / BLOCK_SIZE_Y + 1, 1); dim3 num_threads(BLOCK_SIZE_X, BLOCK_SIZE_Y, 1); uchar4* h_color; uchar4* d_color[2]; cudaSetDevice(0); cudaMallocHost((void**)&h_color, sizeof(uchar4) * sw * sh); for (int i = 0; i < 2; i++) { cudaSetDevice(i); cudaMalloc((void**)&d_color[i], sizeof(uchar4) * sw * sh / 2); } unsigned int timer; cutCreateTimer(&timer); cutStartTimer(timer); for (int i = 0; i < 2; i++) { cudaSetDevice(i); mandelbrot<<<num_blocks, num_threads>>>(t, l + (w / 2) * i, w / 2, h, sw / 2, sh, max_iter, th, d_color[i]); } cudaThreadSynchronize(); cutStopTimer(timer); cout << "Time: " << cutGetTimerValue(timer) / 1000 << endl; for (int i = 0; i < 2; i++) cudaMemcpy(h_color + sw * sh / 2 * i, d_color[i], sizeof(uchar4) * sw * sh / 2, cudaMemcpyDefault); fstream fs(outfile.c_str(), ios_base::out); fs << sw << endl; fs << sh << endl; for (int i = 0; i < sw * sh; i++) fs << h_color[i].x << h_color[i].y << h_color[i].z; fs.close(); for (int i = 0; i < 2; i++) { cudaSetDevice(i); cudaFree(d_color[i]); } cudaSetDevice(0); cudaFreeHost(h_color); } int main(int argc, char* argv[]) { string outfile("mandelbrot.dat"); if (argc >= 2) outfile = argv[1]; // for UVA cudaSetDevice(0); cudaDeviceEnablePeerAccess(1, 0); cudaSetDevice(1); cudaDeviceEnablePeerAccess(0, 0); // write_mandelbrot(output filename, // top position, // left position, // width, // height, // screen width, // screen height, // max number of iterations [256] // threshold [2.0]); //write_mandelbrot(outfile, -1.0f, -1.5f, 2.0f, 2.0f, 480, 480, 100); //write_mandelbrot(outfile, 0.57f, 0.34f, 0.036f, 0.027f, 6400, 4800, 1000); write_mandelbrot(outfile, -0.005f, -0.005f, 0.01f, 0.01f, 4800, 4800, 10000); return 0; }

出力ファイルの画像への変換は複数のGPUでマンデルブロ集合を並列計算の記事に記載されているread_mandelbrot.pyが利用できる。

続きを読む...

2011年5月13日金曜日

C++: Base64エンコード・デコードとURLエンコード

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

覚書。Base64エンコード・デコードを行う関数とURLエンコードを行う関数をどんな環境でもすぐに参照できるように書いておく。ただし、短いコードのままにしておきたかったので、まったくもって安全ではない。入力前にチェックを入れるとか、関数内で制限を掛けるとか、引数に最大文字数を渡すとか、適当な処理を入れること。よく分からない場合は使わないのが吉。

#include <iostream> #include <iomanip> #include <sstream> #include <algorithm> #include <cstring> using namespace std; int base64encode(const char* p, char* buf) { const char b64[] = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/="; const int in_size = strlen(p); if (in_size == 0 || p == buf) return 0; int cnt = 0; for (int i = 0; i < (in_size - 1) / 3 + 1; i++) { const unsigned char* c = reinterpret_cast<const unsigned char*>(&p[i*3]); for (int j = 0; j < 4; j++) if (in_size >= i * 3 + j) buf[cnt++] = b64[(c[j-1]<<(6-j*2) | c[j]>>(j*2+2)) & 0x3f]; else buf[cnt++] = '='; } buf[cnt] = '\0'; return cnt; } int base64decode(const char* p, char* buf) { char b64[128] = {}; for (char n = 0, *w = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/="; *w; b64[*w++] = n++ & 0x3f); char c[4]; int i = 0, j, k; while (*p && *p != '=') { for (j = 0; j < 4 && *p != '='; j++) c[j] = b64[*p++]; for (k = 0; k < j - 1; k++) buf[i++] = c[k]<<(k<<1)+2 | c[k+1]>>(2-k<<1); } buf[i] = '\0'; return i; } string urlencode(const string& text) { const string url_letters("-.0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ_abcdefghijklmnopqrstuvwxyz~"); stringstream ss; for (string::const_iterator p = text.begin(); p != text.end(); ++p) { if (binary_search(url_letters.begin(), url_letters.end(), *p)) ss << *p; else ss << "%" << setw(2) << setfill('0') << hex << static_cast<int>(*p); } return ss.str(); } int main() { // Base64エンコードでは入力文字数の4/3倍ほどの領域が必要. // 文字数が少なかったり改行を入れるともっと必要. char buf[256]; // 256は安全ではない. base64encode("test", buf); // (buf, buf)はダメ. cout << buf << endl; // dGVzdA== base64decode(buf, buf); // 邪悪. cout << buf << endl; // test cout << urlencode("user@test.com") << endl; // user%40test.com return 0; }

続きを読む...

2011年4月24日日曜日

タンパク質分子において水素原子の位置が異なるアミノ酸構造の同型判定

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

今回は興味を持ってもらえる人が限定される記事だと思うけど、以前作ったコードを腐らせてしまうのも勿体無いので敢えて晒してみる。

生体内のタンパク質構造をX線回折やNMRで解析したデータの多くはPDB (Protein Data Bank)データとして登録される。以下に示すように原子タグ、原子の通し番号、原子名、アミノ酸名、アミノ酸の通し番号、原子の座標などがテキストで書かれている。そして、機能の解明などのために分子動力学(MD; molecular dynamics)計算プログラムなどでそれらのPDBデータは利用される。MD計算プログラムは世に沢山あるが、自分がよく使っていたのはAmberというソフトウェアだ。

PDBデータの一部

ATOM 1439 N ILE 212 52.408 -3.243 59.013 1.00 40.03 N ATOM 1440 CA ILE 212 52.511 -3.151 57.556 1.00 33.72 C ATOM 1441 C ILE 212 51.464 -2.243 56.911 1.00 34.65 C ATOM 1442 O ILE 212 51.733 -1.619 55.886 1.00 38.03 O ATOM 1443 CB ILE 212 52.482 -4.572 56.921 1.00 38.88 C ATOM 1444 CG1 ILE 212 53.615 -4.703 55.919 1.00 40.75 C ATOM 1445 CG2 ILE 212 51.129 -4.901 56.281 1.00 39.55 C ATOM 1446 CD1 ILE 212 54.931 -4.446 56.550 1.00 46.31 C ATOM 1447 N HIS 213 50.308 -2.142 57.570 1.00 36.88 N ATOM 1448 CA HIS 213 49.160 -1.334 57.169 1.00 37.60 C ATOM 1449 C HIS 213 49.523 0.121 56.879 1.00 39.63 C ATOM 1450 O HIS 213 48.887 0.773 56.028 1.00 39.63 O ATOM 1451 CB HIS 213 48.115 -1.389 58.289 1.00 44.39 C ATOM 1452 CG HIS 213 46.795 -0.746 57.947 1.00 55.55 C ATOM 1453 ND1 HIS 213 46.364 -0.536 56.651 1.00 54.19 N ATOM 1454 CD2 HIS 213 45.800 -0.306 58.748 1.00 53.97 C ATOM 1455 CE1 HIS 213 45.166 0.002 56.670 1.00 51.63 C ATOM 1456 NE2 HIS 213 44.787 0.157 57.930 1.00 44.85 N ATOM 1457 N CYS 214 50.530 0.634 57.584 1.00 37.10 N ATOM 1458 CA CYS 214 50.986 2.011 57.385 1.00 37.27 C ATOM 1459 C CYS 214 51.504 2.255 55.949 1.00 38.90 C ATOM 1460 O CYS 214 51.699 3.400 55.518 1.00 39.15 O ATOM 1461 CB CYS 214 52.037 2.389 58.443 1.00 36.90 C ATOM 1462 SG CYS 214 53.736 1.793 58.215 1.00 34.05 S

タンパク質はアミノ酸の集合体である。アミノ酸の分子構造はよく知られているが、一部のアミノ酸については水素が結合する位置が異なることがある。例えばヒスチジン(HIS)というアミノ酸では、以下の3種が存在することが知られている(図中の3文字記号はAmberで使われる名称)。


ここでタンパク質に対するMD計算について考えてみる。MD計算ではアミノ酸ごとに計算のためのデータベースを持っており、それを利用して計算される。なので、同じヒスチジンでも水素の付加の位置によって異なるデータが使われる。しかし、PDBデータの多くは水素が付加されておらず(X線結晶構造解析では水素の位置は特定できない)、計算前にアミノ酸の配置や想定されるpHなどから水素を付加させる必要がある。そして、ヒスチジン毎に異なる位置に水素が付加されるわけだ。

異なる位置に水素が付加されたヒスチジンに対して、正しい計算データをMD計算プログラムに伝えるためにそのアミノ酸の名前を変更する必要があるのだが、これが結構な手間となる。例えば、1000個以上のタンパク質構造からヒスチジン分子をビューアで確認してその水素の位置からアミノ酸の名前を付け替えるという作業なんて誰もやりたくないだろう。

そこで、水素の位置の異なるアミノ酸を自動的に判定してアミノ酸の名前を付け直すプログラムを作った。このコードでは部分グラフ同型判定アルゴリズムを利用している。使用したアルゴリズムをここで説明するのは面倒なので、興味のある方はV. Lipets, N. Vanetik, and E. Gudes. Subsea: an efficient heuristic algorithm for subgraph isomorphism. Data Min. Knowl. Disc. 19. 320-350 (2009)あたりを参考にしてほしい。

使い方は次の通り。以下に示す例では、input.pdbが対象とするPDBデータで、名前がHISとなっているアミノ酸について判定する。判定するためのリファレンス構造は3種類で、HID.pdb、HIE.pdb, HIP.pdbの3つのPDB構造に対して行う。構造と一致したら、それぞれアミノ酸名をHID、HIE、HIPに置き換える。最終的に修正した構造をoutput.pdbとして出力する。

% matching_residue input.pdb output.pdb HIS HID HID.pdb HIE HIE.pdb HIP HIP.pdb

標準出力には以下のような表示が出る。この表示の場合、入力PDBファイルには6個のHISがあり、それぞれ17, 66, 132, 146, 167, 261番目のアミノ酸で、HIEとHIPに変換されたことを示している。ここでは判定するための対象としてヒスチジンを例に取っているが、アスパラギン酸リシンでも構わないし、アミノ酸でなくてもいい。

Number of HIS: 6 HIS 17 -> HIE HIS 66 -> HIE HIS 132 -> HIE HIS 146 -> HIP HIS 167 -> HIE HIS 261 -> HIP

今回に処理については、必ずしも部分グラフとして判定する必要はないのだが、今後の拡張のため、例えばアミノ酸の一部から正しい構造を判別する場合などについても考慮して部分グラフとして扱えるようにしている。

最後に今回のプログラムのソースコードを示す。

matching_residue.cpp

// matching_residue.cpp by nox, 2011.04.18 #include <iostream> #include <fstream> #include <sstream> #include <string> #include <vector> #include <map> #include <algorithm> #include <utility> #include <iomanip> #include <limits> #include <cmath> using namespace std; struct Node { vector<Node*> nb; // Node list of neighbor atoms ([Node]) int idx; // index of the atom in the residue (int) int d; // index of traverse history (int) vector<int> t; // neighbors of traverse history ([int]) string atom; // atom type (string) int res_id; // residue id of the PDB (int) vector<double> crd; // coordinate of the atom ([float,float,float]) Node(int idx_, string atom_, int res_id_, vector<double>& crd_) : idx(idx_), d(-1), atom(atom_), res_id(res_id_), crd(crd_) { } }; struct PdbData { string atom; int res_id; vector<double> crd; PdbData(string line) { stringstream ss; double d; atom = line.substr(76, 2); atom.erase(0, atom.find_first_not_of(' ')); atom.erase(atom.find_last_not_of(' ') + 1); ss << line.substr(20, 6); ss >> res_id; ss.str(""); ss.clear(); ss << line.substr(30, 8); ss >> d; crd.push_back(d); ss.str(""); ss.clear(); ss << line.substr(38, 8); ss >> d; crd.push_back(d); ss.str(""); ss.clear(); ss << line.substr(46, 8); ss >> d; crd.push_back(d); ss.str(""); ss.clear(); } }; // comparison function for sorting of traverse history class LessNeighborD { public: bool operator()(const Node* x, const Node* y) { int x_min_d, y_min_d; x_min_d = y_min_d = numeric_limits<int>::max(); for (vector<Node*>::const_iterator p = x->nb.begin(); p != x->nb.end(); ++p) { if ((*p)->d == -1) continue; x_min_d = min(x_min_d, (*p)->d); } for (vector<Node*>::const_iterator p = y->nb.begin(); p != y->nb.end(); ++p) { if ((*p)->d == -1) continue; y_min_d = min(y_min_d, (*p)->d); } return x_min_d < y_min_d; } }; class GreaterNbSize { public: bool operator()(const Node* x, const Node* y) { return x->nb.size() > y->nb.size(); } }; class LessD { public: bool operator()(const Node* x, const Node* y) { return x->d < y->d; } }; class Matching { private: // visit time of traverse history int vtime; // mapping between target and the replacing residue ((Node,Node)) vector<pair<Node*, Node*> > mapping; map<pair<string, string>, double> bond_dist; public: Matching(); double get_bond_dist(const Node* v1, const Node* v2); double distance(const Node* v1, const Node* v2); void connect(vector<Node*>& nodes); void traverse_history(Node* node); void clear_traverse_history(vector<Node*>& nodes); vector<pair<string, vector<string> > > parse_pdb(string inpdb); vector<vector<Node*> > create_nodes_list(vector<pair<string, vector<string> > >& all_lines, vector<string>& rep_res); bool match(Node* g, vector<Node*>& h, int depth); }; Matching::Matching() : vtime(0) { // table of max bond distances bond_dist[make_pair("C", "C")] = 2.0; bond_dist[make_pair("C", "N")] = 1.8; bond_dist[make_pair("C", "O")] = 1.8; bond_dist[make_pair("C", "H")] = 1.2; bond_dist[make_pair("N", "N")] = 1.8; bond_dist[make_pair("N", "O")] = 1.8; bond_dist[make_pair("H", "N")] = 1.2; bond_dist[make_pair("O", "O")] = 1.8; bond_dist[make_pair("H", "O")] = 1.2; } double Matching::get_bond_dist(const Node* v1, const Node* v2) { pair<string, string> atom_pair(make_pair(min(v1->atom, v2->atom), max(v1->atom, v2->atom))); if (bond_dist.find(atom_pair) != bond_dist.end()) return bond_dist[atom_pair]; return 0.0; } double Matching::distance(const Node* v1, const Node* v2) { double d = 0.0; for (int i = 0; i < 3; i++) { double dv = v1->crd[i] - v2->crd[i]; d += dv * dv; } return sqrt(d); } // create connect list for neighbor atoms void Matching::connect(vector<Node*>& nodes) { for (vector<Node*>::iterator p = nodes.begin(); p != nodes.end(); ++p) { for (vector<Node*>::iterator q = nodes.begin(); q != nodes.end(); ++q) { if (*p == *q || (*p)->idx > (*q)->idx) continue; if (distance(*p, *q) < get_bond_dist(*p, *q)) { (*p)->nb.push_back(*q); (*q)->nb.push_back(*p); } } } } // traverse history, similar to canonical labeling void Matching::traverse_history(Node* node) { node->d = vtime; vtime++; vector<Node*> nb; for (vector<Node*>::iterator p = node->nb.begin(); p != node->nb.end(); ++p) { if ((*p)->d == -1) nb.push_back(*p); else node->t.push_back((*p)->d); } sort(node->t.begin(), node->t.end()); sort(nb.begin(), nb.end(), LessNeighborD()); for (vector<Node*>::iterator p = nb.begin(); p != nb.end(); ++p) if ((*p)->d == -1) traverse_history(*p); } void Matching::clear_traverse_history(vector<Node*>& nodes) { for (vector<Node*>::iterator p = nodes.begin(); p != nodes.end(); ++p) { (*p)->d = -1; (*p)->t.clear(); } mapping.clear(); } vector<pair<string, vector<string> > > Matching::parse_pdb(string inpdb) { int res_id, prev_res_id = -1; vector<string> lines; vector<pair<string, vector<string> > > all_lines; string res, line; stringstream ss; fstream fs(inpdb.c_str(), ios_base::in); while (getline(fs, line)) { ss.str(""); ss.clear(); if (line.length() > 20 && line.substr(0, 4) == "ATOM") { ss << line.substr(20, 6); ss >> res_id; } else res_id = -1; if (prev_res_id == res_id) lines.push_back(line); else { if (!lines.empty()) all_lines.push_back(make_pair(res, lines)); lines.clear(); lines.push_back(line); } if (line.length() > 20 && line.substr(0, 4) == "ATOM") { res = line.substr(17, 3); res.erase(0, res.find_first_not_of(' ')); res.erase(res.find_last_not_of(' ') + 1); } else res = ""; prev_res_id = res_id; } all_lines.push_back(make_pair(res, lines)); return all_lines; } vector<vector<Node*> > Matching::create_nodes_list(vector<pair<string, vector<string> > >& all_lines, vector<string>& rep_res) { vector<vector<Node*> > nodes_list; vector<vector<PdbData> > pdb_data; for (vector<pair<string, vector<string> > >::iterator p = all_lines.begin(); p != all_lines.end(); ++p) { string res(p->first); vector<PdbData> pdbs; for (vector<string>::iterator q = p->second.begin(); q != p->second.end(); ++q) { // for l in lines if (q->length() > 78 && q->substr(0, 4) == "ATOM" && find(rep_res.begin(), rep_res.end(), res) != rep_res.end()) pdbs.push_back(PdbData(*q)); } pdb_data.push_back(pdbs); } for (vector<vector<PdbData> >::iterator p = pdb_data.begin(); p != pdb_data.end(); ++p) { if (p->empty()) continue; vector<Node*> nodes; int idx = 0; for (vector<PdbData>::iterator q = p->begin(); q != p->end(); ++q) nodes.push_back(new Node(idx++, q->atom, q->res_id, q->crd)); connect(nodes); vtime = 0; traverse_history(nodes[0]); sort(nodes.begin(), nodes.end(), LessD()); nodes_list.push_back(nodes); } return nodes_list; } bool Matching::match(Node* g, vector<Node*>& h, int depth) { for (vector<Node*>::iterator p = g->nb.begin(); p != g->nb.end(); ++p) if ((*p)->d != -1) g->t.push_back((*p)->d); if (h[depth]->t.size() != g->t.size() || h[depth]->atom != g->atom) { g->t.clear(); return false; } g->d = depth; sort(g->t.begin(), g->t.end()); if (g->t == h[depth]->t) { mapping.push_back(make_pair(h[depth], g)); if (mapping.size() == h.size()) return true; vector<Node*> nb(mapping[h[depth+1]->t[0]].second->nb); sort(nb.begin(), nb.end()); vector<Node*>::iterator end = nb.end(); for (int i = 1; i < h[depth+1]->t.size(); i++) { vector<Node*> nb2(mapping[h[depth+1]->t[i]].second->nb); sort(nb2.begin(), nb2.end()); end = set_intersection(nb.begin(), end, nb2.begin(), nb2.end(), nb.begin()); } nb.erase(end, nb.end()); sort(nb.begin(), nb.end(), GreaterNbSize()); for (vector<Node*>::iterator p = nb.begin(); p != nb.end(); ++p) if ((*p)->d == -1 && match(*p, h, depth + 1)) return true; mapping.pop_back(); } g->d = -1; g->t.clear(); return false; } int main(int argc, char* argv[]) { string inpdb, outpdb, target; vector<pair<string, string> > res_pdb_list; if (argc < 6 || (argc & 1)) { cerr << "Usage: " << argv[0] << " input_pdb output_pdb target_res replace_res1 replace_res1_pdb [replace_res2 replace_res2_pdb...]" << endl; cerr << " ex.) " << argv[0] << " input.pdb output.pdb HIS HID HID.pdb HIE HIE.pdb HIP HIP.pdb" << endl; exit(1); } inpdb = argv[1]; outpdb = argv[2]; target = argv[3]; for (int i = 4; i < argc; i += 2) res_pdb_list.push_back(make_pair(argv[i], argv[i+1])); Matching m; // calculate nodes of replacing residues map<string, vector<Node*> > rep_res; for (vector<pair<string, string> >::iterator p = res_pdb_list.begin(); p != res_pdb_list.end(); ++p) { vector<pair<string, vector<string> > > all_lines(m.parse_pdb(p->second)); vector<string> res; res.push_back(p->first); vector<vector<Node*> > nodes_list(m.create_nodes_list(all_lines, res)); rep_res[p->first] = nodes_list[0]; } // calculate nodes of target residues in the PDB vector<pair<string, vector<string> > > all_lines(m.parse_pdb(inpdb)); vector<string> res; res.push_back(target); vector<vector<Node*> > nodes_list(m.create_nodes_list(all_lines, res)); cout << "Number of " << target << ": " << nodes_list.size() << endl; // matching residues vector<string> result; bool is_found = false; for (vector<vector<Node*> >::iterator p = nodes_list.begin(); p != nodes_list.end(); ++p) { is_found = false; for (map<string, vector<Node*> >::iterator q = rep_res.begin(); q != rep_res.end(); ++q) { if (p->size() != q->second.size()) continue; m.clear_traverse_history(*p); for (vector<Node*>::iterator r = p->begin(); r != p->end(); ++r) { if (m.match(*r, q->second, 0)) { cout << target << " " << setw(4) << (*r)->res_id << " -> " << q->first << endl; result.push_back(q->first); is_found = true; break; } } if (is_found) break; } if (!is_found) { cout << target << " " << setw(4) << (*p)[0]->res_id << " -> Not found." << endl; result.push_back(target); } } // write PDB fstream fs(outpdb.c_str(), ios_base::out); int idx = 0; for (vector<pair<string, vector<string> > >::iterator p = all_lines.begin(); p != all_lines.end(); ++p) { if (p->first == target) { for (vector<string>::iterator q = p->second.begin(); q != p->second.end(); ++q) { size_t pos = q->find(target); if (pos != string::npos) q->replace(pos, target.size(), result[idx]); fs << *q << endl; } idx++; } else { for (vector<string>::iterator q = p->second.begin(); q != p->second.end(); ++q) fs << *q << endl; } } return 0; }

続きを読む...

2011年3月22日火曜日

東北地方太平洋沖地震 可視化地図

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

東日本大震災(東北関東大震災)により被害を受けられた方々には心よりお見舞い申し上げます。

以前、日本全国コンビニ店舗分布地図Processingで作ったが、今回の地震の震源とその大きさについて視覚的に表した地図を作成した。少しでも何かの役に立てば良いのだが。

震源とその大きさを日本地図上に円として表現している。円の大きさはマグニチュードと震源の深さに依存している(これは大まかな目安であり、揺れた地域や震度・エネルギーに正確に対応しているわけではない)。また、震源の深さによって円の色が変わり、浅ければ赤に、深ければ緑に近づく。海の色は時刻を現しており、正午が最も明るい青で表現され、深夜であれば真っ黒となる。データについては、日本気象協会が提供している3月9日から20日までの地震情報を利用させて頂いた。

操作方法は、スペースキーでポーズ、'N'キーで3時間進め、'B'キーで3時間戻る。'+'キーで時間の進みが速くなり、'-'キーで遅くなる。マウスの左クリックでドラッグすれば地図を動かすことができ、右クリックで上下にドラッグすれば拡大・縮小となる。



可視化地図から、3月11日14時46分のマグニチュード9.0の地震の起こる二日ほど前から、三陸沖で中規模の地震が頻発していることが分かる。ところが、M9.0の地震の数時間前からは揺れが全て途絶え、嵐の前の静けさの様相となる。そして、M9.0の地震の直後からはひっきりなしに揺れが続き、数日してやや沈静化してきている。ただ、M6程度の地震は未だに起こっているので、安心するにはまだ早いかもしれない。最近の地震では福島県沖や茨城県沖に多く、ややと南下している様子が分かる。

以下、ソースコードを示す。

earthquake.pde

/* earthquake.pde by nox, 2011.3.22 */ String convData = "quake.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 offsetX = 0.0; float offsetY = 0.0; int centerX = 0; int centerY = 0; float zoom = 1.0; int mxStart, myStart, myPos; int startTime, currentTime, stopTime; int month_, day_, hour_, minute_; boolean is_paused = false; int speed = 10; public void setup() { mapShape = loadShape(japanMap); smooth(); loop(); size(int(mapShape.width), int(mapShape.height), JAVA2D); // 539, 563 shapeMode(CORNER); mapShape.disableStyle(); font = loadFont("CourierNewPSMT-14.vlw"); textMode(SCREEN); textFont(font); readData(); startTime = millis(); loop(); } public void draw() { if (!mousePressed) currentTime = millis() - startTime; int currentDate = currentTime + (hour_ * 60 + minute_) * speed; int mon_ = month_; int d_ = day_ + currentDate / (24 * 60 * speed); int h_ = currentDate / (60 * speed) % 24; int m_ = currentDate / speed % 60; if (d_ >= 21 || d_ <= 8 || h_ < 0 || m_ < 0) startTime = millis(); background(color(0, 0, (720 - abs(h_ * 60 + m_ - 720)) * 0.1)); noStroke(); fill(32); shape(mapShape, int(offsetX * zoom) + centerX, int(offsetY * zoom) + centerY); fill(127); String timeFormat = nf(mon_, 2) + "/" + nf(d_, 2) + " " + nf(h_, 2) + ":" + nf(m_, 2) + " speed: " + nf(10.0 / speed, 1, 2); text(timeFormat, 15, 20); for (int i = 0; i < placeCount; i++) places[i].draw(); } 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() { if (!is_paused) { mxStart = int(mouseX - offsetX * zoom); myStart = int(mouseY - offsetY * zoom); myPos = mouseY; stopTime(); loop(); } } void mouseReleased() { if (!is_paused) startTime(); } void mouseDragged() { if (is_paused) return; if (mouseButton == LEFT) { offsetX = (mouseX - mxStart) / zoom; offsetY = (mouseY - myStart) / zoom; } else if (mouseButton == RIGHT) { if (mouseY - myPos > 4) { float zoomOut = 9.8 / 10.0; mapShape.scale(zoomOut); zoom *= zoomOut; centerX = int(width * (1.0 - zoom) * 0.5); centerY = int(height * (1.0 - zoom) * 0.5); myPos = mouseY; } else if (mouseY - myPos < -4) { float zoom_in = 10.0 / 9.8; mapShape.scale(zoom_in); zoom *= zoom_in; centerX = int(width * (1.0 - zoom) * 0.5); centerY = int(height * (1.0 - zoom) * 0.5); myPos = mouseY; } } } void startTime() { startTime += millis() - stopTime; loop(); } void stopTime() { noLoop(); stopTime = millis(); } void keyPressed() { if (key == ' ') { if (is_paused) { startTime(); is_paused = false; } else { stopTime(); is_paused = true; } } else if (key == 'n' || key == 'N') startTime -= 180 * speed; else if (key == 'b' || key == 'B') startTime += 180 * speed; else if (key == '+' && speed > 3) speed -= 1; else if (key == '-' && speed < 30) speed += 1; } void readData() { new Slurper(); } Place parsePlace(String line) { String pieces[] = split(line, '\t'); int mins = int(pieces[0]); int mon = int(pieces[1]); int d = int(pieces[2]); int h = int(pieces[3]); int m = int(pieces[4]); float x = float(pieces[5]); float y = float(pieces[6]); float magnitude = float(pieces[7]); float depth = float(pieces[8]); String name = pieces[9]; return new Place(mins, mon, d, h, m, x, y, magnitude, depth, name); }

Place.pde

class Place { int mins; int mon, d, h, m; float x, y; float magnitude; float depth; String name; public Place(int mins, int mon, int d, int h, int m, float x, float y, float magnitude, float depth, String name) { this.mins = mins; this.mon = mon; this.d = d; this.h = h; this.m = m; this.x = x; this.y = y; this.magnitude = magnitude; this.depth = depth; this.name = name; } public void draw() { int xx = (int)TX(this.x) + int(offsetX * zoom); int yy = (int)TY(this.y) + int(offsetY * zoom); int a = mins * speed + 255 - currentTime; if (magnitude == 0.0 || a < 0 || a > 255) return; float r = pow(5.0, (magnitude / 2.7)) / sqrt(depth / 10 + 1); float m = r * ((255 - a) / 50.0) * zoom; fill(color(255 - this.depth * 2, this.depth * 2, 0), a); ellipse(xx, yy, m, m); } }

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]; line = reader.readLine(); places[placeCount++] = parsePlace(line); month_ = places[0].mon; day_ = places[0].d; hour_ = places[0].h; minute_ = places[0].m; while ((line = reader.readLine()) != null) places[placeCount++] = parsePlace(line); } catch (IOException e) { e.printStackTrace(); } } }

続きを読む...

2011年3月7日月曜日

開発エンジニアの世界へ

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

昨年は開発エンジニアが転職するという話題が多かったように思う。それに触発されたというわけではないが、自分も昨年11月に開発エンジニアとして転職した。そしてつい先日、試用期間の3ヶ月が過ぎた。これで何かが変わったというわけではないが、正式に採用されたということで、転職して思うことなどをブログに書くことにした。転職先については敢えて名前を出さないが、主にマルチコアやGPGPUなどの並列処理を扱うベンチャー企業といえばピンとくる人もいるかも知れない。

自分はこれまでずっとアカデミックな研究所に所属していたのだけど、少々居座りすぎてしまったことは否めない。基本的に研究者は流動的であることが推奨される。そして、所属していた研究所は所長でさえも任期制で、研究員は普段から職探しをするのが当たり前になっている。しかし、アカデミックなポストは空きも少なく、枠が一人のところに100人200人と応募してくることはザラで、終身雇用の職を得るのはかなり難しい。

そんな中、自分はほとんど職探しをしていなかった。第一の理由として、自分が専攻した薬学の分野でバリバリとコードを書くことと研究が両立できるか疑問だったことがある。ただでさえ狭き枠にそのような特殊な環境を求めるのが非常に厳しいものに思えたからだ。第二の理由として、所属していた研究所ではメインでHPCを扱っており薬学に関連した研究を遂行するためのプログラムを自由に書くことができ、それが求められていたことがある。相手の要求と自分の希望が合致していた。

そうは言っても、ここ最近の政治的判断による不透明さや年齢に伴うマネージメントの必要性など、当たり前だが要求は厳しくなってくる。また、地理的な面でも大きく変わることになったし、こちらの事情で現在の住処を離れたくないということもあった。研究所で一緒に仕事をしていた人たちには大変お世話になっていたので少々迷うところもあったのだが、最終的に転職することにした。

自分の仕事選びの基準は単純で「楽しめるかどうか」だけだ。楽しめればどんなにきつい仕事でも苦にならないし、逆にきつければきついほどチャレンジングであり、やりがいを感じられる。特に開発エンジニアの方なら、自分の技術で困難な問題を克服する達成感、創造する喜び、知識を得る楽しみなどについて分かってもらえるのではないかと思う。それに、企業としても生産性が上がるわけだから、まさにWin-Winの関係だ。もし、エンジニアでこれらの喜びや楽しみを見いだせない人がいたとしたら、何故、極端な専門職であるエンジニアを選んだのか知りたいところだ。

今の職場は自分自身の意志による最初のエントリーだったが、有難い事にそこで働けることになった。何年も前から知っていて、とても技術力のある企業だと思ったからここを選んだ。そして、そう思ったことは間違っておらず、想像していたより遥かに楽しんでいる自分がいる。大企業とは違って小回りはきくし、比較的やりたい事をやらせてもらえるし(もちろん責任は伴う)、別プロジェクトの開発者とは同じフロアにいるので交流し易い。そしてなにより、この職場ではプログラミングコンテスト優勝者やICPC関係者、ネットで時々名前を見かけるような技術力の高い人たちがいて話を聞くだけでも楽しい。示し合わせたわけでもないのに自分を含めて数人、先日のBoost.勉強会に参加している事実を以てしても、技術に対して真摯であることが伺える(参加はしなかったがライブ配信を見ていたという人たちもいたようだ)。それに、個人だと到底いじることができないような特殊な環境に触れることができるのも良い(まあ、特殊な環境という意味では研究所にいたころも複数のスーパーコンピュータをいじっていたわけで、それはそれで楽しかったけど)。ついでに、報酬や待遇も悪くないと思う。

最近では若い人たちを積極的に採用しているようで会社説明会も開いている。プログラミング(特に並列処理)に興味を持っている方に参加してもらえると嬉しい。

続きを読む...

2011年2月15日火曜日

複数のGPUでマンデルブロ集合を並列計算

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

範囲(-0.005, -0.005, 0.005, 0.005)のマンデルブロ集合を描画サイズ4,800x4,800、繰り返しの上限10,000として計算させてみた。C++でコードを書いてXeon X5650 1コアで走らせたところ942.97秒かかった(SIMD最適化はしていない)。TBBを利用した12論理コアの並列計算では88.68秒だった。次いでCUDAを使ってTeslaで走らせてみたが7.45秒まで短縮された。Teslaパない。そして贅沢にも2枚のTeslaをスレッドで並列化して使ってみたら4.26秒で計算できた。実にXeon 1コアの220倍の速度が出たわけだ。カリカリチューンをしなくてもこの程度の高速化ができるということが重要だ。

因みにマンデルブロ集合の上記の範囲は数値が発散しない集合部分であり、描画させてみても真っ黒なので面白くはない。これは並列化した際に計算が偏らないようにしたかったのでこのような範囲を指定している。プログラムで出力させたファイルは色の生データなのでそのままでは表示できない。一応、画像表示するためのPythonコードも付け加えておく。

以下に、GPU2枚を利用してマンデルブロ集合を計算するCUDAコードを示しておく。

mandelbrot_thread.cu

// Mandelbrot set using GPGPU by nox, 2011.02.12 // nvcc -lcutil_x86_64 -arch sm_13 -use_fast_math -prec-sqrt=false -keep -L ~/NVIDIA_GPU_Computing_SDK/C/lib -I ~/NVIDIA_GPU_Computing_SDK/C/common/inc -g mandelbrot_thread.cu -o mandelbrot_thread #include <iostream> #include <fstream> #include <cutil_inline.h> #include <multithreading.h> using namespace std; const int BLOCK_SIZE_X = 16; const int BLOCK_SIZE_Y = 16; struct MandelbrotParam { int dev; // GPU device id uchar4* h_color; // image buffer float t, l, w, h; // top, left, width, height int sw, sh; // screen width, screen height int max_iter; // max number of iterators float th; // threshold }; __device__ uchar4 coloring(int n, int b) { int d = n % b; d *= 256 / b; int m = static_cast<int>(d / 42.667f); switch (m) { case 0: return make_uchar4(0, 6 * d, 255, 0); case 1: return make_uchar4(0, 255, 255 - 6 * (d - 43), 0); case 2: return make_uchar4(6 * (d - 86), 255, 0, 0); case 3: return make_uchar4(255, 255 - 6 * (d - 129), 0, 0); case 4: return make_uchar4(255, 0, 6 * (d - 171), 0); case 5: return make_uchar4(255 - 6 * (d - 214), 0, 255, 0); default: return make_uchar4(0, 0, 0, 0); } } __global__ void mandelbrot(float t, float l, float w, float h, int sw, int sh, int max_iter, float th, uchar4* d_color) { int ix = blockIdx.x * blockDim.x + threadIdx.x; int iy = blockIdx.y * blockDim.y + threadIdx.y; if (ix >= sw || iy >= sh) return; float ci = t + (static_cast<float>(iy) / sh) * h; float cr = l + (static_cast<float>(ix) / sw) * w; float zi = 0.0f; float zr = 0.0f; float zrzi, zr2, zi2; for (int i = 0; i < max_iter; i++) { zrzi = zr * zi; zr2 = zr * zr; zi2 = zi * zi; zr = zr2 - zi2 + cr; zi = zrzi + zrzi + ci; if (zi2 + zr2 >= th) { d_color[ix*sh+iy] = coloring(i, 256); return; } } d_color[ix*sh+iy] = make_uchar4(0, 0, 0, 0); } CUT_THREADPROC mandelbrot_thread(void* args) { MandelbrotParam* param = (MandelbrotParam*)args; cudaSetDevice(param->dev); uchar4* h_color = param->h_color; float t = param->t; float l = param->l; float w = param->w; float h = param->h; int sw = param->sw; int sh = param->sh; int max_iter = param->max_iter; float th = param->th; dim3 num_blocks((sw - 1) / BLOCK_SIZE_X + 1, (sh - 1) / BLOCK_SIZE_Y + 1, 1); dim3 num_threads(BLOCK_SIZE_X, BLOCK_SIZE_Y, 1); uchar4* d_color; cudaMalloc((void**)&d_color, sizeof(uchar4) * sw * sh); mandelbrot<<<num_blocks, num_threads>>>(t, l, w, h, sw, sh, max_iter, th, d_color); cudaThreadSynchronize(); cudaMemcpy(h_color, d_color, sizeof(uchar4) * sw * sh, cudaMemcpyDeviceToHost); cudaFree(d_color); CUT_THREADEND; } // outfile : output filename // t : top position // l : left position // w : width // h : height // sw : screen width // sh : screen height // max_iter : max number of iterations [256] // th : threshold [2.0] void write_mandelbrot(const string outfile, float t, float l, float w, float h, int sw, int sh, int max_iter=256, float th=2.0f) { MandelbrotParam param[2]; uchar4* h_color; cudaMallocHost((void**)&h_color, sizeof(uchar4) * sw * sh); // using 2 GPUs for (int i = 0; i < 2; i++) { param[i].dev = i; param[i].h_color = h_color; param[i].t = t; param[i].l = l; param[i].w = w / 2.0f; param[i].h = h; param[i].sw = sw / 2; param[i].sh = sh; param[i].max_iter = max_iter; param[i].th = th; } param[1].h_color = h_color + sw * sh / 2; param[1].l = l + w / 2.0f; unsigned int timer; cutCreateTimer(&timer); cutStartTimer(timer); CUTThread threads[2]; for (int i = 0; i < 2; i++) threads[i] = cutStartThread((CUT_THREADROUTINE)mandelbrot_thread, (void*)&param[i]); cutWaitForThreads(threads, 2); cutStopTimer(timer); cout << "Time: " << cutGetTimerValue(timer) / 1000 << endl; fstream fs(outfile.c_str(), ios_base::out); fs << sw << endl; fs << sh << endl; for (int i = 0; i < sw * sh; i++) fs << h_color[i].x << h_color[i].y << h_color[i].z; fs.close(); cudaFreeHost(h_color); } int main(int argc, char* argv[]) { string outfile("mandelbrot.dat"); if (argc >= 2) outfile = argv[1]; // write_mandelbrot(output filename, // top position, // left position, // width, // height, // screen width, // screen height, // max number of iterations [256] // threshold [2.0]); //write_mandelbrot(outfile, -1.0f, -1.5f, 2.0f, 2.0f, 480, 480, 100); //write_mandelbrot(outfile, 0.57f, 0.34f, 0.036f, 0.027f, 6400, 4800, 1000); write_mandelbrot(outfile, -0.005f, -0.005f, 0.01f, 0.01f, 4800, 4800, 10000); return 0; }

表示用のPythonコード。

read_mandelbrot.py

#!/usr/bin/env python import sys, os, Image def main(args): if len(args) < 2: print >>sys.stderr, "%s datafile [imagefile]" % os.path.basename(args[0]) sys.exit() datafile = args[1] imagefile = args[2] if len(args) > 2 else "" f = file(datafile) sw = int(f.readline()) sh = int(f.readline()) data = map(ord, f.read()) im = Image.new("RGB", (sw, sh)) for x in range(sw): for y in range(sh): pos = (x * sh + y) * 3; im.putpixel((x, sh - y - 1), tuple(data[pos:pos+3])) if imagefile: im.save(imagefile) im.show() if __name__ == "__main__": main(sys.argv)

続きを読む...