Array C in PyArray

Question

Sto scrivendo un'estensione C Python senza utilizzare Cython.

Desidero assegnare un doppio array in C, utilizzarlo in una funzione interna (che si trova in Fortran) e restituirlo. Faccio notare che l'interfaccia C-Fortran funziona perfettamente in 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); 
} 

ho il debug questo codice e ottengo un errore di segmentazione quando provo a generare alfa su un. Fino a lì tutto funziona bene. La funzione drecur_ funziona e ottengo lo stesso problema se viene rimosso.

Ora, qual è il modo standard per definire un PyArray attorno ai dati C? Ho trovato documentazione ma nessun buon esempio. Inoltre, che dire della perdita di memoria? È corretto per INCREF prima del ritorno in modo che l'istanza di alpha e beta siano preservate? Che dire della deallocazione quando non sono più necessari?

EDIT finalmente ho proprio con l'approccio trovata in 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 */ 
} 

Sentitevi liberi di commentare questo e grazie per le risposte.

Answer 1

Un problema potrebbe essere che gli array (a, b) devono durare almeno quanto l'array numpy che lo contiene. Hai creato i tuoi array in ambito locale in modo che vengano distrutti quando lasci il metodo.

Provare ad avere python che assegna l'array (ad esempio utilizzando PyArray_SimpleNew), copiare il contenuto in esso e passarlo a un puntatore. Si potrebbe anche voler usare boost::python per aver cura di questi dettagli, se costruire contro boost è un'opzione.

Answer 2

Quindi le prime cose che sembrano sospette, è che l'array a e b si trova nell'ambito locale della funzione. Ciò significa che dopo il ritorno si otterrà un accesso illegale alla memoria.

Così si dovrebbe dichiarare le matrici con

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

Quindi è necessario assicurarsi che la memoria viene successivamente liberato da l'oggetto è stato creato. Vedere questa citazione della documentazione:

PyObject PyArray_SimpleNewFromData (int nd, npy_intp affievolisce, int TypeNum, dati void *)

A volte, si vuole avvolgere memoria allocata altrove in un ndarray oggetto per l'uso a valle. Questa routine rende semplice farlo. I primi tre argomenti sono gli stessi di PyArray_SimpleNew, l'argomento finale è un puntatore a un blocco di memoria contigua che il ndarray dovrebbe usare come buffer di dati che sarà interpretato in modo contiguo in stile C. Viene restituito un nuovo riferimento a un narray, ma ndarray non ne possiede i dati. Quando questo ndarray è deallocato, il puntatore non verrà liberato.

È necessario assicurarsi che la memoria fornita non venga liberata mentre l'array restituito è esistente. Il modo più semplice per gestirlo è se i dati provengono da un altro oggetto Python conteggiato di riferimento. Il conteggio dei riferimenti su questo oggetto deve essere aumentato dopo il passaggio del puntatore e il membro di base del parametro ndarray restituito deve puntare all'oggetto Python che possiede i dati.Quindi, quando il narray viene deallocato, il membro di base sarà DECREF appropriato. Se si desidera liberare la memoria non appena viene deallocato ndarray, è sufficiente impostare il flag OWNDATA sul ndarray restituito.

riguarda la seconda domanda la Py_INCREF(alpha); è generalmente necessaria solo se si intende mantenere il riferimento in una variabile globale o un membro della classe. Ma dal momento che si esegue il wrapping di una funzione, non è necessario. Purtroppo potrebbe essere che la funzione PyArray_SimpleNewFromData non imposti il ​​contatore di riferimento su 1, se questo fosse il caso dovresti aumentarlo a 1. Spero che sia comprensibile;).