python - C数组与NumPy数组

Question

在性能（代数运算，查找，缓存等）方面，C数组（可以显示为C数组或cython.view.array [Cython数组]或上述两个的内存视图）之间有区别）和NumPy数组（在Cython中应该没有Python开销）

编辑：

我应该提到的是，在NumPy数组中使用Cython进行了静态类型化，而dtype是NumPy编译时datype（例如，cdef np.int_t或cdef np.float32_t），在C情况下的类型是C等效项（< cc>和cdef int_t）

编辑2：

以下是Cython Memoryview documentation中的示例，以进一步说明我的问题：

from cython.view cimport array as cvarray
import numpy as np

# Memoryview on a NumPy array
narr = np.arange(27, dtype=np.dtype("i")).reshape((3, 3, 3))
cdef int [:, :, :] narr_view = narr

# Memoryview on a C array
cdef int carr[3][3][3]
cdef int [:, :, :] carr_view = carr

# Memoryview on a Cython array
cyarr = cvarray(shape=(3, 3, 3), itemsize=sizeof(int), format="i")
cdef int [:, :, :] cyarr_view = cyarr

坚持 cdef float vs C array vs Cython array有什么区别吗？

Answer 1

我对此的知识仍不完善，但这可能会有所帮助。
我运行了一些非正式的基准测试，以表明每种数组类型都有哪些用处，并对我发现的内容很感兴趣。

尽管这些数组类型在许多方面有所不同，但是如果您要对大型数组进行繁重的计算，则您应该能够从其中的任何一种中获得相似的性能，因为逐项访问在整体上应该大致相同。

NumPy数组是使用Python的C API实现的Python对象。
NumPy数组确实提供C级的API，但是不能独立于Python解释器创建它们。
由于NumPy和SciPy中提供了所有不同的数组操作例程，因此它们特别有用。

Cython内存视图也是Python对象，但是它是Cython扩展类型。
它似乎不是为纯Python设计的，因为它不是可以直接从Python导入的Cython的一部分，但是您可以从Cython函数将视图返回给Python。
您可以在https://github.com/cython/cython/blob/master/Cython/Utility/MemoryView.pyx查看实现

C数组是C语言中的本机类型。
它像指针一样被索引，但是数组和指针是不同的。
在http://c-faq.com/aryptr/index.html上对此进行了很好的讨论
它们可以在堆栈上分配，并且C编译器更容易优化它们，但是在Cython之外访问它们将更加困难。
我知道您可以从由其他程序动态分配的内存中创建一个NumPy数组，但是用这种方法似乎要困难得多。
Travis Oliphant在http://blog.enthought.com/python/numpy-arrays-with-pre-allocated-memory/上发布了此示例
如果您使用C数组或指针在程序中进行临时存储，则它们应该对您来说很好。
它们对于切片或任何其他类型的矢量化计算将不那么方便，因为您必须自己使用显式循环来做所有事情，但是它们应更快地分配和释放，并应为速度提供良好的基准。

Cython还提供了一个数组类。
看起来它是为内部使用而设计的。
实例是在复制memoryview时创建的。
见http://docs.cython.org/src/userguide/memoryviews.html#view-cython-arrays

在Cython中，您还可以分配内存并索引一个指针，以将分配的内存视为数组。
见http://docs.cython.org/src/tutorial/memory_allocation.html

这里有一些基准测试显示了索引大型数组的性能。
这是Cython文件。

from numpy cimport ndarray as ar, uint64_t
cimport cython
import numpy as np

@cython.boundscheck(False)
@cython.wraparound(False)
def ndarr_time(uint64_t n=1000000, uint64_t size=10000):
    cdef:
        ar[uint64_t] A = np.empty(n, dtype=np.uint64)
        uint64_t i, j
    for i in range(n):
        for j in range(size):
            A[j] = n

def carr_time(uint64_t n=1000000):
    cdef:
        ar[uint64_t] A = np.empty(n, dtype=np.uint64)
        uint64_t AC[10000]
        uint64_t a
        int i, j
    for i in range(n):
        for j in range(10000):
            AC[j] = n

@cython.boundscheck(False)
@cython.wraparound(False)
def ptr_time(uint64_t n=1000000, uint64_t size=10000):
    cdef:
        ar[uint64_t] A = np.empty(n, dtype=np.uint64)
        uint64_t* AP = &A[0]
        uint64_t a
        int i, j
    for i in range(n):
        for j in range(size):
            AP[j] = n

@cython.boundscheck(False)
@cython.wraparound(False)
def view_time(uint64_t n=1000000, uint64_t size=10000):
    cdef:
        ar[uint64_t] A = np.empty(n, dtype=np.uint64)
        uint64_t[:] AV = A
        uint64_t i, j
    for i in range(n):
        for j in range(size):
            AV[j] = n

