Acelerar cálculos com matrizes numpy

Question

Eu tenho duas matrizes. Ambos são preenchidos com zeros e uns. Um deles é um grande problema (3000 x 2000 elementos), e o outro é (20 x 20) elementos menores. Eu estou fazendo 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á sendo dolorosamente lento. Estou fazendo algo errado? Existe uma maneira inteligente para fazer este trabalho mais rápido?

edit: Basicamente, eu sou, para cada um (x, y) na grande matriz, verificando todos os pixels de ambos grande matriz ea pequena matriz em torno (x, y) para ver se eles são 1. Se eles são 1 , então eu definir esse valor em newMatrix. Estou fazendo uma espécie de detecção de colisão.

Answer 1

Não consigo pensar em um par de otimizações lá - Como você está usando 4 python nested "para" declarações, você é tão lento como você pode ser.

Eu não consigo descobrir exatamente o que você está procurando - mas para uma coisa, se a densidade seu grande matriz "1" s é baixa, certamente você pode usar "qualquer" função de python em fatias de bigMtarix para verificar rapidamente se existem elementos definido lá - você pode obter um aumento de velocidade várias vezes 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
        (...) 

Neste ponto, se ainda precisa interagir em cada elemento, você faz um outro par de índices de andar cada posição dentro da etapa - mas eu acho que você teve a idéia

.

Além de usar operações numéricas internas como esta "qualquer" uso, você certamente poderia adicionar algum código de controle de fluxo para quebrar-off a (b, a) loop quando o primeiro pixel correspondente for encontrado. (Como, inserindo uma declaração "break" dentro de sua última "se" e outro par if..break para o loop "b".

Eu realmente não pode descobrir exatamente o que a sua intenção é - então não posso dar-lhe um código mais specifc.

Answer 2

O código de exemplo não faz sentido, mas a descrição de seus sons problema como você está tentando fazer uma convolução 2D de uma pequena bitarray sobre o grande bitarray. Há uma função convolve2d no pacote scipy.signal que faz exatamente isso. Só não convolve2d(bigMatrix, smallMatrix) para obter o resultado. Infelizmente a implementação scipy não tem um caso especial para matrizes de booleanos para que a convolução completa é bastante lento. Aqui está uma função que tira proveito do fato de que as matrizes contêm somente uns e zeros:

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

Na minha máquina ele é executado em menos de 9 segundos para um 3000x2000 por 20x20 convolução. O tempo de execução depende do número de uns na matriz mais pequena, sendo 20 ms por cada elemento diferente de zero.

Answer 3

Se os seus bits são realmente embalado 8 por byte / 32 por int, e você pode reduzir o seu smallmatrix para 20x16,
em seguida, tente o seguinte, aqui por uma única linha.
(newMatrix[x, y] = 1 quando qualquer pouco da 20x16 torno x, y é 1 ?? O que você está realmente procurando?)

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