2010-12-12 4 views
9

J'ai une fonction C pour normaliser les rangées d'un tableau dans l'espace de journal (cela empêche le sous-dénombrement numérique).Comment comptabiliser un tableau contigu à une colonne lors de l'extension de numpy avec C

Le prototype de mon C-fonction est la suivante:

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

Vous pouvez voir que cela prend un pointeur sur un tableau et il modifie en place. Le code C suppose bien entendu que les données sont sauvegardées en tant que matrice C-contiguë, c'est-à-dire contiguë à une ligne.

J'enveloppe la fonction comme suit à l'aide cython (importations et cdef extern from omis):

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 

La plupart du temps numpy-réseaux sont rangées contiguës et la fonction fonctionne très bien. Cependant, si un tableau numérique a déjà été transposé, les données ne sont pas copiées mais seulement une nouvelle vue dans les données est retournée. Dans ce cas, ma fonction échoue car le tableau n'est plus contigu à la ligne.

Je peux contourner ce problème en définissant le tableau d'avoir pour Fortran contiguës, de sorte qu'après transposition sera C-contiguë:

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

Il est évident que c'est très sujette à l'erreur et l'utilisateur doit prendre Veillez à ce que le tableau soit dans le bon ordre, ce qui est quelque chose dont l'utilisateur ne devrait pas avoir besoin dans Python.

Quelle est la meilleure façon de faire en sorte que cela fonctionne avec les tableaux contigus à la fois en ligne et en colonne? Je suppose que la vérification des commandes de tableaux dans Cython est la bonne façon de procéder. Bien sûr, je préfère une solution qui ne nécessite pas de copier les données dans un nouveau tableau, mais je suppose presque que c'est nécessaire.

Répondre

7

Si vous souhaitez prendre en charge des tableaux en ordre C et Fortran sans jamais les copier, votre fonction C doit être suffisamment flexible pour prendre en charge les deux ordres. Cela peut être réalisé en faisant passer les progrès du tableau NumPy à la fonction C: Changer le prototype à

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

et l'appel Cython à

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 

Ensuite, remplacer toutes les occurrences de mat[row*ncol + col] dans votre C code par mat[row*rowstride + col*colstride].

+0

Cette réponse de 2010 est-elle toujours à jour ou existe-t-il un meilleur moyen d'y parvenir maintenant? –

+0

@larsmans: Je ne sais pas exactement ce que vous entendez par "ceci".L'écriture d'une fonction C capable de traiter à la fois les tableaux bidimensionnels Fortran-contigus et C-contigus fonctionne toujours de cette façon, si c'est ce que vous voulez. Si tout va bien que vos tableaux soient copiés, il existe (et a été en 2010) d'autres solutions. –

2

Dans ce cas, vous ne voulez vraiment créer une copie du tableau d'entrée (qui peut être vue sur un réel tableau) avec ordre de rangée contiguë garantie. Vous pouvez y parvenir avec quelque chose comme ceci:

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

Pensez aussi à jeter un oeil à la array interface exacte de Numpy (il fait partie C aussi).

0

+1 à Sven, dont la réponse résout le Gotcha (bien, m'a fait) que dstack retourne un tableau 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) 
Questions connexes