スキップしてメイン コンテンツに移動

投稿

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

今回は興味を持ってもらえる人が限定される記事だと思うけど、以前作ったコードを腐らせてしまうのも勿体無いので敢えて晒してみる。 生体内のタンパク質構造を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 ...

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

東日本大震災(東北関東大震災)により被害を受けられた方々には心よりお見舞い申し上げます。 以前、 日本全国コンビニ店舗分布地図 を Processing で作ったが、今回の地震の震源とその大きさについて視覚的に表した地図を作成した。少しでも何かの役に立てば良いのだが。 震源とその大きさを 日本地図 上に円として表現している。円の大きさはマグニチュードと震源の深さに依存している(これは大まかな目安であり、揺れた地域や震度・エネルギーに正確に対応しているわけではない)。また、震源の深さによって円の色が変わり、浅ければ赤に、深ければ緑に近づく。海の色は時刻を現しており、正午が最も明るい青で表現され、深夜であれば真っ黒となる。データについては、 日本気象協会 が提供している3月9日から20日までの 地震情報 を利用させて頂いた。 操作方法は、スペースキーでポーズ、'N'キーで3時間進め、'B'キーで3時間戻る。'+'キーで時間の進みが速くなり、'-'キーで遅くなる。マウスの左クリックでドラッグすれば地図を動かすことができ、右クリックで上下にドラッグすれば拡大・縮小となる。 This browser does not have a Java Plug-in. Get the latest Java Plug-in here. 可視化地図から、3月11日14時46分のマグニチュード9.0の地震の起こる二日ほど前から、三陸沖で中規模の地震が頻発していることが分かる。ところが、M9.0の地震の数時間前からは揺れが全て途絶え、嵐の前の静けさの様相となる。そして、M9.0の地震の直後からはひっきりなしに揺れが続き、数日してやや沈静化してきている。ただ、M6程度の地震は未だに起こっているので、安心するにはまだ早いかもしれない。最近の地震では福島県沖や茨城県沖に多く、ややと南下している様子が分かる。 以下、ソースコードを示す。 earthquake.pde /* earthquake.pde by nox, 2011.3.22 */ String convData ...

開発エンジニアの世界へ

昨年は開発エンジニアが転職するという話題が多かったように思う。それに触発されたというわけではないが、自分も昨年11月に開発エンジニアとして転職した。そしてつい先日、試用期間の3ヶ月が過ぎた。これで何かが変わったというわけではないが、正式に採用されたということで、転職して思うことなどをブログに書くことにした。転職先については敢えて名前を出さないが、主にマルチコアやGPGPUなどの並列処理を扱うベンチャー企業といえばピンとくる人もいるかも知れない。 自分はこれまでずっとアカデミックな研究所に所属していたのだけど、少々居座りすぎてしまったことは否めない。基本的に研究者は流動的であることが推奨される。そして、所属していた研究所は所長でさえも任期制で、研究員は普段から職探しをするのが当たり前になっている。しかし、アカデミックなポストは空きも少なく、枠が一人のところに100人200人と応募してくることはザラで、終身雇用の職を得るのはかなり難しい。 そんな中、自分はほとんど職探しをしていなかった。第一の理由として、自分が専攻した薬学の分野でバリバリとコードを書くことと研究が両立できるか疑問だったことがある。ただでさえ狭き枠にそのような特殊な環境を求めるのが非常に厳しいものに思えたからだ。第二の理由として、所属していた研究所ではメインでHPCを扱っており薬学に関連した研究を遂行するためのプログラムを自由に書くことができ、それが求められていたことがある。相手の要求と自分の希望が合致していた。 そうは言っても、ここ最近の政治的判断による不透明さや年齢に伴うマネージメントの必要性など、当たり前だが要求は厳しくなってくる。また、地理的な面でも大きく変わることになったし、こちらの事情で現在の住処を離れたくないということもあった。研究所で一緒に仕事をしていた人たちには大変お世話になっていたので少々迷うところもあったのだが、最終的に転職することにした。 自分の仕事選びの基準は単純で「楽しめるかどうか」だけだ。楽しめればどんなにきつい仕事でも苦にならないし、逆にきつければきついほどチャレンジングであり、やりがいを感じられる。特に開発エンジニアの方なら、自分の技術で困難な問題を克服する達成感、創造する喜び、知識を得る楽しみなどについて分かってもらえるのではないかと思う。それに、企業と...

