¿Ofrece Cython una forma razonablemente fácil y eficiente de iterar los arreglos de Numpy como si fueran planos?

Digamos que quiero implementar Numpy’s

x[:] += 1 

en cython. Yo podria escribir

 @cython.boundscheck(False) @cython.wraparoundcheck(False) def add1(np.ndarray[np.float32_t, ndim=1] x): cdef unsigned long i for i in range(len(x)): x[i] += 1 

Sin embargo, esto solo funciona con ndim = 1 . Podría usar

 add1(x.reshape(-1)) 

pero esto solo funciona con x contiguas.

¿Ofrece Cython una forma razonablemente fácil y eficiente de iterar los arreglos de Numpy como si fueran planos?

( Reimplementar esta operación en particular en Cython no tiene sentido, ya que el código de Numpy anterior debería ser tan rápido como sea posible, solo uso esto como un simple ejemplo )

ACTUALIZAR:

Comparé las soluciones propuestas:

 @cython.boundscheck(False) @cython.wraparound(False) def add1_flat(np.ndarray x): cdef unsigned long i for i in range(x.size): x.flat[i] += 1 @cython.boundscheck(False) @cython.wraparound(False) def add1_nditer(np.ndarray x): it = np.nditer([x], op_flags=[['readwrite']]) for i in it: i[...] += 1 

La segunda función requiere import numpy as np además de cimport . Los resultados son:

 a = np.zeros((1000, 1000)) b = a[100:-100, 100:-100] %timeit b[:] += 1 1000 loops, best of 3: 1.31 ms per loop %timeit add1_flat(b) 1 loops, best of 3: 316 ms per loop %timeit add1_nditer(b) 1 loops, best of 3: 1.11 s per loop 

Por lo tanto, son 300 y 1000 veces más lentos que Numpy.

ACTUALIZACIÓN 2:

La versión add11 usa un bucle for dentro de un bucle for , por lo que no itera la matriz como si fuera plana. Sin embargo, es tan rápido como Numpy en este caso:

 %timeit add1.add11(b) 1000 loops, best of 3: 1.39 ms per loop 

Por otro lado, add1_unravel , una de las soluciones propuestas, no puede modificar el contenido de b .

http://docs.scipy.org/doc/numpy/reference/arrays.nditer.html

Es un buen tutorial sobre el uso de nditer . Termina con una versión de cython . nditer está destinado a ser el iterador de matrices de uso múltiple en el código de nivel numpy.

También hay buenos ejemplos de matrices en la página de vista de memoria de Cython:

http://cython.readthedocs.io/en/latest/src/userguide/memoryviews.html

http://docs.scipy.org/doc/numpy/reference/c-api.iterator.html

El búfer de datos de un ndarray es un búfer plano. Por lo tanto, independientemente de la forma y los pasos de la matriz, puede iterar sobre todo el búfer en forma de un puntero c plano. Pero cosas como nditer y memoryview se ocupan de los detalles del tamaño del elemento. Por lo tanto, en la encoding de nivel c , en realidad es más fácil recorrer todos los elementos de forma plana que paso a paso, al pasar por filas se debe tener en cuenta el paso de la fila.

Esto se ejecuta en Python, y creo que se traducirá bien en cython (no tengo esa configuración en mi máquina en este momento):

 import numpy as np def add1(x): it = np.nditer([x], op_flags=[['readwrite']]) for i in it: i[...] += 1 return it.operands[0] x = np.arange(10).reshape(2,5) y = add1(x) print(x) print(y) 

https://github.com/hpaulj/numpy-einsum/blob/master/sop.pyx es un script de sum de productos que escribí hace un tiempo para simular einsum .

El núcleo de su cálculo w = sop(x,y) es:

 it = np.nditer(ops, flags, op_flags, op_axes=op_axes, order=order) it.operands[nop][...] = 0 it.reset() for xarr, yarr, warr in it: x = xarr y = yarr w = warr size = x.shape[0] for i in range(size): w[i] = w[i] + x[i] * y[i] return it.operands[nop] 

===================

Al copiar ideas del documento nditer.html , obtuve una versión de add1 que es solo la mitad de la velocidad del número nativo a+1 . El ingenuo nditer (arriba) no es mucho más rápido en cython que en python . Gran parte de la aceleración puede deberse al external loop .

 @cython.boundscheck(False) def add11(arg): cdef np.ndarray[double] x cdef int size cdef double value it = np.nditer([arg], flags=['external_loop','buffered'], op_flags=[['readwrite']]) for xarr in it: x = xarr size = x.shape[0] for i in range(size): x[i] = x[i]+1.0 return it.operands[0] 

También codifiqué este nditer en python con una impresión de size , y encontré que iteraba en su b con 78 bloques de tamaño 8192, que es un tamaño de búfer, no una característica de b su diseño de datos.

 In [15]: a = np.zeros((1000, 1000)) In [16]: b = a[100:-100, 100:-100] In [17]: timeit add1.add11(b) 100 loops, best of 3: 4.48 ms per loop In [18]: timeit b[:] += 1 100 loops, best of 3: 8.76 ms per loop In [19]: timeit add1.add1(b) # for the unbuffered nditer 1 loop, best of 3: 3.1 s per loop In [21]: timeit add1.add11(a) 100 loops, best of 3: 5.44 ms per loop 

Puede intentar usar el atributo flat del ndarray , que proporciona un iterador sobre el objeto de matriz aplanada. Siempre itera en orden mayor de C, con el último índice variando el más rápido. Algo como:

 for i in range(x.size): x.flat[i] += 1 

Con numpy.ravel

numpy.ravel (a, orden = ‘C’)

Devuelve una matriz aplanada contigua.

Se devuelve una matriz 1-D, que contiene los elementos de la entrada. Solo se hace una copia si es necesario.

 @cython.boundscheck(False) @cython.wraparound(False) def add1_ravel(np.ndarray xs): cdef unsigned long i cdef double[::1] aview narray = xs.ravel() aview = narray for i in range(aview.shape[0]): aview[i] += 1 # return xs if the memory is shared if not narray.flags['OWNDATA'] or np.may_share_memory(xs, narray): return xs # otherwise new array reshaped shape = tuple(xs.shape[i] for i in range(xs.ndim)) return narray.reshape(shape)