在Python中访问包含数组成员的C结构体：使用SWIG？Cython？ctypes？

Question

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

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应该作为3D数组可访问）。

简单地用SWIG包装可以很好地给出结构体，但是所有的double*成员都是<Swig Object of type 'double *' at 0x348c8a0>，这使它们无用。我尝试过使用NumPy SWIG接口文件进行操作，但无法使任何类型映射（例如(DATA_TYPE* INPLACE_ARRAY1, int DIM1)）工作（我认为在这种情况下不可能匹配，但如果有证据证明我错了，我会很高兴）。

我猜想我需要手动编写初始化NumPy数组的代码作为这些成员的PyArrayObject，并扩展我的结构体以使它们在Python中可访问？这看起来像是很多工作。有人能否看到使用SWIG更好的方法？如果这些可以使事情变得更容易，那么更改结构或返回它的方法也是可行的。

另外，我看了看cython和ctypes。它们是否更适合我想要实现的东西？我没有使用过cython，所以无法评估它的包装能力。对于ctypes，我可以大致想象如何实现，但这意味着手动编写我原本希望由相当自动化的包装器为我完成的内容。

非常感谢任何建议！

- kynan

---

如何运行这段代码？我尝试使用 get_element(1,2,3,6)，但出现了错误 ValueError: 数组太大；arr.size * arr.dtype.itemsize 大于最大可能的大小。 - Pygirl

Answer 1

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空间中与之进行交互。

- fabrizioM

Answer 2

使用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 对象中提取数据更方便。第一个 typemap 还消除了需要传入 element 类型对象的需要，因此我可以完全将 element 结构体隐藏在 Python 用户之外。

Python 接口的最终形式如下:

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

- kynan

---

你能在这里展示完整的接口文件吗？我没有得到预期的结果。 - xgdgsc

@xgdgsc 你可以在这里找到它。 - kynan

谢谢！现在我的Python程序在访问其中一个返回的结果时崩溃了，造成这种情况的原因是什么？ - xgdgsc

@xgdgsc，你的程序崩溃可能有很多原因，例如如果C函数返回的结构体没有正确初始化，但我不认为这是讨论的地方。此外，我不再使用SWIG，转而使用Cython（请参见Fabrizio的答案）。 - kynan

Answer 3

可以看看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进行转换。

- Adam

---

double *dn等的typemap的问题在于没有办法访问大小信息。不过，element *的typemap解决了这个问题，请参见我的答案。 - kynan

Answer 4

你可以编写帮助函数，接受一个“element *”并返回你所寻找的元素：

double element_get_weight(const element *elt, unsigned n) {
    assert(n < elt->ngi);  /* or similar */
    return elt->weight[n];
}

如果您需要修改和读取，当然需要单独的“getter”和“setter”。

SWIG 应该能够轻松地包装所有这些并将它们暴露给 Python。

性能可能不是很好，但可能不会比其他选择差。

- Nemo

---

这将是我的最后选择，因为添加这些getter和setter非常麻烦，而且从Python中使用它们更加麻烦。性能并不是真正的问题。 - kynan

Answer 5

使用ctypes创建的等效模块如下所示：

from ctypes import *
from numpy import *

lib = cdll.LoadLibrary("_get_element.so")

class ELEMENT(Structure):
    _fields_ = [("dim", c_int),
                ("vertices", c_int),
                ("quadrature_degree", c_int),
                ("polynomial_degree", c_int),
                ("ngi", c_int),
                ("quadrature_familiy", c_int),
                ("weight", POINTER(c_double)),
                ("l", POINTER(c_double)),
                ("n", POINTER(c_double)),
                ("dn", POINTER(c_double))]

cget_element = lib.get_element
cget_element.argtypes = [c_int, c_int, c_int, c_int, POINTER(ELEMENT)]
cget_element.restype = None

def get_element(dim, vertices, quad_degree, poly_degree):
    e = ELEMENT()
    cget_element(dim, vertices, quad_degree, poly_degree, byref(e))
    weight = asarray([e.weight[i] for i in xrange(e.ngi)], dtype=float64)
    l = asarray([e.l[i] for i in xrange(e.ngi*e.dim)], dtype=float64).reshape((e.ngi,e.dim))
    n = asarray([e.n[i] for i in xrange(e.ngi*e.vertices)], dtype=float64).reshape((e.ngi,e.vertices))
    dn = asarray([e.dn[i] for i in xrange(e.ngi*e.vertices*e.dim)], dtype=float64).reshape((e.ngi,e.vertices,e.dim))
    return weight, l, n, dn

- kynan