Python / numpy problema complicado rebanar

Question

Tengo un problema con algunas cosas numpy. Necesito una matriz numpy a comportarse de una manera inusual, devolviendo una rebanada como una vista de los datos que he rodajas, no una copia. Así que aquí está un ejemplo de lo que quiero hacer:

Supongamos que tenemos una matriz simple como esto:

a = array([1, 0, 0, 0])

Me gustaría actualizar entradas consecutivas de la matriz (que se mueve de izquierda a derecha) con la entrada anterior de la matriz, utilizando la sintaxis siguiente:

a[1:] = a[0:3]

Esto tendría el siguiente resultado:

a = array([1, 1, 1, 1])

O algo así:

a[1:] = 2*a[:3]
# a = [1,2,4,8]

Para ilustrar aún más quiero que los siguientes tipos de comportamiento:

for i in range(len(a)):
    if i == 0 or i+1 == len(a): continue
    a[i+1] = a[i]

A excepción Quiero que la velocidad de numpy.

El comportamiento predeterminado de numpy es tomar una copia de la rebanada, así que lo que en realidad me da es la siguiente:

a = array([1, 1, 0, 0])

Ya tiene esta matriz como una subclase de la ndarray, por lo que puede realizar otros cambios en él si es necesario, sólo necesito el corte en el lado derecho de ser actualizada continuamente a medida que se actualiza la rebanada en la mano izquierda lado.

¿Estoy soñando? Es posible esta magia?

Actualización: Todo esto es porque estoy tratando de utilizar Gauss-Seidel iteración para resolver un problema de álgebra lineal, más o menos. Es un caso especial que implica funciones armónicas, yo estaba tratando de evitar entrar en esto porque no es realmente necesario y es probable que las cosas confundir aún más, pero aquí va.

El algoritmo es el siguiente:

while not converged:
    for i in range(len(u[:,0])):
        for j in range(len(u[0,:])):
            # skip over boundary entries, i,j == 0 or len(u)
            u[i,j] = 0.25*(u[i-1,j] + u[i+1,j] + u[i, j-1] + u[i,j+1])

¿verdad? Pero se puede hacer de dos maneras, Jacobi implica la actualización de cada elemento con sus vecinos sin tener en cuenta las actualizaciones que ya haya hecho hasta que los ciclos de bucle while, que hacerlo en bucles que copia la matriz a continuación, actualizar una matriz a partir de la matriz copiado. Sin embargo Gauss-Seidel utiliza la información ya ha actualizado para cada uno de los i-1 y las entradas j-1, por lo tanto no hay necesidad de una copia, el bucle debe esencialmente 'sabe' desde la matriz ha sido re-evaluado después de cada actualización solo elemento . Es decir, cada vez que llamamos a una entrada como u [i-1, j] o T [i, j-1] la información calculada en el circuito anterior a estar allí.

quiero reemplazar esta situación bucle anidado lento y feo, con una bonita línea limpia de código usando rebanar numpy:

u[1:-1,1:-1] = 0.25(u[:-2,1:-1] + u[2:,1:-1] + u[1:-1,:-2] + u[1:-1,2:])

Sin embargo, el resultado es Jacobi iteración porque cuando se toma una rebanada: [u:, - 2,1: -1] copiar los datos, por lo que el corte no está al tanto de las actualizaciones hechas. Ahora numpy todavía bucles ¿verdad? Su no paralelo es sólo una manera más rápida de bucle que se parece a una operación paralela en Python. Quiero aprovechar este comportamiento por tipo de piratería numpy para devolver un puntero en lugar de una copia cuando tomo una rebanada. ¿Derecho? A continuación, cada bucles numpy tiempo, esa rebanada voluntad 'actualización' o en realidad sólo la réplica lo que sucedió en la actualización. Para hacer esto necesito rebanadas en ambos lados de la matriz para ser punteros.

De todos modos, si hay alguna persona realmente muy inteligente por ahí que impresionante, pero he más o menos a mí mismo renunció a creer que la única respuesta es un bucle en C.

Answer 1

respuesta tardía, pero esto resultó en Google, así que probablemente señalar que el documento OP quería. Su problema es claro: cuando se utiliza rebanadas NumPy, los temporales se crean. Envuelve el código de una llamada rápida a weave.blitz para deshacerse de los provisionales y tienen el comportamiento de su necesidad.

Leer la sección de PerformancePython tutorial para más detalles.

Answer 2

acumulan está diseñado para hacer lo que parece querer; Es decir, para proprigate una operación a lo largo de una matriz. He aquí un ejemplo:

from numpy import *

a = array([1,0,0,0])
a[1:] = add.accumulate(a[0:3])
# a = [1, 1, 1, 1]

b = array([1,1,1,1])
b[1:] = multiply.accumulate(2*b[0:3])
# b = [1 2 4 8]

Otra manera de hacer esto es especificar explícitamente la matriz resultado que la matriz de entrada. He aquí un ejemplo:

c = array([2,0,0,0])
multiply(c[:3], c[:3], c[1:])
# c = [  2   4  16 256]

Answer 3

