pythonとfortranのI/O関連まとめ

# はじめに

python<->fortranの入出力について、インターネット上でまとまった記述が発見しにくいので、わかっている範囲でまとめることにした。

余談だが、python同士のデータ入出力には`np.savez`と`np.load`がおすすめです。

# pythonの出力をfortranで読み込む
## numpyを用いる方法
以下のような手順でデータを出力する
1. numpy配列を準備
2. `reshape`でfortranに合うように変換
3. `astype`でバイナリに変換
4. `tofile`で出力

3次元の倍精度実数をlittle endianで出力するpythonのコードの例

```python
import numpy as np

ix, jx, kx = 3, 4, 5
endian='<'
a = np.ones((ix,jx,kx))
a.reshape([ix*jx*kx],order='F').astype(endian+'d').tofile('test.dat')
```
これをfortranで読み込むには以下のようにする。使っている計算機がlittle endianだとする。fortranでの`open`には`access='stream'`が必要となる

```fortran
program test
  implicit none

  integer :: idf
  integer, parameter :: ix = 3, jx = 4, kx = 5
  real(KIND(0.d0)), dimension(ix,jx,kx) :: a

  open(newunit=idf, file='test.dat',form='unformatted',access='stream')
  read(idf) a
  close(idf)
  
  stop
end program test
```

big endianのデータをやりとりしたいときはpythonのコードで

```python
endian='>'
```

と変えれば良い。この`endian`変数は`astype`でバイナリへの変換方法を制御している

fortranコードでは

```fortran
  open(newunit=idf, file='test.dat',form='unformatted',access='stream',convert='big_endian')
```

とする。単精度で出力・入力したいときはpythonコードを以下のように書き換える

```python
import numpy as np

ix, jx, kx = 3, 4, 5
endian='<'
a = np.ones((ix,jx,kx),dtype=np.float32)
a.reshape([ix*jx*kx],order='F').astype(endian+'f').tofile('test.dat')
```
`np.ones`と`astype`の部分を単精度を扱えるように変更した。

対応するfortranコードは以下のようになる

```fortran
program test
  implicit none

  integer :: idf
  integer, parameter :: ix = 3, jx = 4, kx = 5
  real(KIND(0.e0)), dimension(ix,jx,kx) :: a

  open(newunit=idf, file='test.dat',form='unformatted',access='stream')
  read(idf) a
  close(idf)
  
  stop
end program test
```

`real(KIND(0.e0))`の部分のみを変更した。

複数の配列を一つのファイルに出力したい場合は以下の手順を踏む
1. 複数の配列を`np.array`でまとめる
2. `reshape`で１次元配列とする。
2. `T`でインデックスの位置を入れ替える。

```python
import numpy as np

ix, jx, kx = 3, 4, 5
endian='<'
a = np.ones((ix,jx,kx),dtype=np.float32)
b = np.ones((ix,jx,kx),dtype=np.float32)*2
c = np.ones((ix,jx,kx),dtype=np.float32)*3
np.array([a,b,c]).T.reshape([3*ix*jx*kx],order='F').astype(endian+'f').tofile('test.dat')
```

読み込むfortranコードは

```fortran
program test
  implicit none

  integer :: i,j,k
  integer :: idf
  integer, parameter :: ix = 3, jx = 4, kx = 5
  real(KIND(0.e0)), dimension(ix,jx,kx) :: a,b,c

  open(newunit=idf, file='test.dat',form='unformatted',access='stream')
  read(idf) a,b,c
  close(idf)
  
  stop
end program test
```

とすれば良い

一方、大きさの違う配列を出力するときは
1. `reshape`を用いて1次元配列とする
2. `np.concatenate`を用いて配列をまとめる
という手順をふめば良い

```python
pythonコードは以下のようにする

import numpy as np

ix, jx, kx = 3, 4, 5
lx, mx = 2, 4
endian='<'
a = np.ones((ix,jx,kx))
b = np.ones((lx,mx))*2

a1 = a.reshape([ix*jx*kx],order='F')
b1 = b.reshape([lx*mx],order='F')
np.concatenate([a1,b1]).astype(endian+'d').tofile('test.dat')
```

fortranから読み込むときはこれまでと変わらないが以下のようにする。

```fortran
program test
  implicit none

  integer :: idf
  integer, parameter :: ix = 3, jx = 4, kx = 5
  integer, parameter :: lx = 2, mx = 4
  real(KIND(0.d0)), dimension(ix,jx,kx) :: a
  real(KIND(0.d0)), dimension(lx,mx) :: b

  open(newunit=idf, file='test.dat',form='unformatted',access='stream')
  read(idf) a,b
  close(idf)
  
  stop
end program test
```

## scipy.io.FortranFileを用いる方法

上記の方法は、pythonの方のコードは簡単だが、fortranの方で`access='stream'`を使わないと自分でデータのヘッダをうまく編集してやらないといけなくなる。
この部分をうまく調整してくれるのが`scipy.io.FortranFile`である。

Pythonコードは以下のようにする

```python
import numpy as np
from scipy.io import FortranFile

ix, jx, kx = 3, 4, 5
endian='<'
a = np.ones((ix,jx,kx))
with FortranFile('test.dat','w') as f:
    f.write_record(a)
```

読み込みのためのFortranコードは

```fortran
program test
  implicit none

  integer :: i,j,k
  integer :: idf
  integer, parameter :: ix = 3, jx = 4, kx = 5
  real(KIND(0.d0)), dimension(ix,jx,kx) :: a

  open(newunit=idf, file='test.dat',form='unformatted')
  read(idf) a
  close(idf)
  
  stop
end program test
```
とすれば良い。

