11

我有一个 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 中的某种数组顺序检查是可行的方法。当然，我更喜欢不需要将数据复制到新数组中的解决方案，但我几乎认为这是必要的。

4

3 回答 3

8

如果您想在不复制的情况下支持 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*ncol + col]C 代码中出现的每个mat[row*rowstride + col*colstride。

于 2010-12-12T15:53:23.690 回答
3

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

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

另外，考虑看看 Numpy 的确切数组接口（也有 C 部分）。

于 2010-12-12T13:41:16.117 回答
0

+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)
于 2011-05-07T13:25:00.510 回答