Sólo tiene que utilizar un bucle. No puedo pensar inmediatamente de cualquier manera de hacer que el operador de parte comportan de la manera que usted está diciendo que usted quiere que, a excepción de tal vez subclasificando array de numpy y reemplazando el método apropiado con algún tipo de Python vudú ... pero lo más importante, la idea de que a[1:] = a[0:3] debe copiar el primer valor de a en los próximos tres ranuras parece completamente sin sentido para mí. Imagino que nadie podía confundir fácilmente otra persona que mira a su código (por lo menos las primeras veces).

Answer 4

No es la lógica correcta. Voy a tratar de utilizar letras para explicarlo.

array = abcd Imagen con a, b, c, d como elementos.
Ahora, medios array[1:] desde el elemento en la posición 1 (a partir de 0) en.
En este caso: bcd y array[0:3] medio del carácter en la posición 0 hasta el tercer personaje (el que está en posición 3-1) en este caso:. 'abc'

Escritura algo como:
array[1:] = array[0:3]

medios: sustituir bcd con abc

Para obtener la salida que desea, ahora en pitón, usted debe usar algo como:

a[1:] = a[0]

Answer 5

Se debe tener algo que ver con la asignación de una rebanada. Los operadores, sin embargo, como usted ya sabe, no siguen el comportamiento esperado:

>>> a = numpy.array([1,0,0,0])
>>> a[1:]+=a[:3]
>>> a
array([1, 1, 1, 1])

Si ya tiene ceros en el problema del mundo real donde su ejemplo lo hace, entonces esto lo resuelve. De lo contrario, a un costo adicional, que se establece en cero, ya sea mediante la multiplicación por cero o asignar a cero, (el que sea más rápido)

editar: Tuve otro pensamiento. Es posible que prefiera esto:

numpy.put(a,[1,2,3],a[:3]) 

Answer 6

Numpy debe comprobar si la matriz de destino es la misma que la matriz de entrada cuando se hace la llamada setkey. Afortunadamente, hay maneras alrededor de él. En primer lugar, he intentado usar numpy.put lugar

In [46]: a = numpy.array([1,0,0,0])

In [47]: numpy.put(a,[1,2,3],a[0:3])

In [48]: a
Out[48]: array([1, 1, 1, 1])

Y a continuación, a partir de la documentación que, di usando flatiters intentarlo ( a.flat )

In [49]: a = numpy.array([1,0,0,0])

In [50]: a.flat[1:] = a[0:3]

In [51]: a
Out[51]: array([1, 1, 1, 1])

Pero esto no resuelve el problema que tenía en mente

In [55]: a = np.array([1,0,0,0])

In [56]: a.flat[1:] = 2*a[0:3]

In [57]: a
Out[57]: array([1, 2, 0, 0])

Esta falla debido a la multiplicación se hace antes de la asignación, no en paralelo como le gustaría.

Numpy está diseñado para la aplicación repetida de la misma operación exacta en paralelo a través de una matriz. Para hacer algo más complicado, a menos que pueda encontrar descomponerlo en términos de funciones como numpy.cumsum y numpy.cumprod, tendrá que recurrir a algo como scipy.weave o escribir la función en C. (Véase página PerfomancePython para más detalles.) (Además, nunca he utilizado la armadura, así que no puedo garantizar que va a hacer lo que quiere .)

Answer 7

Se puede echar un vistazo a np.lib.stride_tricks.

Hay alguna información en estos excelentes diapositivas: http://mentat.za.net/numpy/numpy_advanced_slides/

con stride_tricks a partir de la diapositiva 29.

No estoy del todo claro sobre la cuestión, así que no puede sugerir algo más concreto - a pesar de que probablemente hacerlo en Cython o FORTRAN con f2py o con tejido. Cada vez me gustan más FORTRAN en este momento porque en el momento de agregar todas las anotaciones de tipos requeridos en Cython Creo que termina buscando menos claro que el FORTRAN.

Hay una comparación de estos enfoques aquí:

www. scipy. org / PerformancePython

(no se puede publicar más enlaces como yo soy un nuevo usuario) con un ejemplo que se parece a su caso.

Answer 8

Al final me encontré con el mismo problema que tú. Tuve que recurrir a utilizar Jacobi iteración y tejedora:

 while (iter_n < max_time_steps):
        expr = "field[1:-1, 1:-1] = (field[2:, 1:-1] "\
                                                      "+ field[:-2, 1:-1]+"\
                                                      "field[1:-1, 2:] +"\
                                                      "field[1:-1, :-2] )/4."                                       

        weave.blitz(expr, check_size=0)

         #Toroidal conditions
        field[:,0] = field[:,self.flow.n_x - 2]
        field[:,self.flow.n_x -1] = field[:,1]

        iter_n = iter_n + 1

Funciona y es rápido, pero no es de Gauss-Seidel, por lo que la convergencia puede ser un poco complicado. La única opción de hacer Gauss-Seidel como un bucle tradicional con índices.

Answer 9

Yo sugeriría Cython en lugar de bucle en c. hay fuerza haber alguna manera de lujo numpy de conseguir su ejemplo de trabajo utilizando una gran cantidad de pasos intermedios ... pero ya sabes cómo se escribe en C ya, acaba de escribir un poquito rápido como un Cython la función y dejar que la magia de Cython hacer el resto del trabajo fácil para usted.