C array a PyArray

Question

Escribo una extensión C de Python sin utilizar Cython.

Quiero asignar una matriz doble en C, usarla en una función interna (que está en Fortran) y devolverla. Señalo que la interfaz C-Fortran funciona perfectamente en C.

static PyObject * 
Py_drecur(PyObject *self, PyObject *args) 
{ 
    // INPUT 
    int n; 
    int ipoly; 
    double al; 
    double be; 

    if (!PyArg_ParseTuple(args, "iidd", &n, &ipoly, &al, &be)) 
    return NULL; 

    // OUTPUT 
    int nd = 1; 
    npy_intp dims[] = {n}; 
    double a[n]; 
    double b[n]; 
    int ierr; 

    drecur_(n, ipoly, al, be, a, b, ierr); 

    // Create PyArray 
    PyObject* alpha = PyArray_SimpleNewFromData(nd, dims, NPY_DOUBLE, a); 
    PyObject* beta = PyArray_SimpleNewFromData(nd, dims, NPY_DOUBLE, b); 

    Py_INCREF(alpha); 
    Py_INCREF(beta); 

    return Py_BuildValue("OO", alpha, beta); 
} 

I depurado el código y me da un fallo de segmentación cuando intento crear alfa de una. Hasta allí todo funciona bien. La función drecur_ funciona y obtengo el mismo problema si se elimina.

Ahora, ¿cuál es la forma estándar de definir un PyArray en torno a los datos C? Encontré documentación pero no un buen ejemplo. Además, ¿qué pasa con la pérdida de memoria? ¿Es correcto INCREF antes del retorno para que se conserve la instancia de alfa y beta? ¿Qué pasa con la desasignación cuando ya no son necesarios?

EDITAR finalmente lo hizo bien con el enfoque que se encuentra en NumPy cookbook.

static PyObject * 
Py_drecur(PyObject *self, PyObject *args) 
{ 
    // INPUT 
    int n; 
    int ipoly; 
    double al; 
    double be; 
    double *a, *b; 
    PyArrayObject *alpha, *beta; 

    if (!PyArg_ParseTuple(args, "iidd", &n, &ipoly, &al, &be)) 
    return NULL; 

    // OUTPUT 
    int nd = 1; 
    int dims[2]; 
    dims[0] = n; 
    alpha = (PyArrayObject*) PyArray_FromDims(nd, dims, NPY_DOUBLE); 
    beta = (PyArrayObject*) PyArray_FromDims(nd, dims, NPY_DOUBLE); 
    a = pyvector_to_Carrayptrs(alpha); 
    b = pyvector_to_Carrayptrs(beta); 
    int ierr; 

    drecur_(n, ipoly, al, be, a, b, ierr); 

    return Py_BuildValue("OO", alpha, beta); 
} 

double *pyvector_to_Carrayptrs(PyArrayObject *arrayin) { 
    int n=arrayin->dimensions[0]; 
    return (double *) arrayin->data; /* pointer to arrayin data as double */ 
} 

Siéntase libre de comentar sobre esto y gracias por las respuestas.

Answer 1

Un problema podría ser que sus matrices (a, b) tienen que durar al menos tanto como la matriz numpy que lo contiene. Ha creado sus matrices en ámbito local para que se destruyan cuando abandone el método.

Intenta que python asigne la matriz (por ejemplo, usando PyArray_SimpleNew), copia el contenido en ella y pásala por un puntero. También es posible que desee utilizar boost::python para cuidar estos detalles, si compilar contra refuerzo es una opción.

Answer 2

Por lo tanto, lo primero que parece sospechoso es que su matriz a y b están en el ámbito local de la función. Eso significa que después de la devolución obtendrá un acceso ilegal a la memoria.

Por lo que debe declarar las matrices con

double *a = malloc(n*sizeof(double)); 

Luego hay que asegurarse de que la memoria es luego liberado por el objeto que ha creado. Ver esta cita de la documentación:

PyObject PyArray_SimpleNewFromData (int nd, npy_intp se atenúa, int TypeNum, void * data)

A veces, usted quiere envolver memoria asignada a otra parte en un ndarray objeto para uso corriente abajo. Esta rutina hace que sea sencillo hacerlo. Los primeros tres argumentos son los mismos que en PyArray_SimpleNew, el argumento final es un puntero a un bloque de memoria contigua que el ndarray debe usar ya que es un búfer de datos que se interpretará en forma contigua al estilo C. Se devuelve una nueva referencia a un ndarray, pero el ndarray no será el propietario de sus datos. Cuando este ndarray se desasigna, el puntero no se liberará.

Debe asegurarse de que la memoria proporcionada no se libere mientras exista la matriz devuelta. La forma más fácil de manejar esto es si los datos provienen de otro objeto Python contado por referencia. El recuento de referencias en este objeto se debe aumentar después de que se pase el puntero, y el miembro base del ndarray devuelto debe apuntar al objeto de Python que posee los datos.Luego, cuando el ndarray es desasignado, el miembro base será DECREFADO apropiadamente. Si desea que la memoria se libere tan pronto como se desasigne el ndarray, simplemente configure el indicador OWNDATA en el ndarray devuelto.

Para su segunda cuestión, el Py_INCREF(alpha); es por lo general sólo es necesario si tiene la intención de mantener la referencia en una variable global o un miembro de la clase. Pero como solo está envolviendo una función, no tiene que hacerlo. Lamentablemente podría ser que la función PyArray_SimpleNewFromData no establezca el contador de referencia en 1, si ese fuera el caso, tendría que aumentarlo a 1. Espero que sea comprensible;).