Accélération des calculs avec des matrices numpy

Question

J'ai deux matrices. Les deux sont remplis de zéros et de uns. L'une est une grande (3000 x 2000 éléments) et l'autre est plus petite (20 x 20). Je fais quelque chose comme:

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

Ce qui est terriblement lent. Est-ce que je fais quelque chose de mal? Existe-t-il un moyen intelligent d’accélérer le travail?

edit: Fondamentalement, je suis, pour chaque (x, y) dans la grande matrice, vérifiant tous les pixels de la grande matrice et de la petite matrice autour (x, y) pour voir s’ils sont 1. Si ils sont 1 , alors je mets cette valeur sur newMatrix. Je fais une sorte de détection de collision.

Answer 1

Je peux penser à quelques optimisations là-bas - Comme vous utilisez 4 python imbriqué " for " déclarations, vous êtes aussi lent que vous pouvez être.

Je ne peux pas comprendre exactement ce que vous recherchez - mais, d’une part, si la densité de votre grosse matrice "1" est faible, vous pouvez certainement utiliser "n'importe quel" de python. fonction sur les tranches de bigMtarix pour vérifier rapidement s’il existe des éléments définis - vous pouvez obtenir une augmentation de la vitesse multipliée par plusieurs:

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

À ce stade, si vous avez toujours besoin d'interagir sur chaque élément, vous créez une autre paire d'index pour parcourir chaque position dans l'étape, mais je pense que vous avez compris l'idée.

Hormis l'utilisation d'opérations numériques internes telles que celle-ci, "any". vous pouvez certainement ajouter du code de flux de contrôle pour rompre la boucle (b, a) lorsque le premier pixel correspondant est trouvé. (Par exemple, insérer une déclaration "break" dans votre dernière "if" et une autre paire if..break pour la boucle "b".

Je ne peux vraiment pas comprendre exactement quelle est votre intention - je ne peux donc pas vous donner de code plus spécifique.

Answer 2

Votre exemple de code n’a aucun sens, mais la description de votre problème donne l’impression que vous essayez de faire une convolution 2D d’un petit bitarray sur un gros bitarray. Le paquet scipy.signal contient une fonction convolve2d qui fait exactement cela. Il suffit de faire convolve2d (bigMatrix, smallMatrix) pour obtenir le résultat. Malheureusement, l'implémentation scipy n'a pas de cas particulier pour les tableaux booléens, donc la convolution complète est plutôt lente. Voici une fonction qui tire parti du fait que les tableaux ne contiennent que des uns et des zéros:

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

Sur ma machine, il s'exécute en moins de 9 secondes pour une convolution 3000x2000 x 20x20. La durée d'exécution dépend du nombre de uns dans le tableau le plus petit, soit 20 ms pour chaque élément différent de zéro.

Answer 3

Si vos bits sont vraiment chargés 8 octets / 32 par int, et vous pouvez réduire votre smallMatrix à 20x16,
puis essayez ce qui suit, ici pour une seule rangée.
( newMatrix [x, y] = 1 lorsque n’importe quel bit du 20x16 autour de x, y vaut 1 ?? Que cherchez-vous vraiment?)

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