如何在用C扩展numpy时处理列连续数组

Question

我有一个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数组都是行连续的，所以这个函数运行得很好。然而，如果一个numpy数组之前被转置过，数据并没有被复制，而是返回了一个新的视图。在这种情况下，我的函数就会失败，因为数组不再是行连续的了。

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

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

显然，这样做很容易出错，用户必须确保数组的顺序是正确的，而这本来不应该是Python用户需要关心的事情。

那么，有什么好的方法可以让这个函数同时支持行连续和列连续的数组呢？我想在Cython中进行某种数组顺序的检查可能是个不错的主意。当然，我更希望有一个不需要将数据复制到新数组中的解决方案，但我几乎认为这可能是必要的。

Answer 1

给Sven点赞，他的回答解决了一个让我困惑的问题，那就是dstack这个函数返回的是一个F连续的数组？!

# 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)

Answer 2

在这种情况下，你确实需要创建一个输入数组的副本（这个数组可以是一个真实数组的视图），并确保它的行是连续的。你可以通过类似下面的方式来实现：

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

另外，建议你查看一下Numpy的具体数组接口，里面还有C语言的部分。

Answer 3

如果你想在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

接下来，把你C代码中每次出现的 mat[row*ncol + col] 替换成 mat[row*rowstride + col*colstride]。