Dav*_*ler 3 python arrays fortran cython fortran90
即时寻找一个简明和高效的方式,采取多维片阵列,就这些切片执行标量和矩阵算术,然后最终得到的阵列保存为另一个阵列的切片.
您可以使用以下语法在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.谢谢
好问题.最简洁的答案是不.
答案很长......这取决于你想要什么以及你愿意付出多少努力.没有一个答案,因为问题非常广泛,不是因为操作不可行.
在撰写本文时,没有默认方法可以对cython的内存视图执行fortran-speed算法.在Cython邮件列表上有关于此的多个讨论.这个最近的线程很好地总结了这种情况.
现在,有很多方法可以使用NumPy数组进行快速算术运算.没有一个具有其他功能,但有些可能适合您的需求.
第一:
np.add(a, b, out=c).如果您有更多特定需求,np.einsumscipy的blas/lapack包装器和numpy的数组减少方法提供了更大的灵活性.在大多数情况下,如果a和b是numpy数组,像
a[:,2,:,4] += 2 * b
应该足够了.
即使你在Cython中运行内存视图而不是numpy数组,你也可以这样做:
import numpy as np
np.add(a, b, out=c)
调用像这样的Python函数有一个固定的成本,但是如果数组很大,那就不重要了.
在你的情况下,numpy仍然会慢一两个原因.首先,启动循环操作的固定成本很高,因为numpy会在runtine处理很多类型和形状元数据.对于较大的数组,这种差异消失了,所以,除非你在一个紧密的循环中,否则无关紧要.其次,numpy和fortran的默认内存布局是不同的.这种差异可能只会出现在具有较大阵列的操作中.如果以fortran内存顺序初始化numpy数组并且数组很大,那么numpy应该与fortran相提并论.
现在,对于那些真正处于numpy不够快或者对阵列没有足够精确控制的情况下,还有很多其他选择.以下是一些更有希望的:
cdef在Cython中用函数内部循环编写算术运算.让它接受内存视图作为参数并传入预先分配的输出数组.签名可能看起来像cdef inline void add(double[:,:] a, double[:,:] b, double[:,:] out) nogil:
Memoryviews可以切片而没有Python调用的开销,因此您可以将切片传递给您定义的函数,而不会产生Python调用的开销.这种方法不是最快的,但老实说,它通常都足够好.它还避免了必须手动管理阵列的内存布局所带来的大部分痛苦.这比使用原始C数组做同样的事情更容易,因为内存视图会动态管理它们的内存.
存在几个numpy-aware auto-vecotrizers.考虑尝试numba,长尾小鹦鹉或pythran.这些包主要集中在编译特定于类型的函数版本,这些函数在numpy函数上运行,然后在适用的情况下调用这些例程.我和numba和长尾小鹦鹉的结果好坏参半,并且没有和pythran玩过多.据我所知,numba和长尾小鹦鹉目前没有进行任何形式的静态表达分析,尽管它不一定超出其范围.Pythran看起来很有希望.它的主要卖点是对数组表达式进行静态分析(如Fortran).这些库易于使用,因此可能值得尝试它们以确定它们是否有用.他们特别擅长优化循环.
您还可以使用一些Python级别的表达式优化器.请参见numexpr和Theano.粗略地说,这些包专注于将算术表达式编译为机器代码,而不是整个函数.据我所知,他们没有针对切片进行优化,但是你可以将它们传递给numpy数组.这些最适合优化大型数组上的操作,但它们通常可以为表达式提供比使用numpy数组的直接算法更快的代码.
您可以在fortran中实现该操作并使用Cython将其包装.在fortran90网站上有关于如何使用fortran执行此操作的说明.我在这个答案中写了一个如何做到这一点的例子.你必须自己管理内存布局,但这不是太难,而且会非常快.
您还可以使用各种C++库中的任何一个进行数组操作.不幸的是,其中许多(eigen,armadillo,blaze-lib)目前没有内置支持高维数组,但你所拥有的特定表达式实际上只是一个跨步的二维数组算术,所以它可以正常工作.如果你对这个约束没有问题,那么最初的努力就是为犰狳制作Cython绑定(cy-arma,虽然它们仍然不发达.一个支持n维数组运算的c ++库是libdynd.它也有一个更好的用Cython编写的一组python绑定.用C++实现操作并用Cython将它们包装成Python可能最简单,但是你可以从包装器中取出C++ pxd声明并尝试编写你的函数在Cython级别.
有一个名为ceygen的库,它使用eigen作为后端,将内存视图上的算术运算写为函数调用.我在编译时遇到了麻烦,但从概念上讲,这可能就是你要找的东西.我不确定它是否能很好地处理跨步信息.
你可以在C中自己实现这些操作.唯一一次,当你想将sse内在函数应用于数组时(或者某种优于纯循环的优化).写这样的东西非常耗时,但你可以做到.如果你只是循环,你最好在Cython中编写循环并在内存视图上操作.
您可以将NumPy C API用于ufuncs或gufuncs.这允许您进行标量(或数组)操作并生成将其应用于更通用数组的函数.这些文档现在还不是很完整.
这对于Cython中的基本算术运算不起作用,但是很快将添加到scipy的一个新功能是BLAS和LAPACK的Cython级包装器.这将使您可以进行各种线性代数运算,而且每个函数调用的成本都很低.现在,这是相关的拉取请求.
其中一些将要求您自己管理内存布局,但这并不是非常困难.
如果我在你的鞋子,我会尝试选项1,然后选项4.选项1将相当快,但如果这还不够,选项4可能会更快,但现在你将不得不担心内存你自己的阵列布局.