从.npy文件读取数组到Fortran 90中

Question

我正在使用Python生成一些初始数据，以2D数组的形式呈现，例如“X”，然后使用Fortran对它们进行一些计算。最初，当数组大小约为10,000 x 10,000时，np.savetxt在速度方面表现良好。但是，一旦我开始增加数组的维数，savetxt的速度就会显着变慢。因此，我尝试了np.save，并且这样保存速度更快，但文件以.npy格式保存。如何在Fortran中读取这样的文件以重构原始数组？从我所了解的情况来看，二进制通常可以实现最低的空间消耗和最快的速度。

在Fortran 90中，

open(10,file='/home/X.npy')

read(10,*)X

我收到了以下错误信息：Fortran运行时错误：列表输入项1中的实数有误。

编辑： 在Python中，我这样做，

import numpy as np 
x = np.linspace(-50,50,10)
y = np.linspace(-50,50,10) 
X,Y = np.meshgrid(x,y) 
np.save('X',X) 

然后在Fortran中，我这样做：

program read_unformatted_binary_from_python
    use iso_c_binding
    implicit none

    integer, parameter :: DPR = selected_real_kind(p=15)
    character(len=*), parameter :: filename = 'X.npy'

    integer, parameter :: N = 10
    real(DPR), allocatable, dimension(:, :) :: X

    allocate(X(N, N))


    open(40, file=filename, status='old', access='stream', form='unformatted')
    read(40) X
    close(40)

    write(*,*) X

end program read_unformatted_binary_from_python

输出以1.8758506894003703E-309，1.1711999948023422E+171，5.2274167985502976E-037，8.4474009688929314E+252，2.6514123210345660E+180，9.9215260506210473E+247，2.1620996769994603E+233，7.5805790251297605E-096，3.4756671925988047E-152，6.5549091408576423E-260，-50.000000000000000，-38.888888888888886，-27.777777777777779等内容开头。只有使用.bin格式时才不会出现此错误，而使用np.save则会导致npy格式。

- ThunderFlash

---

Fortran不会自动读取和正确解释任意二进制文件。您必须找到或编写一个特定于该格式的例程。它相对来说有很好的文档，您喜欢的搜索引擎将帮助您追踪其描述。顺便说一下，错误消息表明文件的前4个（或可能是8个）字节无法解释为Fortran“real”。 - High Performance Mark

更糟糕的是，在代码片段中，它已经被打开为格式化读取，因此它试图将数字作为文本读取。这是行不通的。你在谈论二进制，但你没有为二进制（未格式化）I/O打开文件。 - Vladimir F Героям слава

请一定使用标签[tag:fortran]，以便更多的人能看到您的问题。Fortran 90只是一个古老而已过时的版本。对于这个任务来说，Fortran 90是不足够的，您需要从Fortran 2003中使用流I/O。 - Vladimir F Героям слава

对，我没注意到这一点……我也读过 HDF5 在处理大数据量时非常出色，尤其是在高性能计算环境下。你认为它是二进制的一个好替代品吗？ - ThunderFlash

Answer 1

弗拉基米尔F是正确的，您需要在Fortran中获得“原始二进制”文件的“流”访问。这里是一个最小化可运行代码示例：

Python

import numpy as np
A = np.random.rand(10000, 10000)
print(A.sum())
A.tofile('data.bin')

Fortran

program read_unformatted_binary_from_python
    use iso_c_binding
    implicit none

    integer, parameter :: DPR = selected_real_kind(p=15)
    character(len=*), parameter :: filename = 'data.bin'

    integer, parameter :: N = 10000
    real(DPR), allocatable, dimension(:, :) :: dat

    allocate(dat(N, N))

    open(40, file=filename, status='old', access='stream', form='unformatted')
    read(40) dat
    close(40)

    write(*,*) sum(dat)

end program read_unformatted_binary_from_python

我的Fortran示例可能比必要的还要长，因为我使用了许多不同的系统和编译套件，并且也不喜欢大的静态数组（毕竟我是Fortran用户）。

我使用Python 2.7.x、Numpy 13.x和Homebrew GCC 6.3.0_1上的gfortran在MacBook Pro上快速编码，但这应该在所有系统上都能正常工作。

更新： 这里需要特别注意数组的形状和大小。如果dat被分配的比文件中的要大，则流式read应该尝试填充整个数组，遇到EOF符号并发出错误。在Python中，np.fromfile()方法将读取直到EOF然后返回一个具有适当长度的一维数组，即使A最初是多维的。这是因为原始二进制没有元数据，只是来自RAM的连续字节字符串。

