= 高速化の小技 =

<<TableOfContents(2)>>

== frame fileのまとめ読み ==
例として、full frame 1ファイルから3個の3軸磁束計の9チャンネルを読んでみる。
{{{#!highlight python
from gwpy.timeseries import TimeSeries, TimeSeriesDict
source = '/data/full/12677/K-K1_C-1267797600-32.gwf'
channels = ['K1:PEM-MAG_BS_BOOTH_BS_X_OUT_DQ',
            'K1:PEM-MAG_BS_BOOTH_BS_Y_OUT_DQ',
            'K1:PEM-MAG_BS_BOOTH_BS_Z_OUT_DQ',
            'K1:PEM-MAG_EXC_BOOTH_EXC_X_OUT_DQ',
            'K1:PEM-MAG_EXC_BOOTH_EXC_Y_OUT_DQ',
            'K1:PEM-MAG_EXC_BOOTH_EXC_Z_OUT_DQ',
            'K1:PEM-MAG_EYC_BOOTH_EYC_X_OUT_DQ',
            'K1:PEM-MAG_EYC_BOOTH_EYC_Y_OUT_DQ',
            'K1:PEM-MAG_EYC_BOOTH_EYC_Z_OUT_DQ']
for i in range(9):
   data = TimeSeries.read(source, channels[i], format='gwf.lalframe')    
}}}
21.863秒かかった。次に、TimeSeriesDictでまとめて読みだすと、
{{{#!highlight python
from gwpy.timeseries import TimeSeries, TimeSeriesDict
sources = ['/data/full/12677/K-K1_C-1267797600-32.gwf']
channels = ['K1:PEM-MAG_BS_BOOTH_BS_X_OUT_DQ',
            'K1:PEM-MAG_BS_BOOTH_BS_Y_OUT_DQ',
            'K1:PEM-MAG_BS_BOOTH_BS_Z_OUT_DQ',
            'K1:PEM-MAG_EXC_BOOTH_EXC_X_OUT_DQ',
            'K1:PEM-MAG_EXC_BOOTH_EXC_Y_OUT_DQ',
            'K1:PEM-MAG_EXC_BOOTH_EXC_Z_OUT_DQ',
            'K1:PEM-MAG_EYC_BOOTH_EYC_X_OUT_DQ',
            'K1:PEM-MAG_EYC_BOOTH_EYC_Y_OUT_DQ',
            'K1:PEM-MAG_EYC_BOOTH_EYC_Z_OUT_DQ']
data = TimeSeriesDict.read(sources, channels, format='gwf.lalframe')
}}}
5.166秒と短縮された。

次に、読むファイル数を10個にしてみる。31.008秒かかった。TimeSeriesDict.read の引数に ''nproc=コア数'' を入れると並列化され、より速くなる。
{{{#!highlight python
data = TimeSeriesDict.read(sources, channels, format='gwf.lalframe', nproc=10)
}}}
結果は5.485秒。ちなみに、ファイル数が1のときにはnprocを増やしても変わらなかった。


-----

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

