使用 numpy 读取和操作 fortran 有序数组的最佳实践

Question

我正在阅读 ascii 和二进制文件，它们都以 fortran 顺序指定 3 维数组。我想对这些数组执行一些任意操作，然后将它们导出为相同的 ascii 或二进制格式。

我对在我的库中处理这些数组的最佳方法感到困惑。我目前的设计似乎很容易出错，因为如果创建了任何新数组，我必须不断地从默认的 C 顺序重新调整。

当前设计：

我有一些函数可以读取这些文件并返回 numpy 数组。读取函数的行为方式都类似，本质上是读取数据并返回如下内容：

return array.reshape((i, j, k), order='F')

按照我的理解，我将 fortran 顺序的视图返回到原始数组中。

我的代码假定所有数组都按 fortran 顺序排列。这意味着任何可能创建新数组的新操作我确保使用 reshape 将其转换回 fortran 顺序。

这似乎很容易出错，因为我必须密切注意任何创建新数组的操作，并确保将其重塑为 fortran 顺序，因为默认值通常是 C 顺序。

我稍后可能不得不再次将这些数组导出为二进制或 ascii，并且需要维护 fortran 排序。所以，我使用numpy.nditer 以fortran 顺序写出每个元素。

担忧：

• 目前的方法似乎很容易出错，因为我通常按 C 顺序思考。恐怕我总是会因为错过对 reshape 的调用而被咬伤，这会迫使事情按 C 顺序进行。

  • 我希望不必担心数组元素的顺序，除非在读取输入文件或将数据写入输出文件时。

• 当前的方法看起来很混乱，因为索引可以以不同的方式解释，事情可能会变得混乱。

  • 在处理 fortran 数组时，索引的元组顺序是向后的，对吗？
  • 那么，x[(1, 2, 3)] 对于 fortran 数组意味着 k = 1、j = 2 和 i = 3，而对于 C 阶数组，x[(1, 2, 3)] 意味着 k = 3、j = 2、i = 1 正确吗？李>
  • 这意味着我和我的库的用户必须始终以 (k, j, i) 顺序考虑索引，而不是我们 C/Python 程序员通常认为的 (i, j, k)。

问题：

有没有做这种事情的最佳实践？在一个理想的世界中，我想读取 fortran 有序数组，然后在导出到文件之前忘记排序。但是，恐怕我会一直误解索引等。

我已经阅读了我能找到的唯一 numpy 文档，http://docs.scipy.org/doc/numpy/reference/internals.html#multidimensional-array-indexing-order-issues。然而，这个概念对我来说仍然像泥巴一样清晰。也许我只是需要对 numpy 文档的不同解释，http://docs.scipy.org/doc/numpy/reference/internals.html#multidimensional-array-indexing-order-issues。

Answer 1

Numpy 在 Python 级别抽象出 Fortran 排序和 C 排序之间的区别。 （事实上​​，你甚至可以使用 numpy 对 >2d 数组进行其他排序。它们在 python 级别都被同等对待。）

您唯一需要担心 C 与 F 顺序的情况是当您读取/写入磁盘或将数组传递给较低级别​​的函数时。

一个简单的例子

作为一个例子，让我们以 C 顺序和 Fortran 顺序制作一个简单的 3D 数组：

In [1]: import numpy as np

In [2]: c_order = np.arange(27).reshape(3,3,3)

In [3]: f_order = c_order.copy(order='F')

In [4]: c_order
Out[4]: 
array([[[ 0,  1,  2],
        [ 3,  4,  5],
        [ 6,  7,  8]],

       [[ 9, 10, 11],
        [12, 13, 14],
        [15, 16, 17]],

       [[18, 19, 20],
        [21, 22, 23],
        [24, 25, 26]]])

In [5]: f_order
Out[5]: 
array([[[ 0,  1,  2],
        [ 3,  4,  5],
        [ 6,  7,  8]],

       [[ 9, 10, 11],
        [12, 13, 14],
        [15, 16, 17]],

       [[18, 19, 20],
        [21, 22, 23],
        [24, 25, 26]]])

请注意，它们看起来都相同（它们处于我们正在与它们交互的级别）。你怎么知道它们的顺序不同？首先，让我们看一下标志（注意C_CONTIGUOUS vs F_CONTIGUOUS）：

In [6]: c_order.flags
Out[6]: 
  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : False
  WRITEABLE : True
  ALIGNED : True
  UPDATEIFCOPY : False

In [7]: f_order.flags
Out[7]: 
  C_CONTIGUOUS : False
  F_CONTIGUOUS : True
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  UPDATEIFCOPY : False

如果您不信任这些标志，您可以通过查看arr.ravel(order='K') 有效地查看内存顺序。 order='K' 很重要。否则，当您调用arr.ravel() 时，输出将按C 顺序无论数组的内存布局如何。 order='K' 使用内存布局。

In [8]: c_order.ravel(order='K')
Out[8]: 
array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,
       17, 18, 19, 20, 21, 22, 23, 24, 25, 26])