因此，以下Python代码会产生相同的文件：

A = np.random.rand(10000, 10000)
A.tofile('file.whatever')
A.ravel().tofile('file.whatever')
A.reshape((100, 1000, 1000)).tofile('file.whatever')

那个文件可以被读取并重新塑形为:

B = np.fromfile('file.whatever').reshape(A.shape)
B = np.fromfile('file.whatever').reshape((100, 1000, 100, 10))
# or something like
B = np.fromfile('file.whatever') # just a 1D array
B.resize(A.shape)  # resized in-place

在Fortran中，使用流式访问读取整个未知大小的原始文件非常容易，但显然您需要某种用户输入来重新塑造数据：

program read_unformatted_binary_from_python
    use iso_c_binding
    implicit none

    integer, parameter :: DPR = selected_real_kind(p=15)
    character(len=*), parameter :: filename = 'data.bin'
    integer :: N = 10000, bytes, reals, M
    real(DPR), allocatable :: A(:,:), D(:, :), zeros(:)
    real(DPR), allocatable, target :: B(:)
    real(DPR), pointer :: C(:, :)

    allocate(A(N, N))

    open(40, file=filename, status='old', access='stream', form='unformatted')

    read(40) A
    write(*,*) 'sum of A', sum(A)

    inquire(unit=40, size=bytes)
    reals = bytes/8
    allocate(B(reals))

    read(40, pos=1) B
    write(*,*) 'sum of B', sum(B)

    ! now reshape B in-place assuming the user wants the first dimension 
    ! (which would be the outer dimension in Python) to be length 100
    N = 100
    if(mod(reals, N) == 0) then
         M = reals/N
         call C_F_POINTER (C_LOC(B), C, [N, M])
         write(*, *) 'sum of C', sum(C)
         write(*, *) 'shape of C', shape(C)
    else
         write(*,*) 'file size is not divisible by N!, did not point C to B'
    end if

    ! now reshape B out-of-place with first dimension length 9900, and
    ! pad the result so that there is no size mismatch error  
    N = 9900
    M = reals/N
    if(mod(reals, N) > 0) M=M+1

    allocate(D(N, M))
    allocate(zeros(N), source=real(0.0, DPR))
    D = reshape(B, [N, M], pad=zeros)

    write(*,*) 'sum of D', sum(D)
    write(*,*) 'shape of D', shape(D)

    ! obviously you can also work with shape(A) in fortran the same way you
    ! would use A.shape() in Python, if you already knew that A was the
    ! correct shape of the data
    deallocate(D)
    allocate(D, mold=A)
    D = reshape(B, shape(A))
    write(*,*) 'sum of D', sum(D)
    write(*,*) 'shape of D', shape(D)

    ! or, just directly read it in, skipping the whole reshape B part
    read(40, pos=1) D
    write(*,*) 'sum of D', sum(D)

    close(40)

end program read_unformatted_binary_from_python

- CAZT

---

显然，还有一个numpy库可用于Fortran中读取.npy文件。请参阅此Cookbook底部，其中包括libnpy的下载链接： http://scipy-cookbook.readthedocs.io/items/InputOutput.html - CAZT

太棒了！我还没有在我的电脑上尝试过这个，但这是向后兼容的吗？也就是说，“stream”会在F 90中起作用，还是仅在2003及以上版本中可用？libnpy似乎是一个非常好的工具，我会去查看一下。 - ThunderFlash

流是Fortran 2003。我已经说过了，所以它在Fortran 90中不起作用。Fortran 90已经过时了，没有人使用它。有些人使用Fortran 95，但通常还需要使用许多非标准扩展。Fortran 90已经死了。除了一些非常旧的版本外，没有Fortran 90存在，所以不用担心。 - Vladimir F Героям слава

我明白了...另外，当我使用一个维度为10 x 10的数组，其中数字范围从-50到50运行上述代码时，当我写出该数组时，它不仅显示了-50到50，还在标题部分显示了非常大的数字。我该如何消除它们，以便它不会影响进一步的计算？ - ThunderFlash

请看我上面的更新，简而言之：原始二进制文件没有元数据，因此fromfile()只会将数据读入一个连续的一维数组中，但您始终可以使用reshape()或resize()来调整数组。 - CAZT