Fortranでよく使うdirect accessバイナリをPythonで読み書きするためのライブラリFORBIT

気象・海洋のコミュニティでは、Fortranで計算された結果をバイナリファイルに出力し、それをPythonで読み込んでグラフ化するワークフローが非常に一般的です。
出力ファイルのフォーマットは、最近はNetCDFがよく使われていますが、Fortranが標準でサポートする　form="unformatted", access="direct"　のような、Direct access型のバイナリ(以後、乱暴ですが単にバイナリと書いたらこのフォーマットを指します)は今でも健在です。

バイナリは、Fortranではたとえば次のように読みます

open(newunit=unit         ,
     file   ='sample.bin' ,
     action ='read'       ,
     form   ='unformatted',
     access ='direct'     ,
     recl   =kind*nx*ny*nz)

do t = 1, nmonth, 12
    read(unit,rec=t) foo(1:nx,1:ny,1:nz)
    ! any processes
end do

close(unit)

ここでは、monthlyデータで特定の月の場だけを取り出す雑な例です。
1 record=1変数や1時刻として考えれば、recordの更新が入出力変数の変更やデータ時刻の更新などとそのまま結びつくのでシンプルです。

これと同じことをPythonでNumPyを使って書くと、たぶん次のようになります(いくつか未定義のままの変数がありますがご了承ください)

import numpy as np

size  = nx * ny * nz
shape = [nz,ny,nx]
recl  = size * np.dtype(np.float32).itemsize
foo   = np.empty(shape)
file  = open('sample.bin', 'rb')

for t in range(0, nmonth, 12):
    file.seek(t * recl)
    foo[:,:,:] = np.fromfile(file,dtype=np.float32,count=size).reshape(shape)
    # any processes
    
file.close()

Fortranでは、open時に色々な設定を書く代わりに、loop内での処理が簡単です。また、取得される配列の形状は与えたfooに従います。
Pythonでは、open時の設定入力は最低限で、readごとに読み込みデータ数やオフセットを指定したり、reshapeしたりします。

これらの違いは双方にメリットとデメリットがありますが、Fortranの頭の使い方でPythonを書くと少々苦労します。
個人的には、seekに与えるスキップ先をbyteで管理する点が、recordで管理するFortranと違うのでドキドキします。

Fortranの頭でPythonを書きたい

私はFortranがメイン言語でPythonはグラフ描画専用なので、あまりPython的なやり方には慣れません。
そこで、PythonでのバイナリIOを可能な限りFortranの気持ちで実行するためのライブラリforbitを開発しました。

GitHub

PyPI

forbitは、内部ではCythonを経由して実際にFortranの関数を呼び出します。
Fortranが得意なファイルフォーマットはFortranで読み書きする、が実装の肝です。
$ pip install forbit でインストールできます。
これを使うと、たとえば先ほどの例のmonthlyデータの読み込みは次のように書かれます

import numpy as np
import forbit

shape   = [nz,ny,nx]
kind    = 4
record  = 1
recstep = 12
endian  = 'little_endian'
foo     = np.empty(shape, dtype=np.float32)

file = forbit.open('sample.bin',
                   action ='read' ,
                   shape  =shape  ,
                   kind   =kind   ,
                   record =record ,
                   recstep=recstep,
                   endian =endian )

for t in range(nyear):
    foo[:,:,:] = file.read()
    # any processes

file.close()

Fortran的な、open時に全部指定してreadはシンプルな形になりました。
また、Fortranと同じようにrecord orientedなので、頭の使い方を切り替える必要がありません。
action='write' を与えればwriteも似たような書き方で動かせます。

open() のoptional引数に、recl と dtype があります。
recl は、指定しなければ shape と kind から計算されます。指定する場合には、kind*prod(shape) で計算される値以上の整数が必要です。
dtype は、デフォルトでは 'real' が与えられ、浮動小数型のIOを実行します。整数を使いたい場合には 'integer' を指定します。

入出力はndarray限定です。
ライブラリを動かすための最低限のメタデータしか持たないので、xarrayや独自classを用いたwrapが容易にできる設計になっています。
ユーザの用途や趣味に最適化された運用が可能で、forbitはその基盤を担う想定です。

データ型

実際にやりとりするデータ型は、dtype と kind によって決定されます。
執筆時点で対応しているデータ型は次です

dtype	kind	Fortran	NumPy
'real'	4	real32	numpy.float32
'real'	8	real64	numpy.float64
'integer'	2	int16	numpy.int16
'integer'	4	int32	numpy.int32
'integer'	8	int64	numpy.int64

