f2pyを用いてnumpyのデータを解析

numpyはデータの取り扱いは便利だが、Pythonはforループが遅いので微分などをうまく行う方法がない。f2pyで微分の部分のみfortranで書いてやれば、高速化できて良い。微分の操作などは多くの解析で共通なので、一度書いておくと使いまわせて良い。

f2pyの公式マニュアルはこちらにあるので参照されたい。

fortranコード

サブルーチンを用いて関数を書くのが基本と思われる。
基本事項は以下

intentは必ず書く
返り値はintent(out)で管理

例えば、「配列に1を足して返す関数」は以下のように書ける。ファイル名はf2py_test.f90などとしておく

subroutine add(qq,ix,qq_add)
  implicit none
  
  integer :: i
  integer, intent(in) :: ix
  real(8), dimension(ix), intent(in) :: qq
  real(8), dimension(ix), intent(out) :: qq_add
  
  do i = 1,ix
    qq_add(i) = qq(i) + 1
  enddo

  return
end subroutine add

コンパイル

f2py --fcompiler=gfortran -m f2py_test -c --f90flags='-O3' f2py_test.f90

とする。オプション -c は必ずつけて最後にコンパイルするファイル名を書く。それぞれのオプションは以下

--fcompiler=: fotranコンパイラ指定
-m: Pythonで用いるときのモジュール名(ファイル名と同じでなくて良い)
--f90flags=: fortranコンパイラで使うオプション

Pythonからfortranコードの利用

Pythonからf2pyを利用する時は、以下のようにする。

import f2py_test # f2pyでコンパイルしたモジュールをインポート
import numpy as np
ix = 5
qq = np.zeros(ix) # 適当にnumpy配列を定義
qq_add = f2py_test.add(qq,ix) # f2pyで定義した関数を呼び出す

fortranコードで定義したqq_addはintent(out)としてあるので、戻り値と判断され、関数の引数には含めない。
一方、上記のコードでは配列の大きさixも関数の引数として含めたがf2pyでは、配列の大きさは自動的に判定してオプション変数とした上で、自動的にlen(qq)を代入してくれる。よって、

qq_add = f2py_test.add(qq)

としても同様に結果になる。

触れないが触れたいこと

多次元・複数の配列でも基本的に方法は同じ
変数の形が呼び出しで違ってもf2pyがなんとかしてくれる
一つのファイルの中に複数のサブルーチンを書いても問題ない

終わりに

numpyの配列を解析する時に、f2pyを使う心の障壁を低くして解析を高速にしたいと思う