Differences between revisions 11 and 12
Revision 11 as of 2020-03-23 15:45:14
Size: 5433
Comment:
Revision 12 as of 2020-03-24 13:22:40
Size: 6752
Comment:
Deletions are marked like this. Additions are marked like this.
Line 47: Line 47:
解析コードを構築する際、少し書き換えるたびに何度も何度も読みだすことになるが、これが結構時間を食う。

特にグリッチ解析や地震の解析などevent baseの解析の場合は、ごく限られた時間領域のデータを抜き出すために無駄に多くのデータを読む。
解析コードを構築する際、少し書き換えるたびに何度も何度も読みだすことになるが、これが結構時間を食う。<<BR>>
特にグリッチ解析や地震の解析などevent baseの解析の場合は、ごく限られた時間領域のデータを抜き出すために無駄に多くのデータを読む。<<BR>>
Line 94: Line 93:
と打つだけですぐに表示される。 と打つだけですぐに表示される。<<BR>>
Line 99: Line 98:

== Matplotlibを使わずにpyROOTで描画 ==

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

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

{{{#!highlight python
plot = data.plot().show()
}}}
11.86秒

{{{#!highlight python
from ROOT import TGraph
from array import array
n = data.size
x = array('d', data.times.value.tolist())
y = array('d', data.value.tolist())
g = TGraph(n, x, y)
g.Draw('apl')
}}}
4.96秒と、半分以下の時間になった。しかしコードを書く量は倍どころじゃなく増えたし、横軸も読めない。

{{{#!highlight python
from ROOT import TGraph
from array import array
from gwpy.time import from_gps

def gwTGraph(data):
    n = data.size
    x = array('d', (data.times-data.t0).value.tolist())
    y = array('d', data.value.tolist())
    g = TGraph(n, x, y)
    g.SetTitle(str(data.channel))
    g.GetXaxis().SetTitle('Time [seconds] from {0:%Y-%m-%d %H:%M:%S UTC}'.format(from_gps(data.t0.value)) )
    return g
}}}


-----


Line 106: Line 149:
Cythonを使って高速化したはいいけど、毎回コンパイルするのは面倒だ。複数のコードになってしまうのは嫌だ。という場合の小技。
Cythonを使って高速化したはいいけど、毎回コンパイルするのは面倒だ。複数のコードになってしまうのは嫌だ。という場合の小技。<<BR>>
Line 128: Line 170:
と打つだけでコンパイルから実行まで行われる。

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

6行目は .c ファイルを残さないようにしている。
と打つだけでコンパイルから実行まで行われる。<<BR>>
2, 3行目は、コンパイルに失敗した時にスクリプトが止まるようにするおまじない。<<BR>>
6行目は .c ファイルを残さないようにしている。<<BR>>

高速化の小技

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 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()

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


Matplotlibを使わずにpyROOTで描画

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

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

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

11.86秒

   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


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

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