2010-12-12 10 views
9

Ho una funzione C per normalizzare le righe di un array nello spazio log (ciò impedisce l'underflow numerico).Come rendere conto di un array di colonne contigui durante l'estensione di numpy con C

Il prototipo della mia C-funzione è la seguente:

void normalize_logspace_matrix(size_t nrow, size_t ncol, double* mat); 

Si può vedere che ci vuole un puntatore ad un array e lo modifica in posizione. Il codice C ovviamente presuppone che i dati vengano salvati come un array contiguo C, cioè contiguo alla riga.

I avvolgere la funzione come segue utilizzando Cython (importazioni e cdef extern from omesso):

def normalize_logspace(np.ndarray[np.double_t, ndim=2] mat): 
    cdef Py_ssize_t n, d 
    n = mat.shape[0] 
    d = mat.shape[1] 
    normalize_logspace_matrix(n, d, <double*> mat.data) 
    return mat 

maggior parte del tempo NumPy-array sono riga contigua e la funzione funziona bene. Tuttavia, se un array numpy è stato precedentemente trasposto, i dati non vengono copiati ma viene restituita solo una nuova vista nei dati. In questo caso la mia funzione fallisce perché l'array non è più contiguo alla riga.

posso aggirare questo definendo la matrice per avere ordine Fortran contigue, in modo tale che dopo la trasposizione sarà C contigue:

A = np.array([some_func(d) for d in range(D)], order='F').T 
A = normalize_logspace(A) 

Ovviamente questo è molto soggetto ad errori e l'utente deve prendere cura che l'array sia nell'ordine corretto, cosa che l'utente non dovrebbe aver bisogno di preoccuparsi di Python.

Qual è il modo migliore in cui è possibile farlo funzionare con array sia in righe che in colonna? Presumo che una sorta di controllo degli array in Cython sia la strada da percorrere. Certo, preferirei una soluzione che non richieda di copiare i dati in un nuovo array, ma presumo quasi che sia necessario.

risposta

7

Se si desidera supportare gli array in ordine C e Fortran senza mai copiare, la funzione C deve essere sufficientemente flessibile da supportare entrambi gli ordini. Ciò può essere ottenuto facendo passare i passi della matrice NumPy alla funzione C: Cambiare il prototipo per

void normalize_logspace_matrix(size_t nrow, size_t ncol, 
           size_t rowstride, size_t colstride, 
           double* mat); 

e la chiamata Cython a

def normalize_logspace(np.ndarray[np.double_t, ndim=2] mat): 
    cdef Py_ssize_t n, d, rowstride, colstride 
    n = mat.shape[0] 
    d = mat.shape[1] 
    rowstride = mat.strides[0] // mat.itemsize 
    colstride = mat.strides[1] // mat.itemsize 
    normalize_logspace_matrix(n, d, rowstride, colstride, <double*> mat.data) 
    return mat 

Poi, sostituire ogni occorrenza di mat[row*ncol + col] in C codice da mat[row*rowstride + col*colstride].

+0

Questa risposta del 2010 è ancora attuale oppure esiste un modo migliore per ottenere questo risultato ora? –

+0

@larsmans: non so esattamente cosa intendi con "questo".Scrivere una funzione C in grado di gestire sia gli array bidimensionali Fortran-contigui che C-contigui funziona ancora in questo modo, se questo è ciò che si desidera. Se è ok che i tuoi array vengano copiati, ci sono (e sono stati nel 2010) altre soluzioni. –

2

In questo caso è davvero si vuole creare una copia dell'array di ingresso (che può essere una vista su un vero e proprio gamma ) con ordine di fila contigue garantito. È possibile raggiungere questo obiettivo con qualcosa di simile:

a = numpy.array(A, copy=True, order='C') 

Inoltre, prendere in considerazione uno sguardo alla esatto array interface di Numpy (c'è C parte troppo).

0

+1 a Sven, la cui risposta risolve il Gotcha (beh, mi ha fatto) che dstack restituisce un array F_contiguous?!

# don't use dstack to stack a,a,a -> rgb for a C func 

import sys 
import numpy as np 

h = 2 
w = 4 
dim = 3 
exec("\n".join(sys.argv[1:])) # run this.py h= ... 

a = np.arange(h*w, dtype=np.uint8) .reshape((h,w)) 
rgb = np.empty((h,w,dim), dtype=np.uint8) 
rgb[:,:,0] = rgb[:,:,1] = rgb[:,:,2] = a 
print "rgb:", rgb 
print "rgb.flags:", rgb.flags # C_contiguous 
print "rgb.strides:", rgb.strides # (12, 3, 1) 

dstack = np.dstack((a, a, a)) 
print "dstack:", dstack 
print "dstack.flags:", dstack.flags # F_contiguous 
print "dstack.strides:", dstack.strides # (1, 2, 8) 
Problemi correlati