Acelerar los cálculos con matrices numpy

Question

Tengo dos matrices. Ambos están llenos de ceros y unos. Uno es uno grande (3000 x 2000 elementos) y el otro es más pequeño (20 x 20). Estoy haciendo algo como:

newMatrix = (size of bigMatrix), filled with zeros
l = (a constant)

for y in xrange(0, len(bigMatrix[0])):
    for x in xrange(0, len(bigMatrix)):

        for b in xrange(0, len(smallMatrix[0])):
            for a in xrange(0, len(smallMatrix)):

                if (bigMatrix[x, y] == smallMatrix[x + a - l, y + b - l]):
                    newMatrix[x, y] = 1

Que está siendo dolorosamente lento. ¿Estoy haciendo algo mal? ¿Hay alguna forma inteligente de hacer que esto funcione más rápido?

edit: Básicamente, para cada (x, y) en la matriz grande, verifico todos los píxeles de la matriz grande y la matriz pequeña alrededor de (x, y) para ver si son 1. Si son 1 , luego establezco ese valor en newMatrix. Estoy haciendo una especie de detección de colisión.

Answer 1

Puedo pensar en un par de optimizaciones allí - Como estás usando 4 python anidados " para " afirmaciones, eres lo más lento que puedes ser.

No puedo entender exactamente lo que estás buscando - pero, por un lado, si la densidad de su gran matriz de 1 es baja, puede usar Python's cualquier Funciona en las rodajas de bigMtarix para comprobar rápidamente si hay elementos establecidos allí. Podrías obtener un aumento de la velocidad varias veces:

step = len(smallMatrix[0])
for y in xrange(0, len(bigMatrix[0], step)):
    for x in xrange(0, len(bigMatrix), step):
        if not any(bigMatrix[x: x+step, y: y + step]):
            continue
        (...) 

En este punto, si aún necesitas interactuar con cada elemento, haces otro par de índices para recorrer cada posición dentro del paso, pero creo que tienes la idea.

Además de usar operaciones numéricas internas como esta " cualquier " sin duda, puede agregar algún código de flujo de control para interrumpir el ciclo (b, a) cuando se encuentre el primer píxel coincidente. (Por ejemplo, insertando una declaración de "ruptura" dentro de su último " if " y otra pareja if..break para el bucle " b "

Realmente no puedo entender exactamente cuál es tu intención, así que no puedo darte un código más específico.

Answer 2

Su código de ejemplo no tiene sentido, pero la descripción de su problema parece que está tratando de hacer una convolución 2d de un pequeño bitarray sobre el gran bitarray. Hay una función convolve2d en el paquete scipy.signal que hace exactamente esto. Simplemente haga convolve2d (bigMatrix, smallMatrix) para obtener el resultado. Desafortunadamente, la implementación de scipy no tiene un caso especial para arreglos booleanos, por lo que la convolución completa es bastante lenta. Aquí hay una función que aprovecha el hecho de que los arreglos contienen solo unos y ceros:

import numpy as np

def sparse_convolve_of_bools(a, b):
    if a.size < b.size:
        a, b = b, a
    offsets = zip(*np.nonzero(b))
    n = len(offsets)
    dtype = np.byte if n < 128 else np.short if n < 32768 else np.int
    result = np.zeros(np.array(a.shape) + b.shape - (1,1), dtype=dtype)
    for o in offsets:
        result[o[0]:o[0] + a.shape[0], o[1]:o[1] + a.shape[1]] += a
    return result

En mi máquina, se ejecuta en menos de 9 segundos para una convolución de 3000x2000 por 20x20. El tiempo de ejecución depende de la cantidad de unidades en la matriz más pequeña, que son 20 ms por cada elemento distinto de cero.

Answer 3

Si sus bits están realmente empaquetados 8 por byte / 32 por int, y puedes reducir tu smallMatrix a 20x16,
a continuación, intente lo siguiente, aquí para una sola fila.
( newMatrix [x, y] = 1 cuando cualquier bit del 20x16 alrededor de x, y es 1 ?? ¿Qué estás buscando realmente?)

python -m timeit -s '
""" slide 16-bit mask across 32-bit pairs bits[j], bits[j+1] """

import numpy as np

bits = np.zeros( 2000 // 16, np.uint16 )  # 2000 bits
bits[::8] = 1
mask = 32+16
nhit = 16 * [0]

def hit16( bits, mask, nhit ):
    """
        slide 16-bit mask across 32-bit pairs bits[j], bits[j+1]
        bits: long np.array( uint16 )
        mask: 16 bits, int
        out: nhit[j] += 1 where pair & mask != 0
    """
    left = bits[0]
    for b in bits[1:]:
        pair = (left << 16) | b
        if pair:  # np idiom for non-0 words ?
            m = mask
            for j in range(16):
                if pair & m:
                    nhit[j] += 1
                    # hitposition = jb*16 + j
                m <<= 1
        left = b
    # if any(nhit):  print "hit16:", nhit

' \
'
hit16( bits, mask, nhit )
'

# 15 msec per loop, bits[::4] = 1
# 11 msec per loop, bits[::8] = 1
# mac g4 ppc