使用IPython计时这些时间，我们获得

%timeit -n 10 ndarr_time()
%timeit -n 10 carr_time()
%timeit -n 10 ptr_time()
%timeit -n 10 view_time()

10 loops, best of 3: 6.33 s per loop
10 loops, best of 3: 3.12 s per loop
10 loops, best of 3: 6.26 s per loop
10 loops, best of 3: 3.74 s per loop

考虑到根据 Efficiency: arrays vs pointers，数组不太可能比指针快得多，所以这些结果使我感到有些奇怪。
似乎某种编译器优化正在使纯C数组和类型化的内存视图更快。
我试图关闭我的C编译器上的所有优化标志并获得计时

1 loops, best of 3: 25.1 s per loop
1 loops, best of 3: 25.5 s per loop
1 loops, best of 3: 32 s per loop
1 loops, best of 3: 28.4 s per loop

在我看来，逐项访问几乎完全相同，不同之处在于C数组和Cython内存视图似乎更易于编译器优化。

我在一段时间前发现的这两篇博客文章中可以找到对此的更多评论：
http://jakevdp.github.io/blog/2012/08/08/memoryview-benchmarks/
http://jakevdp.github.io/blog/2012/08/16/memoryview-benchmarks-2/

在第二篇博客文章中，他评论了如何内联内存视图切片，它们如何提供类似于指针算法的速度。
我在自己的一些测试中注意到，并不总是显式地内联使用Memory View slice的函数。
作为示例，我将计算数组两行的每种组合的内积。

from numpy cimport ndarray as ar
cimport cython
from numpy import empty

# An inlined dot product
@cython.boundscheck(False)
@cython.wraparound(False)
cdef inline double dot_product(double[:] a, double[:] b, int size):
    cdef int i
    cdef double tot = 0.
    for i in range(size):
        tot += a[i] * b[i]
    return tot

# non-inlined dot-product
@cython.boundscheck(False)
@cython.wraparound(False)
cdef double dot_product_no_inline(double[:] a, double[:] b, int size):
    cdef int i
    cdef double tot = 0.
    for i in range(size):
        tot += a[i] * b[i]
    return tot

# function calling inlined dot product
@cython.boundscheck(False)
@cython.wraparound(False)
def dot_rows_slicing(ar[double,ndim=2] A):
    cdef:
        double[:,:] Aview = A
        ar[double,ndim=2] res = empty((A.shape[0], A.shape[0]))
        int i, j
    for i in range(A.shape[0]):
        for j in range(A.shape[0]):
            res[i,j] = dot_product(Aview[i], Aview[j], A.shape[1])
    return res

# function calling non-inlined version
@cython.boundscheck(False)
@cython.wraparound(False)
def dot_rows_slicing_no_inline(ar[double,ndim=2] A):
    cdef:
        double[:,:] Aview = A
        ar[double,ndim=2] res = empty((A.shape[0], A.shape[0]))
        int i, j
    for i in range(A.shape[0]):
        for j in range(A.shape[0]):
            res[i,j] = dot_product_no_inline(Aview[i], Aview[j], A.shape[1])
    return res

# inlined dot product using numpy arrays
@cython.boundscheck(False)
@cython.boundscheck(False)
cdef inline double ndarr_dot_product(ar[double] a, ar[double] b):
    cdef int i
    cdef double tot = 0.
    for i in range(a.size):
        tot += a[i] * b[i]
    return tot

# non-inlined dot product using numpy arrays
@cython.boundscheck(False)
@cython.boundscheck(False)
cdef double ndarr_dot_product_no_inline(ar[double] a, ar[double] b):
    cdef int i
    cdef double tot = 0.
    for i in range(a.size):
        tot += a[i] * b[i]
    return tot

# function calling inlined numpy array dot product
@cython.boundscheck(False)
@cython.wraparound(False)
def ndarr_dot_rows_slicing(ar[double,ndim=2] A):
    cdef:
        ar[double,ndim=2] res = empty((A.shape[0], A.shape[0]))
        int i, j
    for i in range(A.shape[0]):
        for j in range(A.shape[0]):
            res[i,j] = ndarr_dot_product(A[i], A[j])
    return res

# function calling nun-inlined version for numpy arrays
@cython.boundscheck(False)
@cython.wraparound(False)
def ndarr_dot_rows_slicing_no_inline(ar[double,ndim=2] A):
    cdef:
        ar[double,ndim=2] res = empty((A.shape[0], A.shape[0]))
        int i, j
    for i in range(A.shape[0]):
        for j in range(A.shape[0]):
            res[i,j] = ndarr_dot_product(A[i], A[j])
    return res

