python - Cython中类似于Fortran的数组切片

Question

我正在寻找一种简洁的和高效的方法来获取数组的多维 slice ，对这些 slice 执行标量和矩阵算术，然后最终将所得数组另存为另一个数组中的 slice 。

您可以使用以下语法在fortran中做到这一点:

real*8, dimension(4,4,4,4)               ::  matrix_a
real*8, dimension(4,4)               ::  matrix_b
...
matrix_a(:, 2, :, 4) =  matrix_a(:, 2, :, 4) + (2 * matrix_b(:, :))

我正在尝试在Cython中找到实现此目的的方法。这是我能想到的最好的方法:

cdef double matrix_a[4][4][4][4]
cdef double matrix_b[4][4]
...
cdef int i0, i1
for i0 in range(4):
  for i1 in range(4):
    matrix_a[i0][1][i1][3] += (2 * matrix_b[i0][i1])

显然，您可以使用Numpy数组以非常简洁的方式执行此操作，但这似乎比上面的代码慢了几个数量级。有什么办法可以让我两全其美吗？可能我可以以更快的方式调用Numpy的某种方式？或者也许可以从Cython中利用某种C / C++ API。谢谢

Answer 1

好问题。最简洁的答案是不。

长答案...取决于您想要什么以及您愿意付出多少努力。没有一个答案，因为问题很广泛，而不是因为不可能进行操作。

在撰写本文时，尚无默认方法可对cython的内存 View 执行fortran速度算术。
Cython邮件列表上对此进行了多次讨论。
这个recent thread很好地概括了这种情况。

现在，有很多方法可以使用NumPy数组进行快速运算。没有任何一个具有其他所有功能，但是有些可以满足您的需求。

第一:

NumPy的 vector 化写得很好。类型推断有一些前期开销，但是对于紧密循环中的小数组操作而言，大多数计算时间都是在C语言中经过非常优化的循环中运行的。如果您要进行简单的in-即使在Cython中，也要进行这样的操作，最好的方法是执行np.add(a, b, out=c)这样的操作。如果您有更具体的需求，np.einsum，scipy的blas / lapack包装器和numpy的数组减少方法可提供更大的灵活性。

在大多数情况下，如果 a和 b是numpy数组，则类似

a[:,2,:,4] += 2 * b

应该足够了。

即使您在Cython中并且在内存 View 而不是numpy数组上进行操作，也可以执行以下操作:

import numpy as np
np.add(a, b, out=c)

调用这样的Python函数有固定的成本，但是如果数组很大，那就没关系了。

numpy在您的情况下仍然会变慢有两个原因。
首先，由于numpy在运行程序中处理了很多类型和形状元数据，因此启动循环操作的固定成本很高。
对于较大的阵列，此差异将消失，因此，除非您处于紧密循环内，否则无关紧要。
其次，numpy和fortran的默认内存布局不同。
这种差异可能只会在较大数组的操作中显示出来。
如果按fortran内存顺序初始化numpy数组，并且数组很大，则numpy应该与fortran保持一致。

现在，对于确实处于numpy速度不够快或无法很好地控制数组的情况的人，还有很多其他选择。
以下是一些更有前途的产品:

在Cython的cdef函数内部编写带循环的算术运算。让它接受内存 View 作为参数，并传递预分配的输出数组。
签名可能看起来像cdef inline void add(double[:,:] a, double[:,:] b, double[:,:] out) nogil:可以对Memoryviews进行 slice ，而无需Python调用，因此您可以将 slice 传递给已定义的函数，而不会产生Python调用的开销。
这种方法不是最快的方法，但是老实说，它通常已经足够好了。
它还避免了大多数麻烦，因为它们必须手动管理阵列的内存布局。
这也比使用原始C数组做相同的事情容易，因为memoryviews的内存是动态管理的。

存在几个支持numpy的自动 vector 处理程序。
考虑尝试numba，parakeet或pythran。
这些软件包主要集中于编译对numpy函数进行操作的特定于类型的函数版本，然后在适用时调用这些例程。
我与numba和长尾小鹦鹉的成绩不一，与pythran的比赛也很少。
据我所知，numba和长尾小鹦鹉目前不进行任何形式的静态表达分析，尽管它不一定超出其范围。
Pythran看起来很有希望。
它的主要卖点是对数组表达式(例如Fortran)执行静态分析。
这些库易于使用，因此值得尝试一下它们是否有用。
他们特别擅长优化循环。

您还可以使用一些Python级别的表达式优化器。
参见numexpr和Theano。
概括地说，这些软件包集中于将算术表达式编译为机器代码，而不是整个函数。
据我所知，它们并没有针对 slice 进行优化，但是您可以将其传递给numpy数组的 slice 。
这些是最优化大型数组上操作的最佳方法，但是它们通常可以为表达式提供代码，比使用numpy数组进行直接算术更快。

您可以在fortran中实现该操作，并使用Cython对其进行包装。在fortran90 website上描述了如何使用fortran进行此操作。不久前，我在this answer中编写了一个有关如何执行此操作的示例。
您必须自己管理内存布局，但这并不难，而且会很快。

您还可以使用多种C++库中的任何一种进行数组操作。
不幸的是，它们中的许多(eigen，armadillo，blaze-lib)目前尚不支持高维数组，但是您拥有的特定表达式实际上只是跨距二维数组的算术运算，因此它可以正常工作。
如果您对此约束没问题，那么可以做出一些最初的努力来为犰狳(Cyton绑定(bind))(cy-arma进行开发，尽管它们仍未开发。
确实支持n维数组操作的一个c ++库是libdynd。
它还具有一组用Cython编写的更好开发的python bindings。
在C++中实现操作并通过Cython将其包装为Python可能是最简单的方法，但是您可以从包装中取出C++ pxd声明，然后尝试在Cython级别上编写函数。

有一个名为ceygen的库，该库使用eigen作为后端在内存 View 上以函数调用的形式编写算术运算。
我很难编译它，但是从概念上讲，这可能就是您想要的。
我不确定它是否可以很好地处理大步信息。

您可以用C自己实现这些操作。只有当您想将sse内在函数应用于数组时(或某种比纯循环更好的优化方法)，这种情况才会更好。编写这样的事情非常耗时，但是您可以做到。如果只是循环，最好在Cython中编写循环并在memoryviews上进行操作。

您可以将NumPy C API用于ufuncs或gufuncs。这使您可以执行标量(或数组)操作并生成将其应用于更通用数组的函数。这些文档目前还不是很完整。

这不适用于Cython中的基本算术运算，但是不久后将添加到scipy的新功能之一是用于BLAS和LAPACK的Cython级别包装器。这样一来，您可以进行各种线性代数运算，而无需为每个函数调用花费很少的成本。现在，这是相关的pull request。

其中一些要求您自己管理内存布局，但这并不难。

如果我穿上鞋子，我会尝试选项1，然后尝试选项4。选项1将相当快，但是如果这还不够，选项4可能会更快，尽管现在您必须担心内存问题数组的布局。