Passaggio di array numpy tridimensionale a C

Question

Sto scrivendo un'estensione C al mio programma Python per scopi di velocità, e in esecuzione in un comportamento molto strano cercando di passare in un array numpy 3-dimensionale. Funziona con un array bidimensionale, ma sono sicuro che sto mettendo a punto qualcosa con i puntatori che cercano di farlo funzionare con la 3a dimensione. Ma ecco la parte strana. Se passo in un array 3-D, si blocca con un errore bus . Se (in Python) creo prima la mia variabile come array 2D, quindi la sovrascrivo con un array 3D, funziona perfettamente con. Se la variabile è prima un array vuoto e poi un array 3D, si blocca con un errore di segmento . Com'è possibile che ciò accada?

Inoltre, qualcuno può aiutarmi a far funzionare un array 3D? O dovrei semplicemente rinunciare e passare in un array 2D e rimodellarlo da solo?

Ecco il mio codice C:

static PyObject* func(PyObject* self, PyObject* args) { 
    PyObject *list2_obj; 
    PyObject *list3_obj; 
    if (!PyArg_ParseTuple(args, "OO", &list2_obj, &list3_obj)) 
    return NULL; 

    double **list2; 
    double ***list3; 

    //Create C arrays from numpy objects: 
    int typenum = NPY_DOUBLE; 
    PyArray_Descr *descr; 
    descr = PyArray_DescrFromType(typenum); 
    npy_intp dims[3]; 
    if (PyArray_AsCArray(&list2_obj, (void **)&list2, dims, 2, descr) < 0 || PyArray_AsCArray(&list3_obj, (void ***)&list3, dims, 3, descr) < 0) { 
    PyErr_SetString(PyExc_TypeError, "error converting to c array"); 
    return NULL; 
    } 
    printf("2D: %f, 3D: %f.\n", list2[3][1], list3[1][0][2]); 
} 

Ed ecco il mio codice Python che chiama la funzione di cui sopra:

import cmod, numpy 
l2 = numpy.array([[1.0,2.0,3.0], [4.0,5.0,6.0], [7.0,8.0,9.0], [3.0, 5.0, 0.0]]) 

l3 = numpy.array([[2,7, 1], [6, 3, 9], [1, 10, 13], [4, 2, 6]]) # Line A 
l3 = numpy.array([])            # Line B 

l3 = numpy.array([[[2,7, 1, 11], [6, 3, 9, 12]], 
       [[1, 10, 13, 15], [4, 2, 6, 2]]]) 

cmod.func(l2, l3) 

Quindi, se io commento sia la linea A e B, si blocca con un Errore del bus. Se la linea A è presente, ma la riga B è commentata, viene eseguita correttamente senza errori. Se la Linea B è presente ma la riga A è commentata, stampa i numeri corretti ma Segue errori. Infine, se entrambe le linee sono presenti, stampa anche i numeri corretti e quindi i difetti di Seg. Che diavolo sta succedendo qui?

MODIFICA: Ok. Wow. Quindi stavo usando int in Python ma li chiamavo double in C. E funzionava bene con gli array 1D e 2D. Ma non in 3D. Così ho cambiato la definizione Python di l3 per essere float, e ora funziona tutto in modo fantastico (Grazie mille Bi Rico).

Ma ora, più strano comportamento con Linee A & B! Ora se entrambe le linee sono commentate, il programma funziona. Se la linea B è presente ma A è commentata, funziona, e idem se entrambi sono non commentati. Ma se la linea A è presente e B è commentato, ottengo di nuovo il fantastico errore Bus. Mi piacerebbe davvero evitarli in futuro, quindi qualcuno ha la minima idea del perché la dichiarazione di una variabile Python possa avere questo tipo di impatto?

EDIT 2:. Beh, come folle come questi errori sono, sono tutti dovuti alla matrice NumPy 3-dimensionale passo in Se mi passa solo in 1 o 2-D array, si comporta come previsto, e la manipolazione delle altre variabili Python non fa nulla. Questo mi porta a credere che il problema si trovi da qualche parte nel conteggio dei riferimenti di Python. Nel codice C il conteggio dei riferimenti è diminuito più di quanto dovrebbe per gli array 3-D, e quando quella funzione restituisce Python prova a ripulire gli oggetti e tenta di eliminare un puntatore NULL. Questa è solo la mia ipotesi, e ho provato a Py_INCREF(); tutto quello che potevo pensare senza successo. Credo che mi limiterò a utilizzare una matrice 2D e rimodellare in C.

Answer 1

Ho già menzionato questo in un commento, ma spero che lo svuota un po 'aiuta a renderlo più chiaro.

Quando si lavora con gli array numpy in C, è bene essere espliciti sulla digitazione degli array. In particolare sembra che tu stia dichiarando i tuoi indicatori come double ***list3, ma in questo modo stai creando l3 nel tuo codice Python riceverai un array con il dtype npy_intp (penso). È possibile risolvere questo problema utilizzando esplicitamente il dtype durante la creazione degli array.

import cmod, numpy 
l2 = numpy.array([[1.0,2.0,3.0], 
        [4.0,5.0,6.0], 
        [7.0,8.0,9.0], 
        [3.0, 5.0, 0.0]], dtype="double") 

