Wrap C struct с элементом массива для доступа в python: SWIG? Cython? ctypes?

Я хочу получить доступ к функции C, которая возвращает структуру, содержащую двойные массивы (где длины этих массивов задаются другими членами int структуры) из python. Объявление

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-wrapping дает мне структуру просто отлично, но все члены 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 я могу представить себе, как это сделать, но это означает, что я пишу вручную то, что, как я надеялся, может сделать для меня достаточно автоматическая обертка.

Любые предложения с благодарностью получили!

Ответ 1

Использование SWIG требует типовой карты для всей структуры. Tyepmaps только для членов указателя недостаточно, так как они не имеют контекста, чтобы узнать, какой размер для инициализации массивов 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)

Ответ 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

Ответ 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, которая будет использовать API NumPy C для преобразования.

Ответ 4

Вы всегда можете написать вспомогательные функции, которые принимают "элемент *" и возвращают искомый элемент:

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

Если вам нужно изменить, а также прочитать, вам понадобятся отдельные "геттеры" и "сеттеры", конечно.

SWIG должен иметь возможность легко переносить все эти файлы и выставлять их на Python.

Производительность может быть невелика, но, вероятно, не хуже альтернатив.

Ответ 5

Эквивалент созданному SWIG-модулю с использованием 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