使用 C 扩展 numpy 时如何考虑列连续数组答案

【问题标题】：How to account for column-contiguous array when extending numpy with C使用 C 扩展 numpy 时如何考虑列连续数组
【发布时间】：2010-12-12 06:52:23
【问题描述】：

我有一个 C 函数来规范化日志空间中数组的行（这可以防止数值下溢）。

我的C函数原型如下：

void normalize_logspace_matrix(size_t nrow, size_t ncol, double* mat);

你可以看到它需要一个指向数组的指针并修改它。 C 代码当然假定数据保存为 C 连续数组，即行连续。

我使用 Cython 将函数包装如下（省略导入和 cdef extern from）：

def normalize_logspace(np.ndarray[np.double_t, ndim=2] mat):
    cdef Py_ssize_t n, d
    n = mat.shape[0]
    d = mat.shape[1]
    normalize_logspace_matrix(n, d, <double*> mat.data)
    return mat

大多数时候 numpy-arrays 是行连续的，并且函数工作正常。但是，如果先前已经转置了 numpy-array，则不会复制数据，而只会返回数据的新视图。在这种情况下，我的函数失败，因为数组不再是行连续的。

我可以通过将数组定义为具有 Fortran 连续顺序来解决这个问题，这样在转置后它将是 C 连续的：

A = np.array([some_func(d) for d in range(D)], order='F').T
A = normalize_logspace(A)

显然这很容易出错，用户必须注意数组的顺序是否正确，这是用户在 Python 中不需要关心的事情。

在行连续数组和列连续数组中，如何实现这一点的最佳方法是什么？我认为 Cython 中的某种数组顺序检查是可行的方法。当然，我更喜欢不需要将数据复制到新数组中的解决方案，但我几乎认为这是必要的。

【问题讨论】：

标签： python c numpy cython

【解决方案1】：

如果您想在不复制的情况下支持 C 和 Fortran 顺序的数组，您的 C 函数需要足够灵活以支持这两种顺序。这可以通过将 NumPy 数组的步幅传递给 C 函数来实现：将原型更改为

void normalize_logspace_matrix(size_t nrow, size_t ncol, 
                               size_t rowstride, size_t colstride,
                               double* mat);

Cython 调用

def normalize_logspace(np.ndarray[np.double_t, ndim=2] mat):
    cdef Py_ssize_t n, d, rowstride, colstride
    n = mat.shape[0]
    d = mat.shape[1]
    rowstride = mat.strides[0] // mat.itemsize
    colstride = mat.strides[1] // mat.itemsize
    normalize_logspace_matrix(n, d, rowstride, colstride, <double*> mat.data)
    return mat

然后，用 mat[row*rowstride + col*colstride] 替换 C 代码中出现的每个 mat[row*ncol + col]。

【讨论】：

这个 2010 年的答案是否仍然有效，还是现在有更好的方法来实现这个目标？
@larsmans：我不知道你所说的“这个”到底是什么意思。如果这是您想要的，编写一个可以同时处理 Fortran 连续和 C 连续二维数组的 C 函数仍然可以这样工作。如果您的数组可以被复制，那么还有其他解决方案（并且在 2010 年已经存在）。

【解决方案2】：

在这种情况下，您确实希望创建输入数组（可以是真实数组上的视图）的副本，并保证行连续顺序。您可以通过以下方式实现此目的：

a = numpy.array(A, copy=True, order='C')

另外，请考虑查看 Numpy 的确切 array interface（也有 C 部分）。

【讨论】：

【解决方案3】：

+1 给 Sven，他的回答解决了问题（好吧，明白了）那dstack 返回一个 F_contiguous 数组？！

# don't use dstack to stack a,a,a -> rgb for a C func

import sys
import numpy as np

h = 2
w = 4
dim = 3
exec( "\n".join( sys.argv[1:] ))  # run this.py h= ...

a = np.arange( h*w, dtype=np.uint8 ) .reshape((h,w))
rgb = np.empty( (h,w,dim), dtype=np.uint8 )
rgb[:,:,0] = rgb[:,:,1] = rgb[:,:,2] = a
print "rgb:", rgb
print "rgb.flags:", rgb.flags  # C_contiguous
print "rgb.strides:", rgb.strides  # (12, 3, 1)

dstack = np.dstack(( a, a, a ))
print "dstack:", dstack
print "dstack.flags:", dstack.flags  # F_contiguous
print "dstack.strides:", dstack.strides  # (1, 2, 8)

【讨论】：