用SWIG、cython或ctypes将包含数组成员的C结构封装以供Python访问？

Question

我想从Python中访问一个C语言的函数，这个函数返回一个包含双精度数组的结构体（这些数组的长度由结构体中的其他整型成员给出）。它的声明是：

typedef struct {
  int dim;
  int vertices;
  int quadrature_degree;
  int polynomial_degree;
  int ngi;
  int quadrature_familiy;
  double *weight; /* 1D: ngi */
  double *l;      /* 2D: ngi * dim */
  double *n;      /* 2D: ngi * vertices */
  double *dn;     /* 3D: ngi * vertices * dim */
} element;

extern void get_element(int dim, int vertices, int quad_degree, int poly_degree, element* e);

关键点是，我希望能够将所有的 double* 成员作为正确形状的NumPy数组来访问（也就是说，dn 应该可以作为一个三维数组来访问）。

简单地用SWIG包装这个结构体可以正常得到结构体，但所有的 double* 成员都变成了 <Swig Object of type 'double *' at 0x348c8a0>，这让它们变得没什么用。我尝试过修改NumPy的SWIG接口文件，但没有成功使用像 ( DATA_TYPE* INPLACE_ARRAY1, int DIM1 ) 这样的类型映射（我觉得在这种情况下可能无法匹配，但如果我错了也很乐意被纠正）。

我猜我需要手动初始化这些成员的NumPy数组为 PyArrayObject，然后让SWIG扩展我的结构体，以便在Python中访问它们？这看起来工作量很大。有没有人能想到更好的方法来使用SWIG？如果改变结构体或返回它的方法能让事情更简单，那也是可以的。

另外，我也看了看cython和ctypes。这些会更适合我想要实现的目标吗？我没有使用过cython，所以无法判断它的包装能力。对于ctypes，我大致能想象怎么做，但这意味着我得手动写出我原本希望自动化包装能帮我完成的工作。

任何建议都非常感谢！

Answer 1

看看SWIG的类型映射功能。它允许你为特定类型、特定实例（类型+名称）或者一组参数编写自己的处理代码。我还没有为结构体做过，但我可以给你一个例子，说明如何特别处理一个C函数，它需要一个数组和它的大小：

%typemap(in) (int argc, Descriptor* argv) {
    /* Check if is a list */
    if (PyList_Check($input)) {
        int size = PyList_Size($input);
        $1 = size;
        ...
        $2 = ...;
    }
}

这个代码会接收一对参数int argc, Descriptor* argv（因为参数名是给定的，所以它们也必须匹配），然后把你使用的PyObject传递给你，你可以在这里写任何你需要的C代码来进行转换。比如，你可以为double *dn做一个类型映射，利用NumPy的C API来完成转换。

Answer 2

Cython 的规则：

cdef extern from "the header.h":

ctypedef struct element:
  int dim
  int vertices
  int quadrature_degree
  int polynomial_degree
  int ngi
  int quadrature_familiy
  double *weight
  double *l
  double *n
  double *dn

void get_element(int dim, int vertices, int quad_degree, int poly_degree, element* e)

然后你可以在 Python 的环境中使用它。

Answer 3

使用SWIG的时候，需要为整个结构体定义一个类型映射。仅仅为指针成员定义类型映射是不够的，因为这样做没有上下文，无法知道要用什么大小来初始化NumPy数组。我通过以下的类型映射实现了我的需求（基本上是从numpy.i复制粘贴过来的，然后根据我的需要进行了调整，可能不是很稳健）：

%typemap (in,numinputs=0) element * (element temp) {
  $1 = &temp;
}

%typemap (argout) element * {
  /* weight */
  {
    npy_intp dims[1] = { $1->ngi };
    PyObject * array = PyArray_SimpleNewFromData(1, dims, NPY_DOUBLE, (void*)($1->weight));
    if (!array) SWIG_fail;
    $result = SWIG_Python_AppendOutput($result,array);
  }
  /* l */
  {
    npy_intp dims[2] = { $1->ngi, $1->dim };
    PyObject * array = PyArray_SimpleNewFromData(2, dims, NPY_DOUBLE, (void*)($1->l));
    if (!array) SWIG_fail;
    $result = SWIG_Python_AppendOutput($result,array);
  }
  /* n */
  {
    npy_intp dims[2] = { $1->ngi, $1->vertices };
    PyObject * array = PyArray_SimpleNewFromData(2, dims, NPY_DOUBLE, (void*)($1->n));
    if (!array) SWIG_fail;
    $result = SWIG_Python_AppendOutput($result,array);
  }
  /* dn */
  {
    npy_intp dims[3] = { $1->ngi, $1->vertices, $1->dim };
    PyObject * array = PyArray_SimpleNewFromData(3, dims, NPY_DOUBLE, (void*)($1->dn));
    if (!array) SWIG_fail;
    $result = SWIG_Python_AppendOutput($result,array);
  }
}

这个方法和C语言的函数不同，它返回的是一个包含我想要的数据的NumPy数组元组，这比之后从element对象中提取数据要方便得多。此外，第一个类型映射还消除了传入element类型对象的需要。因此，我可以完全将element结构体隐藏起来，不让Python用户看到。

最终，Python接口看起来是这样的：

weight, l, n, dn = get_element(dim, vertices, quadrature_degree, polynomial_degree)