Python で並列計算 (multiprocessing モジュール) | 複数の引数を取る関数を map() メソッドで並列に走らせる
いまやノートブックでもデュアルコアやクアッドコアが当たり前になってきたので、似たような処理を延々と繰り返すようなデータ解析のプログラムなどは並列化するとかなりの恩恵が得られる。Python ではバージョン2.6から multiprocessing モジュールというのが標準ライブラリに入っており、比較的簡単に並列計算のスクリプトが書けてしまう。
個人的に並列化できると一番うれしいのは、統計量の有意性の検定のためにサロゲートデータを大量に作って、そこから(帰無仮説に基づいた)統計量の分布を推定するときの、それぞれのサロゲートごとの統計量の計算だ。この場合、個々の計算は完全に独立なので、プロセス間の通信などは考える必要がなく、非常に単純に並列化できる。ということで早速スクリプトを書いてみたが、まあ書くだけならすぐ書けるのだが、使い勝手や汎用性を求めると、ちょっとややこしいことを考えなければならなくなってしまった。
とりあえずの要件として、求めたい統計量の計算をするための関数はあらかじめ手元にあるので、これをそのまま使って並列化したいということがある。つまり、並列化のためだけにわざわざ新しい関数を書くというようなめんどくさいことはしたくない。multiprocessing モジュールで並列化する場合、一番単純なのは Process クラスのオブジェクトに並列で計算したい関数を設定して(multiprocessing.Process(func, args))、それを start() メソッドで走らせるというやり方だが、これだと関数の返り値を受け取るのに、親プロセスと子プロセスの間の通信を使わないといけない。そのためにはこれらのプロセスの間にパイプを通して(multiprocessing.Pipe())、子プロセスの関数にこのパイプ経由で計算結果を出力するようにさせないといけないのだが、これには当然関数の書き換えが必要となってくる。
何かもっと楽な方法はないかと調べてみると、Process のかわりに Pool というのを使うと、返り値を楽に取り扱えることが分かった。すなわち、results = multiprocessing.Pool().map(func, arglist) とすると、arglist 内の要素のそれぞれを引数とした func を並列に計算し、その結果をリストとして results に格納してくれる。同時に走るプロセス数は、マシンのコア数を調べて勝手にそれに合わせてくれる(Pool()の引数で指定することも可)。これはとても楽だ。
が、実際に自分の関数を使って並列化しようとしてみると、困ったことに気付いた。統計量の計算をするとき、しばしばいくつかのパラメータを設定する必要があるのだが、手元の関数では、これを複数の引数を渡すことで実装している。ところが、Pool オブジェクトの map() メソッドでは、関数に複数の引数を渡すことができない。リストにして渡せばひとつの引数にまとめられるが、そうすると再び関数の書き換えが必要になってしまう。これは困った・・・
ということで、なんとかして既存の関数を書きかえることなく、複数の引数を渡せるようにしようといろいろ試してみた。いまのところ、以下のようなラッパー関数を使うのが一番簡単で、かつ使いやすいんではないかというところに落ち着いた。
def argwrapper(args): ''' ラッパー関数 ''' return args[0](*args[1:]) def myfunc(a, b): ''' 並列に計算したい関数 ''' return a*b if __name__ == '__main__': from multiprocessing import Pool p = Pool() func_args = [] for a in xrange(1,10): for b in xrange(1,10): func_args.append( (myfunc, a, b) ) results = p.map(argwrapper, func_args) print results
argwrapper() はリストを引数に取り、このリストの最初の要素を関数として実行し、それ以降の要素をその関数に渡す引数とする。main 部分では、計算したい関数とそれへの引数をまとめたリストを作り、それを map() メソッドで argwrapper() に渡している。この方法なら、あらゆる既存の関数に関して、複数の引数を与え map() メソッドで並列に計算させることができるはず。
複数の引数を与えるだけなら map() メソッドのかわりに apply_async() メソッドを使ってもいいのだが、その場合、結果をひとつづつ get() メソッドで取り出さないといけない。結果が直接リストとして返ってくる上記の方法のほうが、個人的には使い勝手がいいと思う。
Windows XP SP3 で「ネットワーク接続」のアイコンが消えてしまった時の対処法
突然、無線・有線LAN共に機能しなくなり、接続設定を確認しようとして「ネットワーク接続」を開いても、「ローカル エリア接続」や「ワイヤレス ネットワーク接続」のアイコンが消えてしまっておりどうしようもない(にもかかわらずデバイスマネージャでは異常なし)という状況。
これまで2度同じ症状に遭遇したのだが、1度目はOSの再インストールで復旧。2度目は以下の方法で再インストールなしで復旧した。
再起動後、デバイスドライバは自動的にインストールされ(インストールCD-ROMなどのメディアは不要)、「ネットワーク接続」のアイコンが復活する。
ただし、セキュリティ・キー等は再度設定し直す必要がある。
また、接続の名称が、「ワイヤレス ネットワーク接続 2」のように、末尾に番号が振られたものに変わる。
↓の Microsoft Support の記述を参考にした。
How to troubleshoot missing network connections icons in Windows Server 2003 and in Windows XP
Windows版Python2.6 で matplotlib.pyplot をインポートしようとすると Memory Error が出る問題の解決法
誰にでも出る問題というわけではなく、標準以外のフォントをインストールしたシステムでのみ出る(可能性のある)問題らしい。
解決法は、C:\Python26\Lib\site-packages\matplotlib\font_manager.py というファイルの206行
if not local:
の前に以下の行を挿入する
local = None
インデントは、206行の内容" if not local:"と同じレベルになるように。
↓のサイトを参考にさせていただいた。
pylabでMemoryErrorが消えなかった原因がようやくわかった!
Python (Scipy) による神経活動データ解析入門・その2 | プロットあれこれ
前回定義した配列 ecog_mV を、いろいろプロットを変えて詳しく調べていく。
プロット範囲の指定
前回の例では2,000秒間ほどのデータがすべてプロットされていたが、より細かい時間スケールの変動を見るため、プロット範囲を短く設定する。
IN [*]: xlim( 0, 1 )
これでデータの最初の1秒間がプロットされる。同様に、
IN [*]: xlim( 1000, 1001 )
とすると、1,000秒から1,001秒までの1秒間がプロットされる。
(y軸の範囲の指定は ylim( y_min, y_max ) で、x軸y軸の同時指定は axis( [x_min, x_max, y_min, y_max] ) で行う。)
0−1秒の時間帯は麻酔剤注入前の覚醒状態(awake period)、1,000−1,001秒は麻酔が十分効いている状態(sleep period)と考えられるが、これらの間の ECoG 信号を比較してみると、前回確認した、変動の振幅の大きさだけでなく、変動の時間スケールも全く異なっていることが分かる。すなわち、0−1秒では1秒間に20回ほどの規則的な振動が見られるが、1,000−1,001秒ではそれが見られず、大振幅のゆっくりとした振動に変化している。
対数プロット
信号の変動の時間スケールを定量的に知りたい場合、特に、特定の周波数の振動の有無を知りたい場合、まず調べるのは信号のパワースペクトル密度である。SciPy(+Matplotlib) では、コマンド一つでスペクトル密度の計算とプロットがいっぺんにできてしまう。
今回は覚醒状態・麻酔状態で別々にスペクトル密度を見たいので、元のデータのうち最初の500秒間を覚醒状態のサンプル、1,000秒から1,500秒までを麻酔状態のサンプルとして別の配列にコピーし、それぞれについてスペクトル密度を計算する。
IN [*]: ecog_awake = ecog_mV[0:500*1000] IN [*]: ecog_sleep = ecog_mV[1000*1000:1500*1000] IN [*]: psd( ecog_awake, Fs=1000, label='awake' ) IN [*]: psd( ecog_sleep, Fs=1000, label='sleep' )
psd に与えているデータ以外の引数は、Fs がサンプリング周波数、label は後で凡例を作るときに使うラベル。
この方法では、縦軸の単位は自動的に dB/Hz でプロットされる。デシベル(dB)は信号強度を対数変換したものなので、つまりこれは実質的に片対数プロットになっている。横軸も対数スケールにするには、
IN [*]: xscale( 'log' )
とする。その結果が以下のグラフ。(ただし上記の設定以外に、grid()でグリッドを、legend()で凡例を表示させている。)
覚醒状態のグラフ(awake)は、右下がりの直線上に 15-20 Hz の周波数帯でピークが乗った形になっている。両対数グラフで直線になるということは、元の信号がいわゆる1/fゆらぎを持っていることを示している。15-20 Hz のピークは、最初に示したグラフで見えた細かい振動に対応しており、ヒトEEGで閉眼安静時に顕著に見られるアルファ波に相当するものと思われる。麻酔状態(sleep)ではこのピークは消失し、代わりに 40-80 Hzあたりに肩のような盛り上がりが見える。
Python (Scipy) による神経活動データ解析入門・その1 | .mat ファイルの読み込みからデータのプロットまで
前提として、Python, iPython, NumPy, SciPy, Matplotlib はインストール済み、Neurotycho のデータはダウンロード済みであるとする。
SciPyの起動
まずはシェルプロンプトから
% ipython -pylab
で iPython(Python の対話的シェル)を起動。"-pylab" を付けることで、MATLAB 風のデータプロットのコマンドをサポートするパッケージ(Matplotlib) を組み込んだモードになる。
立ち上がった iPython のプロンプトから
In [*]: from numpy import * In [*]: from scipy import *
で NumPy, SciPy のモジュールをすべてインポート。これでひと通りの機能が使えるようになる。
データの読み込み
MATLAB形式のデータ(.mat)を読み込むには scipy.io.loadmat() を使う。まずこれをインポートする。
In [*]: from scipy.io import loadmat
Neurotycho > download > sleep task では ECoG データを4つのファイルに分割して公開しているが、まず最初のデータファイル(ECoG_100604-1_32.mat) を読み込んでみる。
In [*]: dataset = loadmat( 'ECoG_1006040-1_32.mat' )
とし、データファイルの内容を辞書オブジェクト(=連想配列) dataset にコピーする。
データの確認
まずおもむろに
In [*]: dataset
とすると、辞書オブジェクト dataset の概要が表示され、このデータが以下の4つの要素を持つことが分かる:
- ラベル 'X': int16型の配列オブジェクト
- ラベル '__globals__': 空のリスト
- ラベル '__header__': ファイル情報を記述した文字列
- ラベル '__version__': バージョン数を記述した文字列
これらのうち、'X' とラベルされた最初の要素がECoGのデータ本体であると思われる。この要素のみを取り出してみる。
In [*]: ecog_data = dataset['X']
以下、この配列 dataset の情報を得るためのメソッドをいくつか紹介していく。まず、配列の次元を知るには dim メソッドを使う。
In [*]: ecog_data.dim Out[*]: 2
となり、2次元配列であることが分かる。各次元の大きさを知るには shape メソッドを使う。
In [*]: ecog_data.shape Out[*]: (32, 2062570)
となり、データファイルの名前や、データに付属している ReadMe.txt の内容から推測すると、第1次元が ECoG のチャンネル、第2次元がデータ点に対応しているらしいと分かる。すなわち、ecog_data[n] を参照すればn番目のチャンネルのデータが取り出せることになる(ただし 0 <= n <= 31 に注意)。
データのプロット
データを入手してまず最初にすべきことは、それを目で見て確認すること。そのためにデータをプロットしてグラフにする。最も単純な方法は
In [*]: plot( ecog_data[0] )
とすれば、配列 ecog_data[0] の要素をそのままプロットできる。だがこれだと、横軸が配列の要素番号となってしまい、時間との対応がつかない。また縦軸の単位もわからない。
Neurotycho > download の記述に、ECoG データの単位はマイクロボルト(uV)、サンプリング周波数は 1kHz だとあるので、まず、縦軸をミリボルト単位でプロットするため、データの単位を変更して新しい配列 ecog_mV に移す。
In [*]: ecog_mV = ecog_data[0] / 1000.0
横軸に関しては、秒単位でプロットするため、以下のようにして横軸のグリッドを表す配列 n_sec を生成する。
In [*]: n_sec = arange( 0, ecog_data[0].shape[0] ) In [*]: n_sec = n_sec / 1000.0
やっていることは、まず 0 から ecog_data[0] の要素数まで、1づつ増える整数の配列を生成し、それをサンプリング周波数で割っている。この配列 n_sec と先程の単位変更したデータ ecog_mV を使うと、
In [*]: plot( n_sec, ecog_mV )
として、横軸が時間(sec)、縦軸が ECoG 信号(mV)のグラフが描ける。軸にラベルを付けるには
In [*]: xlabel( 'Time (sec)' ) In [*]: ylabel( 'ECoG signal (mV)' )
でそれぞれ横軸と縦軸にラベルが付けられる。結果が以下のグラフ。
ReadMe.txt に、測定開始後631秒の時点で麻酔剤を注入したとあるが、その100秒後あたりから ECoG 信号の振幅が増大していることが分かる。また、それ以前・以後の時間帯では信号の強度は概ね定常であることも分かる。
Multitaper method の直観的な説明・その4
前回の続き。
ある行列に関する固有値問題の解(=固有ベクトル)がスペクトルの漏れに関して最適な窓関数となることを見たが、一般的に固有値方程式の解は行列のランクの数だけあるので、窓関数もその数だけ得られることになる。
これら複数の窓関数がどのような形状をしているかは、例えば wikipedia の Multitaper の項で見ることができる。ここにある図では、大きいほうから3つの固有値に対応する固有ベクトルとして得られる窓関数が描かれている。最大固有値に対応する固有ベクトルはハミング窓のような形状、第2固有値では正弦関数のような形状になるが、その先、より小さい固有値に対応する固有ベクトルはゼロ交差の数がひとつづつ増えていき、より振動性の強い関数形となる。
さて、ここでは証明は省略するが、これらの固有値は1より小さい正の値をとることが知られており、またその値は、それぞれの固有値に対応する窓関数が持つパワー全体のうち、あらかじめ設定した(前回の説明を参照)許容周波数範囲内に納まっているパワーの割合を表す。つまり、1に近い固有値に対応する固有ベクトルは(スペクトルの漏れが少ないという意味で)「良い」窓関数であり、0に近い固有値に対応する固有ベクトルは「悪い」窓関数だ、ということだ。
そして、また証明は省略するが、固有値の分布は前述の許容周波数範囲([tex:-W