Accelerare i calcoli con matrici intorpidite

Question

Ho due matrici. Entrambi sono pieni di zeri e uno. Uno è grande (3000 x 2000 elementi) e l'altro è più piccolo (20 x 20) elementi. Sto facendo qualcosa del tipo:

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

Che è dolorosamente lento. Sto facendo qualcosa di sbagliato? Esiste un modo intelligente per rendere questo lavoro più veloce?

modifica: Fondamentalmente lo sono, per ciascuno (x, y) nella matrice grande, controllando tutti i pixel sia della matrice grande che della matrice piccola intorno (x, y) per vedere se sono 1. Se sono 1 , quindi ho impostato quel valore su newMatrix. Sto facendo una sorta di rilevamento delle collisioni.

Answer 1

Posso pensare a un paio di ottimizzazioni lì - Mentre stai usando 4 python nidificati " per " affermazioni, sei più lento che puoi.

Non riesco a capire esattamente cosa stai cercando - ma per prima cosa, se la tua grande matrice "1" di densità è bassa, puoi certamente usare "qualsiasi" di Python funzione sulle sezioni di bigMtarix per verificare rapidamente se ci sono elementi impostati lì - potresti ottenere un aumento della velocità di diverse volte lì:

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
        (...) 

A questo punto, se hai ancora bisogno di interagire su ciascun elemento, fai un'altra coppia di indici per percorrere ogni posizione all'interno del passaggio - ma penso che tu abbia avuto l'idea.

Oltre all'utilizzo di operazioni numeriche interne come questa " qualsiasi " Durante l'utilizzo, potresti sicuramente aggiungere un codice di flusso di controllo per interrompere il ciclo (b, a) quando viene trovato il primo pixel corrispondente. (Ad esempio, inserendo un'istruzione "break" all'interno dell'ultima "if" e un'altra coppia if..break per il ciclo "b".

Non riesco davvero a capire esattamente quale sia il tuo intento, quindi non posso darti un codice più specifico.

Answer 2

Il tuo codice di esempio non ha senso, ma la descrizione del tuo problema sembra che tu stia provando a fare una convoluzione 2D di un piccolo bitarray sul grande bitarray. C'è una funzione convolve2d nel pacchetto scipy.signal che fa esattamente questo. Basta fare convolve2d (bigMatrix, smallMatrix) per ottenere il risultato. Sfortunatamente l'implementazione di scipy non ha un caso speciale per array booleani, quindi la convoluzione completa è piuttosto lenta. Ecco una funzione che sfrutta il fatto che gli array contengono solo uno e zero:

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

Sulla mia macchina gira in meno di 9 secondi per una convoluzione di 3000x2000 per 20x20. Il tempo di esecuzione dipende dal numero di quelli nell'array più piccolo, ovvero 20 ms per ogni elemento diverso da zero.

Answer 3

Se i tuoi bit sono davvero impacchettati 8 per byte / 32 per int, e puoi ridurre il tuo smallMatrix a 20x16,
quindi prova quanto segue, qui per una singola riga.
( newMatrix [x, y] = 1 quando qualsiasi bit del 20x16 attorno a x, y è 1 ?? Cosa stai davvero cercando?)

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