良いもの。悪いもの。

しきい値は閾値(いきち)の慣用読み。以前はしきい値と云っていたが、閾はしきいとは読めないので、それからは閾値(いきち)と云うようになった。調べたところ、どちらも普通に使われるようだ。もっとも、分野によってどちらが良く使われるかは異なるらしい。 慣用読みと云えば、「輸出」は正しくは「しゅしゅつ」だし、「消耗」も「しょうこう」が正しい読み方だ。そう言えば、医学分野では「口腔(こうこう)」を「こうくう」と読む。これは日本解剖学会により定められているらしい。

備忘録:PythonでHTTPのDigest認証を行う をほとんどそのまま。忘れないようにメモ。 Python Challenge のある問題で大変役立った。 import urllib2 top_url = "www.foo.org" url = "www.foo.org/bar.html" username = "user" password = "pass" passman = urllib2.HTTPPasswordMgrWithDefaultRealm() passman.add_password(None, top_url, username, password) authhandler = urllib2.HTTPBasicAuthHandler(passman) opener = urllib2.build_opener(authhandler) urllib2.install_opener(opener) pagehandle = urllib2.urlopen("http://" + url) print pagehandle.info()

分子動力学(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は、巨大配列の固有値を解く場合には有効だが、開発バージョンであるので問題が出るかもしれない。その辺は自己責任で。