In [9]: f_order.ravel(order='K')
Out[9]: 
array([ 0,  9, 18,  3, 12, 21,  6, 15, 24,  1, 10, 19,  4, 13, 22,  7, 16,
       25,  2, 11, 20,  5, 14, 23,  8, 17, 26])

差异实际上表示（并存储）在数组的strides 中。注意c_order 的步幅是(72, 24, 8)，而f_order 的步幅是(8, 24, 72)。

只是为了证明索引的工作方式相同：

In [10]: c_order[0,1,2]
Out[10]: 5

In [11]: f_order[0,1,2]
Out[11]: 5

阅读和写作

您会遇到此问题的主要地方是当您从磁盘读取或写入磁盘时。许多文件格式需要特定的顺序。我猜你正在使用地震数据格式，其中大多数（例如 Geoprobe .vol 的，我认为 Petrel 的卷格式也是如此）本质上是写一个二进制标题，然后是一个 Fortran 排序的 3D 数组到磁盘。

考虑到这一点，我将使用一个小的地震立方体（我的论文中的一些数据的 sn-p）作为示例。

这两个都是uint8s 的二进制数组，形状为 50x100x198。一个是 C 级的，另一个是 Fortran 级的。 c_order.datf_order.dat

阅读它们：

import numpy as np
shape = (50, 100, 198)

c_order = np.fromfile('c_order.dat', dtype=np.uint8).reshape(shape)
f_order = np.fromfile('f_order.dat', dtype=np.uint8).reshape(shape, order='F')

assert np.all(c_order == f_order)

请注意，唯一的区别是将内存布局指定为reshape。两个数组的内存布局仍然不同（reshape 不会复制），但它们在 python 级别的处理方式相同。

只是为了证明文件确实是以不同的顺序编写的：

In [1]: np.fromfile('c_order.dat', dtype=np.uint8)[:10]
Out[1]: array([132, 142, 107, 204,  37,  37, 217,  37,  82,  60], dtype=uint8)

In [2]: np.fromfile('f_order.dat', dtype=np.uint8)[:10]
Out[2]: array([132, 129, 140, 138, 110,  88, 110, 124, 142, 139], dtype=uint8)

让我们可视化结果：

def plot(data):
    fig, axes = plt.subplots(ncols=3)
    for i, ax in enumerate(axes):
        slices = [slice(None), slice(None), slice(None)]
        slices[i] = data.shape[i] // 2
        ax.imshow(data[tuple(slices)].T, cmap='gray_r')
    return fig

plot(c_order).suptitle('C-ordered array')
plot(f_order).suptitle('F-ordered array')
plt.show()

请注意，我们以相同的方式对它们进行索引，并且它们的显示方式相同。

IO 的常见错误

首先，让我们尝试读取 Fortran 排序的文件，就好像它是 C 排序的一样，然后看看结果（使用上面的 plot 函数）：

wrong_order = np.fromfile('f_order.dat', dtype=np.uint8).reshape(shape)
plot(wrong_order)

不太好！

您提到您必须使用“反向”指标。这可能是因为您通过执行以下操作修复了上图中发生的情况（注意反转的形状！）：

c_order = np.fromfile('c_order.dat', dtype=np.uint8).reshape([50,100,198])
rev_f_order = np.fromfile('f_order.dat', dtype=np.uint8).reshape([198,100,50])

让我们想象一下会发生什么：

plot(c_order).suptitle('C-ordered array')
plot(rev_f_order).suptitle('Incorrectly read Fortran-ordered array')

请注意，第一个图最右边的图像（时间片）与第二个图最左边的图像的转置版本相匹配。

同样，print rev_f_order[1,2,3] 和 print c_order[3,2,1] 都产生 140，而以相同的方式索引它们会产生不同的结果。

基本上，这就是反向索引的来源。 Numpy 认为它是一个 C 有序数组具有不同的形状。请注意，如果我们查看标志，它们在内存中都是 C 连续的：

In [24]: rev_f_order.flags
Out[24]: 
  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : False
  WRITEABLE : True
  ALIGNED : True
  UPDATEIFCOPY : False

In [25]: c_order.flags
Out[25]: 
  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : False
  WRITEABLE : True
  ALIGNED : True
  UPDATEIFCOPY : False

这是因为 fortran 有序数组等同于 C 有序数组具有相反的形状。

以 Fortran 顺序写入磁盘

以 Fortran 顺序将 numpy 数组写入磁盘时，还有一个问题。

除非您另外指定，否则数组将按 C 顺序写入，而不管其内存布局如何！ （ndarray.tofile 的文档中有一个明确的说明，但这是一个常见的问题。相反的行为是不正确的，不过，i.m.o.）

因此，无论数组的内存布局如何，要以 Fortran 顺序将其写入磁盘，您需要这样做：

arr.ravel(order='F').tofile('output.dat')

如果您将其编写为 ascii，则同样适用。使用ravel(order='F')，然后写出一维结果。