matplotlib の散布図でグラデーションをかける方法

よく使われている Python の数式処理用ライブラリに NumPy, SciPy があります。そして、これらを使って得た結果をプロットするのに matplotlib が利用されます。これらはどれもフリーかつ多機能ので色々便利です（しかも Python らしからぬ速さ）。グラフを描画するのには gnuplot を使う人も多いと思いますが、今回は matplotlib に話を限ります。

統計をやっていたりすると散布図を描く必要が出てくるのですが、その中でも特にデータの順番をグラフィカルに表すため各データ点の色にグラデーションをかけてプロットしたくなることがあります。matplotlib を使ってこれを達成するにはどうすればいいか、というのが本エントリのお題です。まずは以下のスクリプトをご覧ください。

# colordemo.py by hidehideta
# -*- coding: utf-8 -*-

from pylab import *
import numpy as np

N = 100
theta =  2 * np.pi * np.linspace(0, 1, N)
figure(figsize=(6, 6))
xlabel("x"); ylabel("y")

x1 = 0.5 * np.cos(theta) ** 3
y1 = 0.5 * np.sin(theta) ** 3
color1 = np.linspace(0, 1, N)
scatter(x1, y1, c=color1)

x2 = 2 * np.cos(theta) ** 3
y2 = 2 * np.sin(theta) ** 3
M = (np.array([1, -1, 1, 1]) / np.sqrt(2)).reshape((2, 2))
x2, y2 = M.dot([x2, y2])

color2 = np.zeros((N, 4))
#color2[::, 0] = [i / float(N) for i in range(N)] # R
color2[::, 1] = [i / float(N) for i in range(N)] # G
color2[::, 2] = [i / float(N) for i in range(N)] # B
color2[::, 3] = [1 for i in range(N)] # A: 0...1 <-> 透明...不透明
scatter(x2, y2, c=color2)

savefig("colordemo.png")

scatter() は c= の形で色を指定する文字もしくは array を受けとることができます。特に array を与える場合はプロットするデータと同じ長さにしておく必要があって、しかも各要素の値は 0 から 1 までの範囲に限られています（255 とか入れるとエラーになる）。
上のスクリプトでいうと color1 が 0 から 1 までの N 点であり、図の内側のアステロイドに対応します。この場合は青から赤に向かって自動的にグラデーションをかけてくれます。
color2 は (N, 4) 行列で、各行が RGBA を指定しています。つまり、長さが N で各要素が RGBA の4成分に対応するリストということです。上の例だと G と B を同時に 0 から 1 まで増やし、かつ α 値を 1 に保っているので徐々にシアンが鮮やかになっていきます（こちらは反時計回りに 45 度回転したアステロイドで、右上の頂点がスタート）。

特にグラデーションにこだわりがないのなら上の color1 を使えばいいでしょうし、好みの配色にしたいのなら color2 を使って RGBA を自分で決めれば OK です。同じく頻繁に使う plot メソッドで color=array とするとエラーが出るので、またこっちの場合についても分かったら記事にするかもしれません。

参考 URL: http://matplotlib.sourceforge.net/api/pyplot_api.html#matplotlib.pyplot.scatter

グラフの軸の色だけを変更する Python スクリプト

お久しぶりです。目標を掲げて毎日着実に進もうという決意とともに新年度を迎えております。とりあえず TOEIC 800 点取るのと、C++ への習熟度を高めることを研究以外のがんばりどころとしたいです。

さて突然ですが、研究の成果などをグラフにして表示するとき、軸の色を黒から白に変更したいと思ったことはないでしょうか？僕は研究室のミーティングで専ら黒背景のスライドを使ってトークするのですが、こういうときは図の背景を透明にして軸や凡例文字を白くし、グラフ自体は何か意味のある色を使って表示するのが素敵だとずっと思ってきました。このエントリではこれを叶えてくれる Python スクリプトの説明をします。

右図は gnuplot で適当に出力した sin 関数で、元々は eps 形式だったのを convert -background none -density 300 sin.eps sin.png として png 形式にしたものです。-background オプションをつけることで背景を透明にしています。今やろうとしているのは、赤い線はそのままにした状態で軸と sin(x) という凡例の文字だけを白くすることです。これは PIL を使うと以下に示す通り簡単に出来ます。注意点は、PIL のバージョン 1.1.6 以上を使うこと。詳しくはこのページの load メソッドについてのコメントを見てください。

