首页 文章

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

提问于
浏览 次
10

我有一个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

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*rowstride + col*colstride ]替换C代码中 mat[row*ncol + col] 的每个出现 .

    回复于
  • 2

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

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

    另外，考虑一下Numpy的确切array interface（也有C部分） .

    回复于
  • 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)
    回复于

相关问题