複数の配列を出力するときはNumpyを用いる時よりも容易で

```python
a = np.ones((ix,jx,kx))
b = np.ones((ix,jx,kx))*2
with FortranFile('test.dat','w') as f:
    f.write_record(a,b)
```

と`write_record`の引数に付け足すのみで良い。
fortranの読み込みコードも

```fortran
program test
  implicit none

  integer :: i,j,k
  integer :: idf
  integer, parameter :: ix = 3, jx = 4, kx = 5
  real(KIND(0.d0)), dimension(ix,jx,kx) :: a,b

  open(newunit=idf, file='test.dat',form='unformatted')
  read(idf) a,b
  close(idf)
  
  stop
end program test
```
と読み込む配列を横に並べればいいだけである。

とはいえ、[scipyのウェブサイト](https://docs.scipy.org/doc/scipy/reference/generated/scipy.io.FortranFile.html)を見ると

> Since this is a non-standard file format, whose contents depend on the compiler and the endianness of the machine, caution is advised. Files from gfortran 4.8.0 and gfortran 4.1.2 on x86_64 are known to work.

> Consider using Fortran direct-access files or files from the newer Stream I/O, which can be easily read by numpy.fromfile.

と書かれている。fortranのフォーマットはコンパイラ依存が大きいらしくnumpyを用いた方法が推奨されるようである。

# fortranの出力をpythonで読み込む
## numpyを用いる方法

numpyでfortranのデータを読み込むためには`access='stream'`で出力するのがおすすめである(ヘッダなどを自分で調節すればいいので必須ではない)。

fortranで配列を定義して出力する

```fortran
program test
  implicit none

  integer :: i,j,k
  integer :: idf
  integer, parameter :: ix = 3, jx = 4, kx = 5
  real(KIND(0.d0)), dimension(ix,jx,kx) :: a

  a = 1.d0
  
  open(newunit=idf, file='test.dat',form='unformatted',access='stream')
  write(idf) a
  close(idf)
  
  stop
end program test
```

pythonで読み込むには最も単純には、`np.fromfile` を用いれば良い。

```python
import numpy as np

ix, jx, kx = 3, 4, 5
endian='<'
with open('test.dat','rb') as f:
    a = np.fromfile(f,endian+'d',ix*jx*kx)
    a = a.reshape((ix,jx,kx),order='F')
```

`np.fromfile` を用いる場合でも`np.dtype`を使うとより汎用性の高い読み込みが可能になる。例えば以下のように行う

```python
import numpy as np

ix, jx, kx = 3, 4, 5
endian='<'
dtype = np.dtype([('a',endian+str(ix*jx*kx)+'d')])
with open('test.dat','rb') as f:
    a = np.fromfile(f,dtype=dtype)['a']
a = a.reshape((ix,jx,kx),order='F')
```

複数の配列を取り扱うときは`np.dtype`を用いる方が便利である。
まずfortranで以下のように複数の配列を一つのファイルに出力する。

```fortran
program test
  implicit none

  integer :: i,j,k
  integer :: idf
  integer, parameter :: ix = 3, jx = 4, kx = 5
  real(KIND(0.d0)), dimension(ix,jx,kx) :: a,b

  a = 1.d0
  b = 2.d0
  
  open(newunit=idf, file='test.dat',form='unformatted',access='stream')
  write(idf) a,b
  close(idf)
  
  stop
end program test
```

これをpythonで読み込むには以下のようにする

```python
import numpy as np

ix, jx, kx = 3, 4, 5
endian='<'
dtype = np.dtype([('a',endian+str(ix*jx*kx)+'d'),('b',endian+str(ix*jx*kx)+'d')])
with open('test.dat','rb') as f:
    data = np.fromfile(f,dtype=dtype)
a = data['a'].reshape((ix,jx,kx),order='F')
b = data['b'].reshape((ix,jx,kx),order='F')
```

## scipy.io.FortranFileを用いる方法
`access='stream'`を用いていないfortranデータを読み込むためには`scipy.io.FortranFile`を用いるのが容易である。

fortranの出力は以下のようにする

```fortran
program test
  implicit none

  integer :: i,j,k
  integer :: idf
  integer, parameter :: ix = 3, jx = 4, kx = 5
  real(KIND(0.d0)), dimension(ix,jx,kx) :: a

  a = 1.d0
  
  open(newunit=idf, file='test.dat',form='unformatted')
  write(idf) a
  close(idf)
  
  stop
end program test
```

これを読み込むpythonコードは以下のようにする

```python
import numpy as np
from scipy.io import FortranFile

ix, jx, kx = 3, 4, 5
endian='<'
a = np.ones((ix,jx,kx))
b = np.ones((ix,jx,kx))*2
f = FortranFile('test.dat','r')
a = f.read_record(endian+'('+str(ix)+','+str(jx)+','+str(kx)+')d')
f.close()
```
`np.dtype`を用いて以下のようにしても読み込むことが出来る

```python
import numpy as np
from scipy.io import FortranFile

ix, jx, kx = 3, 4, 5
endian='<'
a = np.ones((ix,jx,kx))
b = np.ones((ix,jx,kx))*2

dtype = np.dtype([('a',endian+str(ix*jx*kx)+'d')])
with FortranFile('test.dat','r') as f:
    a = f.read_record(dtype=dtype)
a = a['a'].reshape(ix,jx,kx,order='F')
```