# b2w.py by hidehideta

from sys import argv
from os import wait
import Image

def black2white(ori):
    if ori.endswith(".eps"):
        import subprocess
        command = "convert -background none -density 300 " \
                  + ori + " " + ori[:-4] + ".png"
        subprocess.Popen(command.split()); wait()
        ori = ori[:-4] + ".png"

    UP = 20 # threshold value to judge "black" pixels
    im = Image.open(ori)
    h, w = im.size
    jm = Image.new(im.mode, im.size)

    # pixel access objects
    pi = im.load()
    pj = jm.load()

    for i in range(h):
        for j in range(w):
            p = pi[i, j]
            if p[0] < UP and p[1] < UP and p[2] < UP:
                pj[i, j] = (255, 255, 255, 255)
            else: pj[i, j] = p

    modified = ori[:-4] + "_b2w.png"
    jm.save(modified, format="png")

if __name__ == "__main__":
    if len(argv) != 2:
        print "Usage: python b2w.py < image file to convert >"
    else:
        ori = argv[1]
        black2white(ori)

ポイントは 21, 22 行目で、load() メソッドを呼ぶことで各ピクセルの値を取得、変更できるようにしています（図は RGBA モードで読み込まれているので各ピクセルごとに４つの要素を持つタプルが割り当てられている）。肝心の黒色は RGB で (0, 0, 0) に対応するので、UP で決まる値程度までは遊びを持たせて 27 行目で判別しています。以下、軸や凡例に対応すると思われるピクセルには白を表す (255, 255, 255, 255) を与え、他はオリジナルのピクセルをコピーして終了。最終的に右図が得られるはずです。入力する画像のフォーマットは PIL がサポートしているものなら何でもよいと思いますが、png 以外試してないので断言できません。なお、eps ファイルとか png ファイルがたくさんあるのを一括処理するには glob.glob() を使うとすごく楽です。

今回はかなり原始的な方法を使いましたが、Scipy や OpenCV を利用したらもっと簡単かつ複雑なことも出来るのかもしれません（もしかして ImageMagick だけで何とかなるのだろうか）。また何か分かったら書きます。

ちなみにこの記事が Syntax Highlighter 初挑戦だったんですけど、これ、どうだろう？幅の細かい調整とかどうやるのか分からなくて困る。。。（今回はこことかのおかげで何とか使えるようにはなりました）

生産性を上げる魔法の道具

Apple Keyboard を購入しました。日本語、テンキー付きでコンピュータと USB 接続するタイプです。もちろん Mac 用に使うのが普通のところを今回は Ubuntu 10.04 で試してみました。

システム -> 設定 -> キーボード -> レイアウト
とたどって、キーボードの形式を"Apple Aluminium Keyboard (JIS)"として

• Alt/Window キーの挙動： Control is mapped to Alt keys, Alt is mapped to Win keys
• CapsLock キーの挙動： Make CapsLock an additional ESC
• Ctrl キーの位置： Make CapsLock and additional  Ctrl

と設定しました。小指ならどちらでも Ctrl キーが使える、逆に CapsLock を使わないようにする、という感じの設定です。どうやら fn キーは使えないようですが大した問題じゃないとは思います。

研究室では一年近く Logicool の製品 を使っていたのですが、自室の Apple Keyboard （こっちは US タイプ）に比べると明らかに鍵打が遅かったのでちょっと思い切って浮気してみました。テンキー付きのこの製品は bluetooth 接続するタイプのものに比べてかちゃかちゃ感が抑えられていて、キーを叩いたときの力がより無駄なく伝わっているような気がします。もちろん Logicool のやつより速く打てます。ストロークが違うんだよね。たぶん。

USB 端子が側面に付属しているのも嬉しいですね。どうも USB メモリは認識されないようですが….…。

ところで、本屋で立ち読みして実に面白いと思ったので『Python クックブック』（第2版）を購入しました。Python の言語仕様に関連する項目からテクニカルに役立つものまで、実に多彩なレシピが載っています。訳は原著者の妙なハイテンションすら伝わってくるような自然さです。単にマニュアル的な本というわけではなく、目下の課題に対して解法を与え、どうすれば Pythonic な書き方が出来るのかごくごく丁寧に説明してくれます。Python の基本的な使い方を知っている人なら何がしかの役立つ情報が得られると思います。第2版の段階では Python 2.4 までを想定して書かれているので、2.5 以降とは状況が違っているような部分もあるのかもしれません。原書はすでに第3版が出ており、2.6 から 3.x までに対応しているようです。

