Differences between revisions 1 and 19 (spanning 18 versions)

高速化の小技

Contents

高速化の小技

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を増やしても変わらなかった。

for文をリスト内包表記に置き換える

pythonはループ処理が非常に遅いことで有名だが、積極的にリスト内包表記で置き換えることでかなり高速化される。
参考：リスト内包表記とif文を組み合わせるとき

慣れないうちは、解析をdevelopする際は(ループ数を減らして)for文で書いて、処理内容が決まった後で内包表記に翻訳するのがおすすめ。

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

解析コードを構築する際、少し書き換えるたびに何度も何度も読みだすことになるが、これが結構時間を食う。
特にグリッチ解析や地震の解析などevent baseの解析の場合は、ごく限られた時間領域のデータを抜き出すために無駄に多くのデータを読む。
といったプロセスがloopで繰り返されることになる。

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

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

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

こうして作ったshelveファイルは非常に軽くなるため、自分のローカルPCにコピーして解析するとこも可能だが、絶対にKAGRAコラボレーター外にはデータを流出しないように注意すること。破ると神罰が下ります。
特にDropBoxやGoogleDriveなどのクラウドストレージを利用している人は、これに同期されたディレクトリに置くと企業にデータを流したことになってしまう。

時系列データをDown samplingする / スペクトログラムの時点でBand limit

見たい周波数帯に比べてサンプリング周波数が無駄に大きいときには、あらかじめTimeSeriesに .resample() をかけて計算コストを下げると良い。

スペクトルのパーセンタイルを計算するときには、まず.spectrogram2()を使うが、その時点でcrop_frequencies() でband limitしておくとその後の処理が軽くなる。

上記のshelveに保存する前にもこれらをやっておくべし。

Matplotlibを使わずにpyROOTで描画

Matplotlibの描画は遅いことで有名であり、これを使っているgwpy.plotも同様に遅い。
代わりにpyROOTのTGraphを使って描画し、速度を比較してみた。

とあるデータ(16384 Hz sampling × 320.0 sec = 5242880点)を読み込んでから、時系列のplotが表示されるまでの時間を測った。

まずgwpyの場合、

   1 plot = data.plot().show()

11.86秒。次にpyROOTの場合、

   1 from ROOT import TGraph
   2 from array import array
   3 n = data.size
   4 x = array('d', data.times.value.tolist())
   5 y = array('d', data.value.tolist())
   6 g = TGraph(n, x, y)
   7 g.Draw('apl')

結果は4.96秒と、半分以下の時間になった。しかしコードを書く量は倍どころじゃなく増えたし、横軸も読めない。
そこで次のような関数を定義しておけば、だいぶ楽&まともになる。

   1 from ROOT import TGraph
   2 from array import array
   3 from gwpy.time import from_gps
   4 
   5 def gwTGraph(data):
   6     n = data.size
   7     x = array('d', (data.times-data.t0).value.tolist())
   8     y = array('d', data.value.tolist())
   9     g = TGraph(n, x, y)
  10     g.SetTitle(str(data.channel))
  11     g.GetXaxis().SetTitle('Time [seconds] from {0:%Y-%m-%d %H:%M:%S UTC}'.format(from_gps(data.t0.value)) )
  12     return g

plot.gw2.png

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に渡すおまじない。

リンク集

C言語は使えるけど最近Pythonを始めた」人向けの記事

KAGRA/Subgroups/PEM/PythonMemoJP/speedup (last edited 2021-08-04 18:33:51 by tatsuki.washimi)