Wrap C struct with array member for access in python: SWIG? cython? ctypes?

Question

I want to access a C function that returns a struct containing double arrays (where the lengths of these arrays is given by other int members of the struct) from python. The declaration is

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);

The important point is I want to be able to access all the double* members as NumPy arrays of the correct shape (i.e. dn should be a accessible as 3D array).

Simply SWIG-wrapping this gives me the struct just fine, but all the double* members are <Swig Object of type 'double *' at 0x348c8a0> which makes them useless. I played around with the NumPy SWIG interface file but couldn't get any of the typemaps like ( DATA_TYPE* INPLACE_ARRAY1, int DIM1 ) to work (I think it's not possible to get them to match in this case but I'd be happy to be proven wrong).

My guess is I'd have to hand code initialization of the NumPy arrays as PyArrayObject for these members and SWIG extend my struct to make them accessible in Python? That looks like a lot of work. Can anyone see a nicer way using SWIG? It would be possible to change the struct or the method returning it if that made things easier.

Alternatively I had a look at cython and ctypes. Would these be better suited for what I'm trying to achieve? I haven't used cython so can't judge it's wrapping capabilities. For ctypes I can roughly imagine how to do it, but it means writing by hand what I had hoped a reasonably automated wrapper could do for me.

Any suggestions gratefully received!

Answer 1

Cython rules:

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)

and then you can interface it, from python space

Answer 2

Using SWIG requires a typemap for the entire struct. Tyepmaps for only the pointer members are not enough, since they don't have the context to know what size to initialize the NumPy arrays with. I managed to get what I wanted with the following typemaps (which was basically copy & paste from numpy.i and adapt to my needs, probably not very robust):

%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);
  }
}

This works different from the C function in that it returns a tuple of NumPy arrays with the data I want, which is more convenient than having to extract it from the element object later. The first typemap furthermore eliminates the need to pass in an object of type element. Hence I can hide the element struct entirely from the python user.

The python interface finally looks like this:

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

Answer 3

Check out SWIG's typemaps. They let you write your own code for handling specific types, specific instances (type+name) or even groups of arguments. I haven't done it for structures, but to specially handle a case where the C function takes an array and its size:

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

That will take the pair of arguments int argc, Descriptor* argv (since the names are provided they have to match as well) and pass you the PyObject used and you write whatever C code you need to do the conversion. You could do a typemap for double *dn that would use the NumPy C API to do the conversion.