まだざっと目を通した段階ですが、すでにいくつかのレシピを今書いているプログラムに反映させました：os.popen() なんて関数が存在するとは知らなかったです（恥ずかしながら）。これまでは C で popen() 使って gnuplot を呼ぶ関数を書いた後 DLL を作って ctypes.cdll で読み込んでいたのですが、os モジュールによっていくつもの手順が省けるようになりました。まあしょぼい例ではありますが、『クックブック』をパラパラと眺めていたときに見つけることが出来て良かったです。今後もお世話になる予感。

Cython 素敵

Cython なるものを使い始めました。これはちょっと一言で表し切るのが難しいのですが、Python とほとんど同じ文法を使いつつも C 言語並みの実行速度を得る事が出来るようなプログラムのことです。実物を見てもらった方が早いと思うので、今日から取り組み始めた Project Euler の第９問目を題材にしてごく簡単に説明します。

問題はこうです：和がちょうど 1000 になるようなピタゴラスの三つ組の積を求めよ。

ひとまず、何も考えずに下のようなスクリプトを書きました。

#pe9_job_normal.py
def job():
    for i in range(1, 1001):
        for j in range(i+1, 1001):
            for k in range(j+1, 1001):
                if i**2 + j**2 == k**2 and i + j + k == 1000:
                    return i*j*k

#pe9_normal.py
import pe9_job_normal
print pe9_job_normal.job()

time python pe9_normal.py とした結果は
real    0m36.493s
user    0m36.172s
sys    0m0.046s
で、体感としてはちょっと待たされたという感じでした。

肝心の Cython を利用したコードは次の通りです。

#pe9_job_cython.pyx
def job():
    cdef int i, j, k
    for i in range(1, 1001):
        for j in range(i+1, 1001):
            for k in range(j+1, 1001):
                if i**2 + j**2 == k**2 and i + j + k == 1000:
                    return i*j*k

#pe9_cython.py
import pyximport; pyximport.install()
import pe9_job_cython
print pe9_job_cython.job()

元のスクリプトとの相違点は３カ所あります。１つはループ計算をする関数が .py ではなく .pyx ファイルに書かれていること。そのファイルの 2 行目では cdef int という C 言語ライクな型宣言がなされています。この .pyx ファイルを import pyximport することで pe9_cython.py にインポートしています。ただこれだけのことですが、実行結果は
real    0m0.381s
user    0m0.255s
sys    0m0.074s
となり、およそ 100 倍の速度を得ています。ループ計算中に頻繁に用いられる変数の型を決めてしまう事で、実行速度が劇的に改善される事が分かりました（むしろ Python はどんだけ遅いんだと思ってしまうくらい）。range 関数を使う場合は上のコードのようにループ変数を cdef することで最適化が行われるそうです。

Cython のドキュメントなどによると元々は Sage を開発する目的で作られたとのことで、確かに Python みたいな書きやすさと C に匹敵する速さがあれば開発スピードも実行速度も申し分ないのだろうと想像は出来ます。表面上は C に似ているとはいえ、コンパイルしなくていいというだけでも心理的な負担はかなり少ないです。

Cython については今日気まぐれで使い始めただけなので、まだ色々分かってません。Python のリストや辞書をどう扱ったらいいのかとか、pure C との使い分けをどうするかとか理解したいです。今のところは Python スクリプトを気楽に高速化するための便利システムという程度の認識ですが、今後はもっと入れ込んだりするのかも知れません。なんか分かったらまた書きます。

Python 関連の記事ばっかりになる悪寒がする orz

＊追記 ＠ Feb. 17
Wikipedia の 関連ページ を見てみたらほとんど同じような内容が書かれていてショック……。先に確認しておけば良かった。

The Python Challenge

前々から気になっていた The Python Challenge に取り組み始めました。これは、問題が出題されているページに隠されている次ページへ行くためのキーワードを、主に Python を使って探し出す一連のクイズ集です。言葉にすると分かりにくいかもしれないですが、例えば、1 問目のページに載ってる画像からある文字列を読み取るとか、そんな感じのものです。宝探しみたいなもの。今のところ全部で 34 問あるうちの 12 問目まで来ました。

