良いもの。悪いもの。

世の中のことをもっと知るにはどうしたら良いだろうと思うときがある。世の中の多くの事柄はログやデータに落とされる。Googleなどの検索サイトは良い例だろう。さて、そのログやデータをどうすれば良いのか? 多くの場合、視覚化が有効な手段となる。 まずは身の回りの日常的なデータやログを何とかしたい。ただ、日常のデータを視覚化するのに数十行以上のコードは書きたくない。まるで息をするかのごとく自然に視覚化を行いたいのだ。そのためには1～2行、長くて数行で済ませることが必要だ。そこでPythonとmatplotlibを使う。加えて、IPythonがあればなお良い。IPythonの導入については以前のブログ記事である IPythonの埋め込みプロットが素晴らしい を参考にして欲しい。 まずは事前にnumpyとmatplotlibをインポートしておく。できればscipyも。 >>> from numpy import * >>> from pylab import * 短いコードで視覚化を行うためには、Pythonの内包表記は必須だ。例えば、5, 2, 1, 5, 8をデータとするグラフを書きたいのならIPythonを使って以下の1行で実現できる。 >>> plot([5, 2, 1, 5, 8]) 数値が行ごとにinput.txtというファイルに書かれていた場合は以下の通り。 >>> plot([int(x) for x in file("input.txt")]) map関数を使えばもっとスマートに書ける。ファイル内の数値を文字列として配列で読み込み、それらの文字列をmap関数によりintで整数に変換する。 >>> plot(map(int, file("input.txt"))) さて、ログが2次元で書かれていた場合はどうするのか。例えば以下のデータがinput.txtに書かれていたとする。 2, 2.5 5, 6.2 6, 3.6 7, 6.3 10, 1.9 この場合、以下のようにすればプロットできる。ファイルから行ごとの文字列を読み出して、それを","を区切りとしてリス...

先日、 Tokyo.SciPy #3 に参加して、SciPyの生みの親である Travis Olipant氏のセッション を聞いた。その時に最近の IPython では埋め込みプロットができることを知ったので早速入れてみたところ、これがとてもクールでカッコよかったのでここで紹介したい。 埋め込みプロットはIPythonのバージョン0.11からできるようになっている。ただ、自分が使っているLinux環境では普通に持ってくるとバージョン0.10が入ってしまうので使えない。そこで、IPythonの公式サイトから最新版のバージョン0.12を持ってきて入れてみた。 ただ、埋め込みプロットができるのはターミナル版ではなくQt版のシェルなので、それに関連するライブラリなども一緒に入れなくてはならない。IPythonの起動時に何々のライブラリが足りないとかいろいろと文句を言われるが、足りないものをyumなりapt-getなりソースコードをコンパイルなりして順次入れていけば使えるようになると思う。 あと、IPythonの設定ファイルの構成が0.10から大きく変わった。以前は .ipython/ 内の pythonrc を利用していたのだが、それがなくなって、代わりに .ipython/profile_default/ 内の ipython_config.py と ipython_qtconsole_config.py が使われるようになった。デフォルトの設定ファイルは以下のコマンドで作ることができる。 $ ipython profile create そして、埋め込みプロットを使うためのQt版IPythonの起動コマンドは以下の通り。 $ ipython qtconsole --pylab=inline これでQt版シェルが立ち上がって埋め込みプロットを表示できるようになる。 せっかくなので実際に埋め込みプロットを使ってみる。ここでは自分のお気に入りのマンデルブロ集合をプロットしてみた。因みに、シェル上での複数行のソースコードの編集も改良されたので入力がかなり楽になった。冒頭にスクリーンショットを貼っておく。 rangex = arange(-2.0, 1.0, 0.01) rangey = arange(-1.0, 1.0, 0.01) (X, ...

Python + RPy + SciPy で階層的クラスタリングを行う。下記のtest.dat(30個の5次元データ)に対して閾値0.5として計算する。 clustering.py test.dat 0.5 標準出力: Data file: test.dat Element file: element.out Output file: data.out Cluster file: cluster.out EPS file: cluster.eps Threshold: 0.500000 (distance) 30 x 30 = 900 Cluster 1 : 5 Cluster 2 : 8 Cluster 3 : 1 Cluster 4 : 1 Cluster 5 : 12 Cluster 6 : 2 Cluster 7 : 1 Total : 30 Number of clusters: 7 test.dat : 30個の5次元データ。 -2.39257 0.39147 0.21229 -0.84501 -0.02255 -2.13082 0.44121 0.12561 -0.91345 -0.12832 -2.23498 0.26590 -0.35356 -0.03484 -0.17437 -2.27747 0.29900 -0.05188 -0.63788 -0.13000 -2.11555 0.39492 0.11465 -0.86242 -0.34207 -2.31947 0.26036 0.47846 -0.71782 0.07270 -2.32941 0.38637 0.44304 -0.75346 -0.04302 -2.30690 0.59296 0.22947 -0.82801 -0.09882 -2.29088 0.53157 0.03204 -0.75087 -0.05802 -2.36682 0.34835 0.23964 -0.72942 -0.09720 -2.11368 0.52013 0.27528 -0.81122 -0.20050 -2.27156 0.21167 0.57165 -0.66193 -0.03302 -2....

