Calculation of grid cell center value by Python
Pythonで格子のセルの中心の値を算出する方法について検討した。結論は，4配列一括演算が最速で，全てリスト演算するとその5倍程の時間がかかった。しかし，リストはわかりやすいのでこれでもよいかもしれない。
1 introduction
数値計算をやっていると格子の中心の値がほしいことがよくある。具体的に説明する。以下の図でノード(青い点)ごとに値をもっている。この4点のノードの値を平均して中心のセルの中心(赤い点)の値がほしい。
numpyとかscipy，pandasのメソッドを調べてみてもそんな便利な関数が見当たらなかった。
例えば，以下ではFortranでそういう関数を自作してf2pyで利用しているという投稿があった（なお，上記の図は以下から引用している）。
python – Resample grid from center coordinates to external (i.e. corner) coordinates – Stack Overflowhttp://stackoverflow.com/questions/22666937/resample-grid-from-center-coordinates-to-external-i-e-corner-coordinates
自力で実装するしかないようだ。実装する際に考えられる方法としては以下がある。
- 素直に2重ループにより4点を平均。
- 4配列の一括演算
それぞれ100万ノードの配列を計算を試してどちらが早いか検討した。以下で実装と検討結果を説明する。なお，結論をいうと2.の4配列の一括演算が最も早かった。ただし，4配列一括演算の約5倍の速度となるが，可読性やメモリの占有を考えると全てリストで演算してもよいと思った。
2 2重ループ
2重ループで実装する場合，numpyの関数の組み合わせていくつかパターンがある。違いは以下だ。
- ループするときにlistで回すかndarrayで回すか
- 平均するときにnumpy.meanを使うかlist演算(sumを使って合計を出しそれを割って算出)するか
これを試すために以下の表のケースのプログラムを用意した。
表 8.1: 2重ループのケース
プログラム名
|
ループ変数
|
平均操作
|
center-np.py
|
ndarray
|
numpy
|
center-loop-array.py
|
ndarray
|
list
|
center-loop-npmen.py
|
list
|
numpy
|
center-loop-alllist.py
|
list
|
list
|
2.1 center-np.py: ループ: ndarray, 平均: np.mean
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
| #!/usr/bin/env python2.7# coding: utf-8# (File name: center-loop.py)# Author: SENOO, Ken# (Last update: 2014-09-23T18:30+09:00)# License: MITimport numpy as npLENG=1000000node=np.arange(LENG).reshape(-1,np.sqrt(LENG))nrow, ncol=node.shapecenter=[np.mean((node[row][col], node[row][col+1],node[row+1][col], node[row+1][col+1]), dtype=np.float32) for row in range(nrow-1) for col in range(ncol-1)]center=np.array(center).reshape(nrow-1, ncol-1) |
14行目で演算している。左上から順に以下のようにz字型に値を列挙し，平均している
UL, UR, LL, LR
ここで1文字目のU/LはUpper(上)/Lower(下)を意味し，2文字目のL/RはLeft(左)/Right(右)を意味する。
なお，np.mean関数を使うときは精度を上げておく。そうしないと，小数点が桁落ちして正しい値が求まらない。
これ以降のプログラムでは12行目に以下のコードでリストに変更するか，
node=node.tolist()
15行目のループでnp.meanをsum /4としているかの違いとなっている。
2.2 center-loop-npmean.py：ループ:list, 平均: np.mean
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
| #!/usr/bin/env python2.7# coding: utf-8# (File name: center-loop-npmean.py)# Author: SENOO, Ken# (Last update: 2014-09-23T18:30+09:00)# License: MITimport numpy as npLENG=1000000node=np.arange(LENG).reshape(-1,np.sqrt(LENG))nrow, ncol=node.shapenode=node.tolist()center=[np.mean((node[row][col], node[row][col+1],node[row+1][col], node[row+1][col+1]), dtype=np.float32) for row in range(nrow-1) for col in range(ncol-1)]center=np.array(center).reshape(nrow-1, ncol-1) |
2.3 center-loop-array.py: ループ: ndarray, 平均: list
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
| #!/usr/bin/env python2.7# coding: utf-8# (File name: center-loop-array.py)# Author: SENOO, Ken# (Last update: 2014-09-23T19:24+09:00)# License: MITimport numpy as npLENG=1000000node=np.arange(LENG).reshape(-1,np.sqrt(LENG))nrow, ncol=node.shapecenter=[sum((node[row][col], node[row][col+1],node[row+1][col], node[row+1][col+1]))/4.0 for row in range(nrow-1) for col in range(ncol-1)]center=np.array(center).reshape(nrow-1, ncol-1) |
2.4 center-list.py: ループ: list，平均: list
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
| #!/usr/bin/env python2.7# coding: utf-8# (File name: center-loop-alllist.py)# Author: SENOO, Ken# (Last update: 2014-09-23T19:11+09:00)# License: MITimport numpy as npLENG=1000000node=np.arange(LENG).reshape(-1,np.sqrt(LENG))nrow, ncol=node.shapenode=node.tolist()center=[sum((node[row][col], node[row][col+1],node[row+1][col], node[row+1][col+1]))/4.0 for row in range(nrow-1) for col in range(ncol-1)]center=np.array(center).reshape(nrow-1, ncol-1) |
3 4配列一括演算
2重ループを使わないもう一つの方法を説明する。この方法は，右上，左下，右下の点を，基準とする左上の点の位置にずらした配列を作り，合計四つの配列を一括して演算する。具体的には以下の手順を取る。
- 4点が重なるように右，左下，右下の点から始まる配列を作る。
- 作成した配列の端っこの部分をnanで埋める。
- この四つの配列を一括で演算して平均化する。
- できた配列の最終行，最終列を除外して目的の配列を得る。
このロジックに基づいたプログラムを以下に示した。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
| #!/usr/bin/env python2.7# coding: utf-8# (File name: center4array.py)# Author: SENOO, Ken# (Last update: 2014-09-23T18:30+09:00)# License: MITimport numpy as npLENG=1000000node=np.arange(LENG).reshape(-1,np.sqrt(LENG))nrow, ncol=node.shapeur=np.column_stack([node[:,1:], np.repeat(np.nan, nrow)])ll=np.vstack([node[1:,:], np.repeat(np.nan, ncol)])lr=np.column_stack([node[1:,1:], np.repeat(np.nan, nrow-1)])lr=np.vstack([lr, np.repeat(np.nan, ncol)])center=np.nanmean(np.dstack([node, ur, ll, lr]), axis=2, dtype=np.float32)[:-1,:-1] |
4 Discussion
用意した合計五つのプログラムを以下のようにtimeコマンドを実行時間を測定した。
time ./center-4array.py
集計結果を以下の表にまとめた。
表 8.2: 格子の中心値を算出するプログラムの実行結果
ファイル名
|
説明
|
所要時間[s]
|
center-4array.pyとの比
|
center-np.py
|
2重ループ。ループ: ndarray，平均: np.mean
|
13.57
|
96.93
|
center-loop-npmean.py
|
2重ループ。ループ:list, 平均: np.mean
|
12.26
|
87.57
|
center-loop-array.py
|
2重ループ。ループ: ndarray, 平均: list
|
5.43
|
38.79
|
center-list.py
|
2重ループ。ループ: list，平均: list
|
0.62
|
4.43
|
center-4array.py
|
4配列一括演算
|
0.14
|
1.00
|
この結果から，4配列を一括演算する方法が最も高速となった。次点は全ての演算をlistで行うプログラムとなった。4配列での一括演算は速度こそ高速であるが，占有メモリとコードの可読性の点ではあまりよくないと思った。そのため，こだわりがなければ全てリストで演算してもよいかもしれない。
center-np.pyとcente-loop-npmean.pyの比較から，ループの最も内側でnumpy.meanを使うと使わない場合の約2.5倍遅くなっている。また，center-list.pyとcenter-loop-array.pyの比較からループを回すときにlistでなくndarrayのまま行うとlistの約10倍程遅くなっている。
このことから，numpyを使うときは可能な限り，配列演算を行い，それが無理でやむなくループを使う場合は，全てリストで演算すべきだといえる。中途半端にnumpyの要素(関数，変数)をループに入れると計算速度が極端に遅くなる。
ここ半年くらいの悩みが解決してすっきりした。