複数の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_SIZ...

2010年から2011年へ

明けましておめでとうございます。 今年は2011年で素数。前回の素数の年は8年前の2003年、次の素数の年は6年後の2017年。だからどうと言うわけでもないのだけど、なんとなく特別な年であるように思えてくる。昨年は自分にとっては大きな節目の年であった。環境も大きく変わって、なんちゃってプログラマから「なんちゃって」を外すことができたのではないかと思う。そして今年は、どんな数にも割り込まれず自分自身を失わない2011年のように、自らの意志をしっかり持って行動する年にしたい。 ところで、冒頭のQRコードだけど、いつものように Android SL4A のソースコードになっている。SL4Aを起動して、メニューからAddを選択し、Scan Barcodeで取り込むと new_year.py が作られる。大したコードではないのだけど、当ブログからのメッセージが入っている。音が出るので注意。 これからもどうぞよろしくお願いします。

ある重さを量るのに使う錘の組み合わせは何通り?

問題 : 1, 5, 10, 50, 100, 500グラムの6種類の錘があります。それぞれの錘はいくつでも使うことができます。ある重さを量るときに使う錘の組み合わせは何通りになりますか。それを求めるプログラムを書いてください。ここで量る重さはグラム単位の整数とし、最大で100キログラムとします。10グラムの重さを量る場合、10グラムの錘が1個、5グラムの錘が2個、5グラムの錘1個と1グラムの錘が5個、1グラムの錘が10個の組み合わせになり、全部で4通りとなります。 以下、解答例。 weights.cpp // 指定した重さを量るのに使う錘の組み合わせの数を求めるプログラム. // 錘は 1g, 5g, 10g, 50g, 100g, 500g の6種類. // 重さ 700,000g 程度までなら桁あふれせずに計算できる. #include <iostream> #include <map> #include <algorithm> #include <cstdlib> using namespace std; const int weights[] = { 1, 5, 10, 50, 100, 500 }; // 錘(昇順). const int N = sizeof(weights) / sizeof(int); // 錘の種類. // 求めた錘の組み合わせの数のメモ. // 下位3ビットは使用する錘(solveのidx)を表す. // 残りのビットは重さ(solveのv)を表す. map<int, unsigned long long> memo; // 指定した重さに対する錘の組み合わせの数を求める. // v : 重さ // idx : 使用する錘 // N-1の場合はweights[0]からweights[N-1]までの錘を使用する(すべての錘). // 2の場合はweights[0]からweights[2]までの錘を使用する(1, 5, 10の錘). // 戻り値: 組み合わせの数 unsigned long long solve(int v, int idx=N-1) { // 既に...

Android端末による自動通訳を一行のPythonコードで実装する

最近このブログではAndroidの SL4A ネタが多いのだけど、それだけ遊べるのだから仕方がない。SL4AのPythonを使えば実質一行で自動通訳プログラムを作ることもできる。日本語で喋った文章が翻訳されて英語の音声で返ってくる。 import android,urllib,urllib2,simplejson;droid=android.Android();droid.ttsSpeak(simplejson.loads(urllib2.urlopen('http://ajax.googleapis.com/ajax/services/language/translate?v=1.0&q=%s&langpair=%s%%7C%s'%(urllib.quote(droid.recognizeSpeech().result.encode('utf-8')),'ja','en')).read())['responseData']['translatedText']) このプログラムでは日本語を英語に変換するが、コード中の en を変更すれば別の外国語に通訳することもできる。例えば fr に変更すればフランス語に通訳される。ただし、通訳した外国語はGoogleによる翻訳の精度に依存するので、実用的に用いるにはもう少し精度が良くなって欲しいところだ。それでもこれだけ簡単に自動通訳プログラムを作れるのだからAndroidとSL4Aは大したものだと思う。 以下にもう少し分かりやすいコードを示しておく。ついでに翻訳結果をダイアログに表示するように変更してある。 translator.py # -*- coding: utf-8 -*- import android,urllib,urllib2,simplejson def translator(text,from_lang,to_lang): url='http://ajax.googleapis.com/ajax/services/language/translate?v=1.0&q=%s&langpair=%s%%7C%s'%(urllib.quote(...