分子動力学(MD)計算で得られた座標データに対して、SciPyで主成分分析をしてみた。通常、座標のトラジェクトリーに対して最小二乗法によるフィッティングを行ってからPCAを行うのだが、ここで示すPCAの仮引数trjに渡すトラジェクトリーはフィッティングを前提にしてある。 実際のコードではフィッティングやら、ARPACKの利用やら、原子の重み付けの方法やらをオプションで引き渡して処理するようになっているのだが、ここでは内容を明確にするために、MDの一般的なPCAに絞ったコードを示す。 MD計算により10個の構造が出力される場合、N原子を持つ一構造のXYZ座標のリスト([N][3])をトラジェクトリー([10][N][3])として渡す(trj)。原子量は要素がN個のリストとなる(mass)。原子量のリストを渡さなければ原子の重み付けは行われず、Essential Dynamics (ED)の計算となる。 以下、コード。 from scipy import * def pca(trj, mass=[], num=2): """trj: 座標データ. / mass: 原子量のリスト(重み付けが必要であれば) / num: 求める主成分の数""" # 座標の平滑化. R = [mat(crd).A1 for crd in trj] NDIM, N = len(R), len(R[0]) # 平均構造の作成. R_mean = zeros(N) for i in range(NDIM): R_mean += R[i] R_mean /= NDIM # 各構造と平均構造の差. for i in range(NDIM): R[i] -= R_mean # 必要なら原子の重み付けを行う. if mass: sq_mass = [sqrt(m) for m in mass] for i in range(NDIM): for j in range(N): R[i][j] *= sq_mass[j // 3] # ...

SciPy で固有値・固有ベクトルを求める場合、scipy.linalg.eigを利用する。SciPy自体、LAPACKを使っており、大抵の場合は問題がない。しかし、巨大な行列での固有値問題を解く場合、メモリなどの関係で計算が破綻する場合がある。FORTRANであれば、 ARPACK が利用されることが多い。 SciPyでは現時点でARPACKの利用はできない。しかし、開発バージョンにおいては、既にARPACKのラッパーが提供されており、疎行列演算のモジュール群であるscipy.splinalgの scipy.splinalg.eigen を利用することになる。 個人的に、分子動力学(MD)計算によるタンパク質座標データの軌跡に対して主成分分析(PCA)を行いたかったのだが、2GB程度のメモリだと、scipy.linalg.eigではメモリオーバーで破綻してしまった。そこで、scipy.splinalg.eigenを利用したところ、問題なく処理することができた。 scipy.splinalg.eigenは、巨大配列の固有値を解く場合には有効だが、開発バージョンであるので問題が出るかもしれない。その辺は自己責任で。

PythonのモジュールであるSciPyを使った行列演算を示す。scipy.linalgではBLASやLAPACKを利用しているので、高速な演算が可能になっている。 >>> from scipy import * >>> a = [[1, 3, 5], [2, 5, 1], [2, 3, 8]] # リスト >>> A = mat(a) # リスト→行列 >>> A.tolist() # 行列→リスト [[1, 3, 5], [2, 5, 1], [2, 3, 8]] >>> A # 行列 matrix([[1, 3, 5], [2, 5, 1], [2, 3, 8]]) >>> A.I # 逆行列 matrix([[-1.48, 0.36, 0.88], [ 0.56, 0.08, -0.36], [ 0.16, -0.12, 0.04]]) >>> A.H # エルミート行列 matrix([[1, 2, 2], [3, 5, 3], [5, 1, 8]]) >>> A.T # 転置行列 matrix([[1, 2, 2], [3, 5, 3], [5, 1, 8]]) >>> A.A # 行列→配列 array([[1, 3, 5], [2, 5, 1], [2, 3, 8]]) >>> A.A1 # 行列→1次元配列 array([1, 3, 5, 2, 5, 1, 2, 3, 8]) >>> A[1] # 2番目の行を取り出す matrix([[2, 5, 1]]) >>> A.T[1] # 2番目の列を取り出す matrix([[3, 5, 3...