l3 = numpy.array([[[2,7, 1, 11], [6, 3, 9, 12]], 
        [[1, 10, 13, 15], [4, 2, 6, 2]]], dtype="double") 

cmod.func(l2, l3) 

Un'altra nota, a causa del modo in cui funziona pitone è quasi impossibile per "Linea A" e "linea B" per avere alcun effetto sul codice C che cosa così mai. So che questo sembra essere in conflitto con la tua esperienza empirica, ma sono abbastanza sicuro su questo punto.

Sono un po 'meno sicuro di ciò, ma in base alla mia esperienza con C, errori di bus e segofault non sono deterministici. Dipende dall'allocazione, dall'allineamento e dagli indirizzi di memoria. In alcune situazioni il codice sembra funzionare benissimo 10 volte e fallisce nell'11a esecuzione anche se non è cambiato nulla.

Avete considerato l'utilizzo di cython? So che non è un'opzione per tutti, ma se si tratta di un'opzione è possibile ottenere velocità quasi C con typed memoryviews.

Answer 2

Secondo http://docs.scipy.org/doc/numpy/reference/c-api.array.html?highlight=pyarray_ascarray#PyArray_AsCArray:

Nota La simulazione di un array di C-stile non è completa per 2-D e 3- d matrici. Ad esempio, gli array simulati di puntatori non possono essere passati a subroutine che prevedono array 2-d e 3-d specifici e definiti staticamente. Per passare a funzioni che richiedono quel tipo di input, è necessario definire staticamente la matrice richiesta e copiare i dati.

Penso che questo significhi che PyArray_AsCArray restituisca un blocco di memoria con i dati in esso in ordine C. Tuttavia, per accedere a tali dati, sono necessarie ulteriori informazioni (vedere http://www.phy225.dept.shef.ac.uk/mediawiki/index.php/Arrays,_dynamic_array_allocation). Ciò può essere ottenuto conoscendo le dimensioni in anticipo, dichiarando un array e quindi copiando i dati nell'ordine corretto. Tuttavia, sospetto che il caso più generale sia più utile: non conosci le dimensioni finché non vengono restituite. Penso che il codice seguente creerà la struttura del puntatore C necessaria per consentire l'indirizzamento dei dati.

static PyObject* func(PyObject* self, PyObject* args) { 
    PyObject *list2_obj; 
    PyObject *list3_obj; 
    if (!PyArg_ParseTuple(args, "OO", &list2_obj, &list3_obj)) return NULL; 

    double **list2; 
    double ***list3; 

    // For the final version 
    double **final_array2; 
    double **final_array2; 

    // For loops 
    int i,j; 

    //Create C arrays from numpy objects: 
    int typenum = NPY_DOUBLE; 
    PyArray_Descr *descr; 
    descr = PyArray_DescrFromType(typenum); 

    // One per array coming back ... 
    npy_intp dims2[2]; 
    npy_intp dims3[3]; 

    if (PyArray_AsCArray(&list2_obj, (void **)&list2, dims2, 2, descr) < 0 || PyArray_AsCArray(&list3_obj, (void ***)&list3, dims3, 3, descr) < 0) { 
     PyErr_SetString(PyExc_TypeError, "error converting to c array"); 
     return NULL; 
    } 

    // Create the pointer arrays needed to access the data 

    // 2D array 
    final_array2 = calloc(dim2[0], sizeof(double *)); 
    for (i=0; i<dim[0]; i++) final_array2[i] = list2 + dim2[1]*sizeof(double); 

    // 2D array 
    final_array3 = calloc(dim3[0], sizeof(double **)); 
    final_array3[0] = calloc(dim3[0]*dim3[1], sizeof(double *)); 
    for (i=0; i<dim[0]; i++) { 
     final_array3[i] = list2 + dim3[1]*sizeof(double *); 
     for (j=0; j<dim[1]; j++) { 
      final_array[i][j] = final_array[i] + dim3[2]*sizeof(double); 
     } 
    } 

    printf("2D: %f, 3D: %f.\n", final_array2[3][1], final_array3[1][0][2]); 
    // Do stuff with the arrays 

    // When ready to complete, free the array access stuff 
    free(final_array2); 

    free(final_array3[0]); 
    free(final_array3); 

    // I would guess you also need to free the stuff allocated by PyArray_AsCArray, if so: 
    free(list2); 
    free(list3); 
} 

non riuscivo a trovare una definizione per npy_intp, la assume sopra di esso è lo stesso di int. In caso contrario, sarà necessario convertire dim2 e dim3 in matrici int prima di eseguire il codice.

Answer 3

Invece di convertire in un array in stile c, di solito accedo agli elementi di array numpy direttamente utilizzando PyArray_GETPTR (vedere http://docs.scipy.org/doc/numpy/reference/c-api.array.html#data-access).

Ad esempio, per accedere a un elemento di una matrice numpy tridimensionale di tipo double use double elem=*((double *)PyArray_GETPTR3(list3_obj,i,j,k)).

Per l'applicazione, è possibile rilevare il numero corretto di dimensioni per ciascun array utilizzando PyArray_NDIM, quindi accedere agli elementi utilizzando la versione appropriata di PyArray_GETPTR.

Answer 4

C'era un bug nel NumPy C-API, che dovrebbe essere fissato Now:

https://github.com/numpy/numpy/pull/5314