ちょいちょい「これはプログラムの問題じゃないな」ってのも混ざっていますが、だいたいどの問題も良い勉強になりますね。正規表現の使い方とか、PIL の扱いとか。ヒントを求めるためにどうしてもページのソースを見るようになるので、自然と HTML の文法とかにも慣れてきます。問題を解いて初めて閲覧可能になる解答ページも非常に充実していて、別解が山のように載っているので（場合によっては他言語での解き方とかも議論されてます。Ruby とか bash とか）感心させられる事しきりです。自分がいかに Pythonic なコードを書けていないか痛感する……。可読性を保つ限りで簡潔に書く、ってのが楽しくも難しいところであります。

ipython, numpy, scipy, matplotlib, PIL (and sympy) の組み合わせが最強すぎて ipython シェルから出て来れません。os(.system), shutil モジュールの存在もその一因かな。

shutil.move のバージョン依存性

Python で標準的に用いることが出来るモジュールの中に shutil というのがあります。これは shell utility のことで、move や copy のようにファイルのコピーやディレクトリに関する処理を行なうシェルコマンドに対応するメソッドが含まれています（なお、mkdir とかは os モジュールに入っている）。

先ほど気づいたのですが、shutil.move の挙動は python のバージョンによって異なるのでメモしておきます。チェックした環境は CentOS 5.6, bash です。以下のスクリプト movetest.py を試してみましょう。

# -*- coding: utf-8 -*-
from os import mkdir
from shutil import move

mkdir("src")
mkdir("dst")
f = open("spam", "w")
f.write("test")
f.close()
move("spam", "dst") # まず、テキストファイル "spam" を dst/ に移す
move("src", "dst") # 次に src/ を dst/ に移す

その結果、バージョン 2.6/2.7 と 3.x では直感的な振る舞いをするのですが、2.4 はちょっと違っています：

• 2.4:
  Traceback (most recent call last):
    File "movetest.py", line 10, in ?
      move("src","dst")
    File "/usr/lib64/python2.4/shutil.py", line 190, in move
      copytree(src, dst, symlinks=True)
    File "/usr/lib64/python2.4/shutil.py", line 111, in copytree
      os.mkdir(dst)
  OSError: [Errno 17] File exists: 'dst' 
• 2.6 / 2.7, 3.x:
  dst/ 以下に src/ と spam が存在する。

2.4 では、移動先のフォルダに中身が入っていると move 出来ないようです。そこで、spam ファイルを dst に移さないで空の状態の dst に対して src を move しようとすると次のようになります。

• 2.4:
  src/ が dst/ にリネームされ、dst/ ディレクトリだけが残る。
• 2.6 / 2.7, 3.x:
  dst/ 以下に src/ が存在する。

2.4 では普通のファイルをリネームするかのようなはたらきをします。

これで分かったように、普通のシェルコマンドと同じ挙動を期待するなら少なくとも 2.6 以降を利用するべきだということでしょう。 2.5 については確認してないので分かりませんが……。

ちょっと時代遅れな記事だったかもしれないですね。

＊追記：
Python 2.5 は 2.4 と同じ振る舞いをするようです（MacOS 10.6 で確認）。2.5 --> 2.6 のタイミングで何かあったんでしたっけ。

コアの ID を指定して並列計算する方法（mvapich2 編）

並列計算をする際、通常は "mpirun -np 4" のように使用するコアの数のみを指定して実行しますが、何らかの理由があってコアの ID まで指定したい場合にはどうすればいいのでしょうか。mvapich2 を利用する場合の解決策を知ったのでここにメモしておきます。方法は以下の通り。

1. 計算を行う際のホスト名を例えば "host" というファイルに書いておきます。書き込むホスト名の先頭にはスペースを、行末には改行記号をつけておく必要があるようですが理由は知りません（試したところどちらかが欠けている場合にはエラーが出て実行出来ませんでした）。
2. mpirun_rsh -np 4 -hostfile ./host MV2_CPU_MAPPING=0:1:2:3 ./a.out と打ち込みます。MV2... が使用するコアを指定するオプションです。

あるジョブを特定のコアがすでに受け持ってくれているような場合で、それらの邪魔をせずに他のジョブを実行させたい時などは役に立つと思います。
■参考 URL（質問者が僕です）：http://goo.gl/stEjh