# Version with explicit looping and item-by-item access.
@cython.boundscheck(False)
@cython.wraparound(False)
def dot_rows_loops(ar[double,ndim=2] A):
    cdef:
        ar[double,ndim=2] res = empty((A.shape[0], A.shape[0]))
        int i, j, k
        double tot
    for i in range(A.shape[0]):
        for j in range(A.shape[0]):
            tot = 0.
            for k in range(A.shape[1]):
                tot += A[i,k] * A[j,k]
            res[i,j] = tot
    return res

计时这些，我们看到了

A = rand(1000, 1000)
%timeit dot_rows_slicing(A)
%timeit dot_rows_slicing_no_inline(A)
%timeit ndarr_dot_rows_slicing(A)
%timeit ndarr_dot_rows_slicing_no_inline(A)
%timeit dot_rows_loops(A)

1 loops, best of 3: 1.02 s per loop
1 loops, best of 3: 1.02 s per loop
1 loops, best of 3: 3.65 s per loop
1 loops, best of 3: 3.66 s per loop
1 loops, best of 3: 1.04 s per loop

具有显式内联的结果与没有内联时一样快。
在这两种情况下，类型化的内存视图都可以与未切片的函数版本进行比较。

在博客文章中，他必须编写一个特定的示例来强制编译器不内联函数。
看起来像样的C编译器（我正在使用MinGW）可以解决这些优化问题，而无需告诉他们内联某些函数。
即使在没有显式内联的情况下，在Cython模块内的函数之间传递数组切片时，Memoryview也会更快。

但是，在这种特殊情况下，即使将循环推至C并不能真正达到通过适当使用矩阵乘法所能达到的速度。
BLAS仍然是做这种事情的最好方法。

%timeit A.dot(A.T)
10 loops, best of 3: 25.7 ms per loop

从NumPy数组自动转换为memoryview，如

cimport cython

@cython.boundscheck(False)
@cython.wraparound(False)
def cysum(double[:] A):
    cdef tot = 0.
    cdef int i
    for i in range(A.size):
        tot += A[i]
    return tot

一个问题是，如果您希望函数返回NumPy数组，则必须使用 np.asarray将内存视图对象再次转换为NumPy数组。
这是相对便宜的操作，因为内存视图符合 http://www.python.org/dev/peps/pep-3118/

结论

在Cython模块内部使用类型化的内存视图似乎是NumPy数组的可行替代方案。
使用内存视图时，数组切片将更快，但为内存视图编写的函数和方法不如为NumPy数组编写的函数和方法少。
如果您不需要调用一堆NumPy数组方法并且想要简单的数组切片，则可以使用内存视图代替NumPy数组。
如果既需要切片又需要给定数组的NumPy功能，则可以创建一个内存视图，该视图指向与NumPy数组相同的内存。
然后，您可以使用视图在函数之间传递切片，并使用数组调用NumPy函数。
该方法仍然受到一定限制，但是如果您使用单个数组进行大多数处理，它将很好地工作。

C数组和/或动态分配的内存块可能对中间计算很有用，但是将它们传递回Python并在其中使用并不容易。
我认为，动态分配多维C数组也比较麻烦。
我知道的最好的方法是分配一个大的内存块，然后使用整数算术对其进行索引，就好像它是多维数组一样。
如果您想动态分配数组，这可能是个问题。
另一方面，对于C数组，分配时间可能要快很多。
其他数组类型的设计速度和便利性都差不多，因此，除非有充分的理由，否则我建议您使用它们。

更新：如@Veedrac的回答中所述，您仍然可以将Cython内存视图传递给大多数NumPy函数。
当您执行此操作时，NumPy通常将不得不创建一个新的NumPy数组对象以无论如何都与内存视图一起使用，因此这会稍微慢一些。
对于大型阵列，其影响可忽略不计。
不管数组大小如何，调用 np.asarray进行内存视图都会相对较快。
但是，为了证明这种效果，这是另一个基准：

Cython文件：

def npy_call_on_view(npy_func, double[:] A, int n):
    cdef int i
    for i in range(n):
        npy_func(A)

def npy_call_on_arr(npy_func, ar[double] A, int n):
    cdef int i
    for i in range(n):
        npy_func(A)

在IPython中：

from numpy.random import rand
A = rand(1)
%timeit npy_call_on_view(np.amin, A, 10000)
%timeit npy_call_on_arr(np.amin, A, 10000)

输出：

10 loops, best of 3: 282 ms per loop
10 loops, best of 3: 35.9 ms per loop

我试图选择一个可以很好地显示这种效果的示例。
除非涉及相对较小的数组上的许多NumPy函数调用，否则这不会花费很多时间。
请记住，无论我们以哪种方式调用NumPy，Python函数调用仍然会发生。

这仅适用于NumPy中的功能。
大多数数组方法不适用于内存视图（某些属性仍然可用，例如 size和 shape和 T）。
例如，带有NumPy数组的 A.dot(A.T)将变为 np.dot(A, A.T)。