2014-06-05 3 views
6

Consultez le code suivant à l'aide des tableaux numpy qui est très lent:Accélérer une boucle numpy en python?

# Intersection of an octree and a trajectory 
def intersection(octree, trajectory): 
    # Initialize numpy arrays 
    ox = octree.get("x") 
    oy = octree.get("y") 
    oz = octree.get("z") 
    oe = octree.get("extent")/2 
    tx = trajectory.get("x") 
    ty = trajectory.get("y") 
    tz = trajectory.get("z") 
    result = np.zeros(np.size(ox)) 
    # Loop over elements 
    for i in range(0, np.size(tx)): 
     for j in range(0, np.size(ox)): 
      if (tx[i] > ox[j]-oe[j] and 
       tx[i] < ox[j]+oe[j] and 
       ty[i] > oy[j]-oe[j] and 
       ty[i] < oy[j]+oe[j] and 
       tz[i] > oz[j]-oe[j] and 
       tz[i] < oz[j]+oe[j]): 
       result[j] += 1 
    # Finalize 
    return result 

Comment réécrire la fonction pour accélérer le calcul? (np.size(tx) == 10000 et np.size(ox) == 100000)

+0

Envisagez-vous également d'utiliser OpenCL? –

+0

Je n'ai pas besoin de performances complètes, je veux juste une accélération brute. – Vincent

+1

Construire un 'scipy.spatial.KDTree' à partir des points tx, ty, tz et ensuite utiliser la recherche du plus proche voisin dans la norme infinity pour chaque point dans ox, oy, oz pour voir s'il y a un point assez proche. –

Répondre

6

Vous allouent 10000 listes de taille 100000. La première chose à faire serait de cesser d'utiliser range pour la boucle j imbriquée et utiliser la version du générateur xrange à la place. Cela vous permettra d'économiser du temps et de l'espace allouer toutes ces listes.

Le prochain serait d'utiliser les opérations vectorisés:

for i in xrange(0, np.size(tx)): 
    index = (ox-oe < tx[i]) & (ox+oe > tx[i]) & (oy-oe < ty[i]) & (oy+oe > ty[i]) & (oz-oe < tz[i]) & (oz+oe > tz[i]) 
    result[index] += 1 
+0

Changement fait, mais il ne fournit pas une énorme accélération – Vincent

+0

Donnez-moi une seconde, il y a le prochain à venir :) –

+0

Ok, le La forme vectorisée est en (peut nécessiter une vérification des conditions), mais comme vous avez modifié votre question en réponse à ma réponse initiale, une partie de la réponse est devenue inutile même si ce serait encore une bonne première étape en cas de données non cahoteuses. :) –

0

Je pense que cela donnera le même résultat pour la double boucle et être plus rapide:

for j in xrange(np.size(ox)): 
    result[j] += sum(abs(tx-ox[j])<oe[j] & abs(ty-oy[j])<oe[j] & abs(tz-oz[j])<oe[j]) 

pour obtenir: 1) les boucles réordonner (c.-à-échange t hem) qui est valide puisque rien ne change dans les boucles; 2) tirez result[j] en dehors de la boucle i; 3) convertir tous les t>ox-oe and t<ox+oe en abs(t-ox)<oe (même si cela ne peut pas être une énorme accélération, il est plus facile à lire).

Comme vous n'avez pas de code exécutable, et que je n'ai pas voulu faire de test pour cela, je ne suis pas sûr à 100% que c'est correct.

Questions connexes