それにしても、高い負荷がかかっているコアを自動的に避けて複数のジョブを実行するってのはバッチ処理とかしないと実現出来ない事なんだろうか？

おひさ

あけましておめでとうございます、と言う機会はとうに逃してしまいました。お久しぶりです。とりあえず元気でやってます。ブログを更新するのもいつ以来か忘れてしまいましたが、雑多な近況なんかを書いておきます。

■バーゲンに行ったよ
もう２５にもなるというのに、相変わらずお年玉をもらっています。両親と、おばあちゃんとから。さすがに社会的にはどうかと思ってはいるんですけどねー、潤いが欲しいのは誰でも同じ。
そんなわけで、今月の初めころに脱非リア充を目論む友達に連れられてバーゲンに行ってきました。梅田の HEP とかそこら辺だったかな？春まで使えそうなジャケットと手提げのカバンを買ったら諭吉さんが数人出兵してしまって、潤いは一瞬で失われてしまいましたとさ。実用品が手に入ったのと、そもそも友達と買い物に長時間行ったこと自体が楽しかったので後悔は無し。■

■TOEFL はどうなったのか
そういえば昨年11月下旬に受けた TOEFL ITP の結果が来てるので公開してみます。

Listening Comprehension: 49
Structure & Written Expression: 57
Reading Comprehension: 63
Total Score: 563

Total Score は各項目の単純な和にはならないのですね。TOEFL ITP には二つの Level があるそうですが、今回のは Level 1 の方でした（Level 2 の方が簡単）。満点は660点で、ペーパー版と同等とのことです。この資料（PDFが表示されるのでご注意を）によれば ITP で563点というのは TOEIC でいうところの 750-800 点あたりに対応するので、まあ悪くはない出来だったと思います。昔から苦手なリスニングはやっぱり点が良くないですね。普段から論文を読んだりしてるので文章読解はそれなりになるし、日常的に英語を聞く必要性は無い環境にいるので当然といえば当然ですが、改善したいところではあります。

前の記事にも書いた通り、しばらくは TOEFL 対策みたいなことをせずにどれだけ点数を上げられるのか試していこうかなと考えてます。まあそんな悠長なこと言っているうちは点なんて上がらないか……。■

■NVIDIA Tesla 導入
僕が今やっている研究はタンパク質の構造変化に関与しています。分子動力学によりアプローチしているため、とにかく計算は速ければ速いほど良いという状況にあります。分子構造の可視化とか、二次構造含量の分析とか、何をするにしてもまず系のトラジェクトリデータを計算する必要があるからです。

１２月末までの段階では Xeon X5690x12 で計算していたのですが、ボスの協力もあって今月から NVIDIA の Tesla M2090 を2機導入しました。これは、いま使っている Amber という MD パッケージが CUDA を使った高速計算にも対応したのを利用するためです。ベンチマークページによれば4,5倍は速くなるとのことでしたが、実際に今研究している系で試してみると7倍くらい高速化されてました（詳細は省きますが、NTP 条件下で 40,000 原子系を扱うと 23ns/day くらい。以前は良くて 3.5ns/day だった）。Ubuntu 10.04, OpenMPI を使っていたのを CentOS 5.6, MVAPICH2 にしたり、intel MKL を使うように設定し直したのでそもそも基本的なセットアップが違いますが、それにしても目覚ましい進歩だということには変わりありません。Amber のインストールに関してはちょっと苦労した点もあったので、別記事にまとめておきます。

僕は未熟な Python 使いなのですが、CUDA API である PyCUDA というのが存在するのでそのうち使ってみたいです。まずは普通の並列計算をきちんと修めてからですけどね……。Python 大好きです。最近は Bioinformatics Programming Using Python をちょこちょこ読んで色々勉強してます。XML とかの扱いを知りたくて。■

■処女論文が近い？
そろそろ論文書こうか、という話がボスとの間で持ち上がってきています。もうD1も終わりだしね、時期的には遅すぎるくらいなんですが、ようやく結果がまとまって来たので何とかしようということです。これが初論文。
どうなんだろう、データの解釈が圧倒的に足りていない気がする。先行研究との関連付けなんかも全く不十分だし、結果が割りと面白いだけに justification みたいなことをしっかりやらないと絶対 reject されると思う。これからが勝負です。とはいえ今月末が第一稿の締めきりということになっているorz■