そうだPythonを勉強しよう(NumPy編)

Python

http://d.hatena.ne.jp/yukichanko/20110215/1297784078

あいもかわらず、この続き。

とりあえず、NumPyを使う

Pure Python版で、

  ey = [0.0] * (Const.NX+1)
  hz = [0.0] * Const.NX

としている定義を単にNumPy Arrayに置き換えます。

  import numpy as np
  ey = np.array([0.0] * (Const.NX+1))
  hz = np.array([0.0] * Const.NX)

メインのループはそのままです。

def calc_ey(ey, hz, dx, dt):
  """ 電界(Ey)更新 """
  for i in range(1, Const.NX):
    ey[i] = ey[i] - dt / (Const.PERMITTIVITY * dx[i]) * (hz[i] - hz[i-1])

def calc_hz(hz, ey, dx, dt):
  """ 磁界(Hz)更新 """
  for i in range(0, Const.NX):
    hz[i] = hz[i] - dt / (Const.PERMEABILITY * dx[i]) * (ey[i+1] - ey[i])

とすると、スピードはどうなるか。。。当然スピードはあがりません。
結果は、

言語 CPU Time Memory
C 0.91sec(1.0) 416KB
Python(Pure List) 72.9sec(80.1) 3.29MB
Python(NumPy Loop) 476sec(523) 7.86MB

。。。えっ( ꒪⌓꒪) いや、予想していたとは言え、Pure ListとNumPy Loop版でこんなに遅くなるとは!!

正直、予想外。Pure Listと同等だと思っていました。話の流れとしては、単にリストをNumPy Arrayにしてもダメ。Pythonのループが遅いから、NumPy内のループに変えましょう、と言うのが話の流れなわけだけど、ここまで遅くなるというのは、Pythonのループの問題だけじゃないね。

PythonのループとNumPyのスカラアクセスが相性悪いのか、それともわざと非効率なインプリにして、ユーザーにアラートを与えているのか？

Pythonのループをなくす

ここまで遅いと、さすがにループの中でNumPyは使えません。このループをなくしましよう。どうするか。。。スライスを使用して、ベクトル演算に変えます。

def calc_ey(ey, hz, dx, dt):
  """ 電界(Ey)更新 """
  ey[1:Const.NX] = ey[1:Const.NX] \
      - dt / (Const.PERMITTIVITY * dx[1:Const.NX]) * (hz[1:Const.NX] - hz[0:Const.NX-1])

def calc_hz(hz, ey, dx, dt):
  """ 磁界(Hz)更新 """
  hz[0:Const.NX] = hz[0:Const.NX] \
      - dt / (Const.PERMEABILITY * dx[0:Const.NX]) * (ey[1:Const.NX+1] - ey[0:Const.NX])

はてさて、実力は。。。

言語 CPU Time Memory
C 0.91sec(1.0) 416KB
Python(Pure List) 72.9sec(80.1) 3.29MB
Python(NumPy Loop) 476sec(523) 7.86MB
Python(NumPy Vector) 5.53sec(6.07) 7.90MB

うん。だいぶん使えるぐらいになりました。まだ、Cとの差はありますが、それでもPythonのフレキシビリティを考えると、許容できるぐらいの速度差にはなったのではないでしょうか？

(NumPyを使ってもCほどは速くならないよ、と言う話を"Python Scripting for Computational Science"の中で見たような気がするのですが、線を引いてなかったので、どこかわからず。見つけたら、参照してみます。)

次は

SWIGを使って、電磁界更新部をCのコードに置き換えます。