最初から dtype=numpy.float32 のような与え方をすれば良いのでは…という検討もしましたが、Fortranの世界はkind parameterが記述していると思っているので、kind parameterを中心に据えた世界観を優先しました。
また、私の周りではバイナリファイルの大半がreal型なので、わざわざ指定しなくてもデフォルトの挙動でrealが指定されてほしかったのも大きな理由です。

Record

read/writeの実行ごとに指定されるrecordは、Cython classが保持している forbit.__record に従います。
この値は、forbit.open 時に与えた record 引数と recstep によって決定され、read/writeするごとにrecord=record+recstepのような計算が起こります。
すなわち、record=1, recstep=12 で初期化した場合、1回目のreadでは1 record目が、2回目、3回目…はrecord=13, 25, …となります。
それぞれの時刻に3変数(たとえばu, v, t)が書かれているバイナリを読む時には、recstep=3 でひとつの変数の時系列を取り出せます。

このルールから外れた複雑なrecord指定をしたい場合には、forbit.reset_record() で次のrecordを直接指定できます。
引数は newRecord と increment があり、newRecord に値を与えた場合には与えた値がそのまま次のrecordとして使用されます。
increment に値を与えると、現在のrecordにその値が加算され、使用されます。
両方が指定された場合は、newRecord が優先され、increment は無視されます。

また、次のアクセス先を知りたい時には、forbit.get_record() でrecordを取得できます。

Endian

forbit.open() に与える endian 引数は、'little_endian', 'big_endian', 'native' に対応しています。
Fortranの open() にある convert 引数に、この文字列をそのまま与える設計です。
バイナリはendianの処理が面倒ですが、この引数を適切に与えることで正しくデータを扱うことが可能です。
ただし、Fortranでは convert 引数は規格に定義されておらず、コンパイラ依存があるので、後述するようにコンパイラをデフォルトから変更する場合には注意が必要です。

内部処理

※オタク以外は読み飛ばしてください
内部処理は至ってシンプルです。Fortranで書かれた関数をCythonでwrapしています。

forbit.open() の kind および dtype 引数で与えられたデータ型のやり取りをするため、open() では概念的には次のような処理を行います

if (dtype == 'real' and kind == 4):
    if (action == 'read'):
        self.read  = self.read_r4      # for reading real32
        self.write = self.write_err    # for printing error message
    elif (action == 'write'):
        self.read  = self.read_err     # for printng error message
        self.write = self.write_r4     # for writing real32
    elif (action == 'readwrite'):
        self.read  = self.read_r4      # for reading real32
        self.write = self.write_r4     # for writing real32
elif (dtype == 'real' and kind == 8):
    ...

これにより、コンストラクタで read() と write() が適切に振り分けられるので、読み書き実行時の条件分岐が最低限になり、若干の高速化が期待できます。
また、ユーザはopen時にだけデータ型を考えれば、以後は read_r4 のような名前ではなく、共通の read を呼び出すだけなので、データ型と関数名の対応に関する不安から解放されます。
Fortranのgeneric interfaceのようなイメージですね。
なお、上記のコードはあくまで概念の例で、実際にはlist型とbit操作を活用したシンプルな選択アルゴリズムを組んでいます。何度も分岐を書くのは趣味じゃないので。

コンパイラ

pipでインストールした場合には、デフォルトでgfortranとgccが使われます。
GitHubからcloneすれば、src/Makefile を編集することで好みのコンパイラへの変更が可能です。

ただし、Fortranの open 文でendian指定に convert 引数を使用するので、これをサポートしたコンパイラが必要です。
また、recl はbyteを単位とした指定になっているので、ifxなどの格子数で指定するコンパイラに変更する場合にも注意しなければなりません。ifxを使う場合には、-assume byterecl をコンパイラオプションに加える必要があります。

速度

私の環境でのベンチマークによると、numpy.fromfile() および ndarray.tofile() に比べforbitが若干速く読み書きしました。
ただし、ストレージのHDD/SSDの別やコンパイラ、NumPyのバージョンに依存する要素が大きいと思われるので、明言しません。
forbitの性能はgfortranによるFortran IOに依存します。
少なくとも、このライブラリの使用が速度面で足を引っ張る心配はないはずです。

最後に

初めての自作ライブラリPyPI公開と、初めてのQiita執筆でした。
バグやご意見等あれば、コメントやIssueにて。
緊張したー