高速化の小技

frame fileのまとめ読み

例として、full frame 1ファイルから3個の3軸磁束計の9チャンネルを読んでみる。

   1 from gwpy.timeseries import TimeSeries, TimeSeriesDict
   2 source = '/data/full/12677/K-K1_C-1267797600-32.gwf'
   3 channels = ['K1:PEM-MAG_BS_BOOTH_BS_X_OUT_DQ',
   4             'K1:PEM-MAG_BS_BOOTH_BS_Y_OUT_DQ',
   5             'K1:PEM-MAG_BS_BOOTH_BS_Z_OUT_DQ',
   6             'K1:PEM-MAG_EXC_BOOTH_EXC_X_OUT_DQ',
   7             'K1:PEM-MAG_EXC_BOOTH_EXC_Y_OUT_DQ',
   8             'K1:PEM-MAG_EXC_BOOTH_EXC_Z_OUT_DQ',
   9             'K1:PEM-MAG_EYC_BOOTH_EYC_X_OUT_DQ',
  10             'K1:PEM-MAG_EYC_BOOTH_EYC_Y_OUT_DQ',
  11             'K1:PEM-MAG_EYC_BOOTH_EYC_Z_OUT_DQ']
  12 for i in range(9):
  13    data = TimeSeries.read(source, channels[i], format='gwf.lalframe')    

21.863秒かかった。次に、TimeSeriesDictでまとめて読みだすと、

   1 from gwpy.timeseries import TimeSeries, TimeSeriesDict
   2 sources = ['/data/full/12677/K-K1_C-1267797600-32.gwf']
   3 channels = ['K1:PEM-MAG_BS_BOOTH_BS_X_OUT_DQ',
   4             'K1:PEM-MAG_BS_BOOTH_BS_Y_OUT_DQ',
   5             'K1:PEM-MAG_BS_BOOTH_BS_Z_OUT_DQ',
   6             'K1:PEM-MAG_EXC_BOOTH_EXC_X_OUT_DQ',
   7             'K1:PEM-MAG_EXC_BOOTH_EXC_Y_OUT_DQ',
   8             'K1:PEM-MAG_EXC_BOOTH_EXC_Z_OUT_DQ',
   9             'K1:PEM-MAG_EYC_BOOTH_EYC_X_OUT_DQ',
  10             'K1:PEM-MAG_EYC_BOOTH_EYC_Y_OUT_DQ',
  11             'K1:PEM-MAG_EYC_BOOTH_EYC_Z_OUT_DQ']
  12 data = TimeSeriesDict.read(sources, channels, format='gwf.lalframe')

5.166秒と短縮された。

次に、読むファイル数を10個にしてみる。31.008秒かかった。TimeSeriesDict.read の引数に nproc=コア数 を入れると並列化され、より速くなる。

   1 data = TimeSeriesDict.read(sources, channels, format='gwf.lalframe', nproc=10)

結果は5.485秒。ちなみに、ファイル数が1のときにはnprocを増やしても変わらなかった。


解析コードをDevelopする際 (特にevent baseの解析)

解析コードを構築する際、少し書き換えるたびに何度も何度も読みだすことになるが、これが結構時間を食う。

特にグリッチ解析や地震の解析などevent baseの解析の場合は、ごく限られた時間領域のデータを抜き出すために無駄に多くのデータを読む。

といったプロセスがloopで繰り返されることになる。

このような作業の際は、以下のようにeventごとに必要なチャンネル・時間領域を抜き取ったTimeSeriesDictをlistに詰め、これをpickleもしくはshelveに保存しておくとよい。

   1 def SingleEvent(TriggerTime): #eventの時刻を引数にする
   2    data  = TimeSeriesDict.read(...)
   3    state = StateVector.read(...)
   4    return data, state
   5 
   6 def main():
   7    ...
   8    data  = []
   9    state = [] 
  10 
  11    for t in tqdm(TrigTimes): #予め各eventの時刻をリストにしておく
  12         d, s = SingleEvent(t)
  13         data.append(d)
  14         state.append(s)
  15 
  16    f = shelve.open('./EventList.slv')
  17    f['data']  = data
  18    f['state'] = state
  19    f['TrigTimes'] = TrigTimes
  20    f.close()

このshelveを読むには、以下のようなコードを書けばよい。

   1 import shelve
   2 f = shelve.open('./EventList.slv')
   3 for data, state, t in zip(f['data'], f['state'], f['TrigTime']):
   4    ...#後は好きなように解析

もしくは、「10 eventの時系列見たいなー」といった場合は、インタープリンタで

>>> import shelve 
>>> f = shelve.open('./EventList.slv') 
>>> f['data'][10].plot().show()

と打つだけですぐに表示される。 爆速になるので超絶おススメ。


Cython

Cythonについて調べると、setup.py なるファイルを作れという記事が多くヒットするが、現代ではこれは必要なくなった。(ここを参照)

Cythonのコンパイル -> 実行を一発でできるようにする

Cythonを使って高速化したはいいけど、毎回コンパイルするのは面倒だ。複数のコードになってしまうのは嫌だ。という場合の小技。

home directoryに .cython.sh という名前で以下を書き込んだシェルスクリプトを作る。

   1 #!/usr/bin/bash
   2 set -e
   3 trap 'echo compile error' ERR
   4 code=$1
   5 cythonize -i $code > .cython_log.txt
   6 rm ${code%.pyx}.c
   7 echo 'import '${code%.pyx}  > .a.py
   8 echo ${code%.pyx}'.main()' >> .a.py
   9 python .a.py "${@:2}"

また、 ~/.bashrc に次の一行を書き足す。

   1 alias cythonsh='sh ~/.cython.sh'

これで .bashrc を有効にすれば、ターミナルで

> cythonsh (cythonコード).pyx 引数

と打つだけでコンパイルから実行まで行われる。

2, 3行目は、コンパイルに失敗した時にスクリプトが止まるようにするおまじない。

6行目は .c ファイルを残さないようにしている。

9行目の末尾の "${@:2}" は、引数をpythonに渡すおまじない。