pythonはループ処理が非常に遅いことで有名だが、積極的にリスト内包表記で置き換えることでかなり高速化される。<<BR>>
参考：[[https://qiita.com/kokorinosoba/items/eb72dac6b68fccbac04d|リスト内包表記とif文を組み合わせるとき]]

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


-----

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

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

{{{#!highlight python
import shelve
from tqdm import tqdm
from gwpy.timeseries import TimeSeries, TimeSeriesDict, StateVector

def SingleEvent(TriggerTime): #eventの時刻を引数にする
   data  = TimeSeriesDict.read(...)
   state = StateVector.read(...)
   return data, state

def main():
   ...
   data  = []
   state = [] 

   for t in tqdm(TrigTimes): #予め各eventの時刻をリストにしておく
        d, s = SingleEvent(t)
        data.append(d)
        state.append(s)

   f = shelve.open('./EventList.slv')
   f['data']  = data
   f['state'] = state
   f['TrigTimes'] = TrigTimes
   f.close()
}}}
このshelveを読むには、以下のようなコードを書けばよい。

{{{#!highlight python
import shelve
f = shelve.open('./EventList.slv')
for data, state, t in zip(f['data'], f['state'], f['TrigTime']):
   ...#後は好きなように解析
}}}
「ちょっと10 eventの時系列波形見たいなー」といった場合は、インタープリンタで
{{{
>>> import shelve 
>>> f = shelve.open('./EventList.slv') 
>>> f['data'][10].plot().show()
}}}
と打つだけですぐに表示される。<<BR>>
爆速になるので超絶おススメ。

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


-----

== 時系列データをDown samplingする / スペクトログラムの時点でBand limit ==
見たい周波数帯に比べてサンプリング周波数が無駄に大きいときには、あらかじめTimeSeriesに
[[https://gwpy.github.io/docs/stable/api/gwpy.timeseries.TimeSeries.html#gwpy.timeseries.TimeSeries.resample|.resample()]] 
をかけて計算コストを下げると良い。

[[https://gwpy.github.io/docs/stable/examples/frequencyseries/percentiles.html|スペクトルのパーセンタイルを計算するとき]]には、
まず[[https://gwpy.github.io/docs/stable/api/gwpy.timeseries.TimeSeries.html#gwpy.timeseries.TimeSeries.spectrogram2|.spectrogram2()]]を使うが、
その時点で[[https://gwpy.github.io/docs/stable/api/gwpy.spectrogram.Spectrogram.html#gwpy.spectrogram.Spectrogram.crop_frequencies|crop_frequencies()]]
でband limitしておくとその後の処理が軽くなる。

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


-----

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

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

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

まずgwpyの場合、
{{{#!highlight python
plot = data.plot().show()
}}}
{{attachment:plot_gwpy.png}}

11.86秒。次にpyROOTの場合、

{{{#!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')
}}}

{{attachment:plot_root.png}}

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

{{attachment:plot.gw2.png}}

-----





== Cython ==
Cythonについて調べると、setup.py なるファイルを作れという記事が多くヒットするが、現代ではこれは必要なくなった。([[https://se.miyabikno-jobs.com/cython-compile-cythonize/|ここ]]を参照)


=== Cythonのコンパイル -> 実行を一発でできるようにする ===
Cythonを使って高速化したはいいけど、毎回コンパイルするのは面倒だ。複数のコードになってしまうのは嫌だ。という場合の小技。<<BR>>
home directoryに ''.cython.sh'' という名前で以下を書き込んだシェルスクリプトを作る。 
{{{#!highlight sh
#!/usr/bin/bash
set -e
trap 'echo compile error' ERR
code=$1
cythonize -i $code > .cython_log.txt
rm ${code%.pyx}.c
echo 'import '${code%.pyx}  > .a.py
echo ${code%.pyx}'.main()' >> .a.py
python .a.py "${@:2}"
}}}
また、 ''~/.bashrc'' に次の一行を書き足す。
{{{#!highlight bash
alias cythonsh='sh ~/.cython.sh'
}}}
これで ''.bashrc'' を有効にすれば、ターミナルで
{{{
> cythonsh (cythonコード).pyx 引数
}}}
と打つだけでコンパイルから実行まで行われる。<<BR>>
2, 3行目は、コンパイルに失敗した時にスクリプトが止まるようにするおまじない。<<BR>>
6行目は .c ファイルを残さないようにしている。<<BR>>
9行目の末尾の "${@:2}" は、引数をpythonに渡すおまじない。


-----


== C++でPythonモジュールを作る ==
以下のサイトがわかりやすい
[[https://cpp-learning.com/python_c_api_step1/]]

-----

== リンク集 ==
 * [[https://qiita.com/shaka/items/f180ae4dc945dc7b9066|Pythonを速くしたいときにやったこと]]

C言語は使えるけど最近Pythonを始めた」人向けの記事
 * [[https://qiita.com/jabberwocky0139/items/c3620fb2f011f20a633b|NumPyによる数値計算の高速化 : 基礎]] 
 * [[https://qiita.com/jabberwocky0139/items/a9751d11caa64bc19226|NumPy・SciPyを用いた数値計算の高速化 : 応用その1]]
 * [[http://qiita.com/jabberwocky0139/items/26451d7942777d0001f1|NumPy・SciPyを用いた数値計算の高速化 : 応用その2]]