NumPy：高效读取大型数组

Question

我有一个包含密集的n*m个32位浮点数的二进制文件。最高效的方法是什么，可以将其读入Fortran-排序的numpy数组中？

该文件大小为多GB。我可以控制格式，但必须紧凑（即长度约为4*n*m字节），并且必须易于从非Python代码生成。

编辑：由于数据量很大，所以必须直接生成Fortran-排序矩阵（因此我无法承担创建C-排序矩阵，然后将其转换为单独的Fortran-排序副本的成本）。

- NPE

---

http://www.scipy.org/Cookbook/InputOutput 这个网站是否回答了你的问题？（请参考“二进制文件”部分） - nimrodm

@nimrodm 谢谢。事实上，我已经尝试了其中一些方法。我提出这个问题是希望有人能够前来提供第一手经验或者熟悉 numpy 内部机制并从那个角度给予建议。 - NPE

通常，当我读取非常大的数组到numpy中时，我发现需要事先知道大小，以便预先分配适当的数组来保存数据。你事先知道大小吗？如果不知道，可以尝试使用两次扫描的方法：首先进行扫描以发现数据的大小/维度，然后分配数组，最后将数据读取/解析到该数组中。 - Peter Hansen

@Peter 很好的观点，谢谢。我事先知道大小（我控制数据格式，因此可以将大小写入文件头的一部分）。 - NPE

Answer 1

NumPy提供fromfile()函数来读取二进制数据。

a = numpy.fromfile("filename", dtype=numpy.float32)

将创建一个包含您的数据的一维数组。要将其作为二维Fortran排序的n x m矩阵访问，您可以重塑它：

a = a.reshape((n, m), order="FORTRAN")

[编辑：在这种情况下，reshape()实际上会复制数据（请参见评论）。要避免复制，请使用

a = a.reshape((m, n)).T

感谢Joe Kingtion指出这一点。

但是，说实话，如果您的矩阵有几个千兆字节，我会选择像h5py或PyTables这样的HDF5工具。这两个工具都有FAQ条目将该工具与另一个进行比较。我通常更喜欢h5py，尽管PyTables似乎更常用（两个项目的范围略有不同）。

HDF5文件可以从用于数据分析的大多数编程语言中编写。链接的维基百科文章中的接口列表不完整，例如还有一个R interface。但我实际上不知道您想使用哪种语言来编写数据...

- Sven Marnach

---

@Sven，您能否在您的示例中澄清一些事情。如果我在一个4GB的文件上执行a = numpy.fromfile("filename", dtype=numpy.float32)，然后执行a = a.reshape((n, m), order="FORTRAN")，这样做是否有可能在内存中创建一个4GB的“C”矩阵，然后立即再创建一个4GB的内存副本将其翻转为Fortran格式？ - NPE

@Joe 即使在 C 和 Fortran 格式之间进行翻转（按行或按列存储），这怎么可能呢？ - NPE

@Joe 没问题。我知道改变numpy矩阵的形状不会复制数据。然而，出于性能原因，我关心内存中的存储方式，并且我不能承受两个相同的多千兆字节矩阵的副本，因此我有所有这些问题。 - NPE

@aix - 我假设您正在从磁盘中读取Fortran有序数组。 您会得到一个实际上是ixj Fortran有序数组的平坦数组。 这与jxi C-ordered数组相同，只是它被转置了。 因此，我们将其重塑为jxi，然后将其转置为ixj。 Numpy认为它是C-ordered jxi被视为ixj，但这相当于直接存储为Fortran有序ixj在内存中。 - Joe Kington

@aix - 它在内存中是Fortran顺序...唯一的区别是numpy的标志...如果您要从Fortran访问numpy数组的内存缓冲区，您将直接读取它作为ixj数组。np.fromfile将按照磁盘上存储的顺序读取它。就像我说的，我假设它已经按Fortran顺序存储在磁盘上。如果您想将C顺序数组写入磁盘作为Fortran顺序数组，只需使用a.ravel('F').tofile(fid)。 - Joe Kington

Answer 2

基本上，Numpy将数组存储为平坦向量。多个维度只是由Numpy迭代器使用的不同视图和步幅创建的幻觉。

要了解Numpy内部工作原理的详细但易于理解的解释，请参见优秀的《The Beatiful Code》书籍第19章。

至少Numpy array()和reshape()有一个参数用于C（'C'）、Fortran（'F'）或保留顺序（'A'）。 还请参阅问题如何强制numpy数组顺序为Fortran风格？

默认C索引（行主序）的示例：

>>> a = np.arange(12).reshape(3,4) # <- C order by default
>>> a
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11]])
>>> a[1]
array([4, 5, 6, 7])

>>> a.strides
(32, 8)

使用Fortran顺序（列主序）进行索引：

>>> a = np.arange(12).reshape(3,4, order='F')
>>> a
array([[ 0,  3,  6,  9],
       [ 1,  4,  7, 10],
       [ 2,  5,  8, 11]])
>>> a[1]
array([ 1,  4,  7, 10])

>>> a.strides
(8, 24)

---

另一种视角

此外，您可以始终使用数组的参数T来获取另一种视角：

>>> a = np.arange(12).reshape(3,4, order='C')
>>> a.T
array([[ 0,  4,  8],
       [ 1,  5,  9],
       [ 2,  6, 10],
       [ 3,  7, 11]])

>>> a = np.arange(12).reshape(3,4, order='F')
>>> a.T
array([[ 0,  1,  2],
       [ 3,  4,  5],
       [ 6,  7,  8],
       [ 9, 10, 11]])

---

您也可以手动设置步幅：

>>> a = np.arange(12).reshape(3,4, order='C')
>>> a
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11]])
>>> a.strides
(32, 8)
>>> a.strides = (8, 24)
>>> a
array([[ 0,  3,  6,  9],
       [ 1,  4,  7, 10],
       [ 2,  5,  8, 11]])

- peterhil