2014-06-05 16 views
6

Si consideri il seguente codice usando array NumPy che è molto lento:Accelerare un ciclo di Numpy in 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 

come riscrivere la funzione di velocizzare il calcolo? (np.size(tx) == 10000 e np.size(ox) == 100000)

+0

Considerate anche l'utilizzo di OpenCL? –

+0

Non ho bisogno di prestazioni complete, voglio solo una velocità raw. – Vincent

+1

Costruisci un 'scipy.spatial.KDTree' dai punti tx, ty, tz e poi usa la ricerca del vicino più prossimo nella norma dell'infinito per ogni punto in bue, oy, oz per vedere se c'è qualche punto abbastanza vicino. –

risposta

6

si assegnano 10000 elenchi delle dimensioni 100000. La prima cosa da fare sarebbe quella di smettere di usare range per la nidificato j loop e utilizzare la versione generatore xrange invece. Questo ti farà risparmiare tempo e spazio per allocare tutte queste liste.

Il prossimo potrebbe essere quella di utilizzare le operazioni vectorized:

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

Cambio fatto, ma non fornisce un'enorme accelerazione – Vincent

+0

Dammi un secondo, c'è il prossimo in arrivo :) –

+0

Ok, il è in forma vettoriale (potrebbe essere necessario verificare le condizioni), ma come hai modificato la tua domanda in risposta alla mia risposta originale che parte della risposta è diventata irrilevante, anche se sarebbe ancora un buon passo iniziale da fare in caso di dati non accidentati :) –

0

credo che questo vi darà lo stesso risultato per il doppio loop e essere più veloce:

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]) 

per ottenere: 1) riordinare i loop (vale a dire, di swap t orlo) che è valido poiché nulla cambia all'interno dei loop; 2) tirare result[j] al di fuori del ciclo i; 3) convertire tutti gli t>ox-oe and t<ox+oe in abs(t-ox)<oe (anche se questo potrebbe non essere un enorme aumento della velocità, è più facile da leggere).

Poiché non si dispone di codice eseguibile e non volevo creare un test per questo, non sono sicuro al 100% che sia corretto.

Problemi correlati