Python/NumPy: attuare una somma parziale (ma non abbastanza)

Question

Dato sono due matrici di uguale lunghezza, uno dati Holding, uno che tengono i risultati ma inizialmente impostato a zero, per esempio:

a = numpy.array([1, 0, 0, 1, 0, 1, 0, 0, 1, 1]) 
b = numpy.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0]) 

avevo piace calcolare la somma di tutti i possibili sottoinsiemi di tre elementi adiacenti in a. Se la somma è 0 o 1, i tre elementi corrispondenti in b rimangono invariati; solo se la somma è superiore a 1 sono i tre elementi corrispondenti b impostato a 1, in modo che dopo il calcolo b diventa

array([0, 0, 0, 1, 1, 1, 0, 1, 1, 1]) 

Un semplice ciclo raggiungerà questo:

for x in range(len(a)-2): 
    if a[x:x+3].sum() > 1: 
     b[x:x+3] = 1 

Dopo questo, b ha la forma desiderata

Devo farlo per una grande quantità di dati, quindi la velocità è un problema. C'è un modo più veloce in NumPy per eseguire l'operazione di cui sopra?

(Capisco che questo è simile a una convoluzione, ma non proprio la stessa).

Answer 1

si può iniziare con una convoluzione, scegliere i valori che superano 1, e, infine, utilizza una "dilatazione":

b = numpy.convolve(a, [1, 1, 1], mode="same") > 1 
b = b | numpy.r_[0, b[:-1]] | numpy.r_[b[1:], 0] 

Dal momento che questo evita il ciclo Python, dovrebbe essere più veloce del tuo approccio, ma io non ha fatto i tempi.

Un'alternativa è quella di utilizzare una seconda circonvoluzione per dilatare:

kernel = [1, 1, 1] 
b = numpy.convolve(a, kernel, mode="same") > 1 
b = numpy.convolve(b, kernel, mode="same") > 0 

Se avete SciPy a disposizione, ancora un'altra opzione per la dilatazione è

b = numpy.convolve(a, [1, 1, 1], mode="same") > 1 
b = scipy.ndimage.morphology.binary_dilation(b) 

Edit: Facendo some timings, Ho trovato che questa soluzione sembra essere la più veloce per i grandi array:

b = numpy.convolve(a, kernel) > 1 
b[:-1] |= b[1:] # Shift and "smearing" to the *left* (smearing with b[1:] |= b[:-1] does not work) 
b[:-1] |= b[1:] # … and again! 
b = b[:-2] 

Per un array di un milione di voci, era più di 200 volte più veloce del tuo approccio originale sulla mia macchina. Come sottolineato da EOL nei commenti, questa soluzione potrebbe essere considerata un po 'fragile, tuttavia, poiché dipende dai dettagli di implementazione di NumPy.

Answer 2

È possibile calcolare le somme "convoluzione" in modo efficiente con:

>>> a0 = a[:-2] 
>>> a1 = a[1:-1] 
>>> a2 = a[2:] 
>>> a_large_sum = a0 + a1 + a2 > 1 

Aggiornamento b può allora essere fatto in modo efficiente scrivendo qualcosa che significa "almeno uno dei tre vicini a_large_sum valori è vera" : per la prima volta si estende a_large_sum gamma di nuovo allo stesso numero di elementi come a (a destra, a sinistra ea destra, e poi a sinistra):

>>> a_large_sum_0 = np.hstack([a_large_sum, [False, False]]) 
>>> a_large_sum_1 = np.hstack([[False], a_large_sum, [False]]) 
>>> a_large_sum_2 = np.hstack([[False, False], a_large_sum]) 

È quindi ottiene b in modo efficiente:

>>> b = a_large_sum_0 | a_large_sum_1 | a_large_sum_2 

Questo dà il risultato che si ottiene, ma in modo molto efficiente, attraverso una leva di NumPy cicli veloci interne.

PS: Questo approccio è essenzialmente lo stesso della prima soluzione di Sven, ma è molto più pedonale dell'elegante codice di Sven; è altrettanto veloce, tuttavia. La seconda soluzione di Sven (doppia convolve()) è ancora più elegante, e è due volte più veloce.

Answer 3

Ti potrebbe piacere anche dare un'occhiata a NumPy's stride_tricks. Utilizzando impostazione tempi di Sven (vedi link in risposta di Sven), ho scoperto che per il (molto) grandi array, questo è anche un modo veloce per fare quello che vuoi (vale a dire con la tua definizione di a):

shape = (len(a)-2,3) 
strides = a.strides+a.strides 
a_strided = numpy.lib.stride_tricks.as_strided(a, shape=shape, strides=strides) 
b = np.r_[numpy.sum(a_strided, axis=-1) > 1, False, False] 
b[2:] |= b[1:-1] | b[:-2] 

Dopo la modifica (vedi i commenti sotto) non è più il modo più veloce.

Questo crea una vista particolarmente graduale sull'array originale. I dati in a non vengono copiati, ma vengono semplicemente visualizzati in un nuovo modo. Vogliamo fondamentalmente creare un nuovo array in cui l'ultimo indice contiene i sotto-array che vogliamo sommare (cioè i tre elementi che volete sommare). In questo modo, possiamo facilmente sommare alla fine con l'ultimo comando.

L'ultimo elemento di questa nuova forma deve pertanto essere 3, e il primo elemento sarà la lunghezza della vecchia a meno 2 (perché possiamo riassumere solo fino al -2 nd elemento).

L'elenco delle espansioni contiene le falcate, in byte, che il nuovo array a_strided deve eseguire per raggiungere l'elemento successivo in ciascuna dimensione della forma. Se si impostano questi uguali, significa che a_strided[0,1] e a_strided[1,0] saranno entrambi a[1], che è esattamente ciò che vogliamo. In un array normale questo non sarebbe il caso (il primo passo sarebbe "dimensione della prima dimensione della lunghezza della matrice-prima dimensione (= forma [0])"), ma in questo caso possiamo farne buon uso.

Non sono sicuro se ho spiegato tutto questo molto bene, ma stampo appena a_strided e vedrai quale sia il risultato e quanto sia facile questa operazione.