lunghezza scoperta di sequenze di valori identici in una matrice NumPy (run-length encoding)

Question

In un programma pylab (che potrebbe probabilmente essere un programma MATLAB pure) Ho una matrice NumPy di ​​numeri che rappresentano distanze: d[t] rappresenta l'distanza al momento t (e il periodo di tempo dei miei dati è len(d) unità di tempo).

Gli eventi a cui sono interessato sono quando la distanza è inferiore a una determinata soglia e voglio calcolare la durata di questi eventi. È facile ottenere una serie di booleani con b = d<threshold e il problema si riduce al calcolo della sequenza delle lunghezze delle parole True-only in b. Ma non so come farlo in modo efficiente (cioè usando le primitive numpy), e ho fatto ricorso a walk the array e per fare il rilevamento del cambio manuale (cioè inizializzare contatore quando il valore passa da False a True, aumentare il contatore fintanto che il valore è True e invia il contatore alla sequenza quando il valore torna a False). Ma questo è tremendamente lento.

Come rilevare in modo efficiente quel tipo di sequenze in array numpy?

Di seguito è riportato un codice python che illustra il mio problema: il quarto punto prende un tempo molto lungo per apparire (se non, aumentare la dimensione della matrice)

from pylab import * 

threshold = 7 

print '.' 
d = 10*rand(10000000) 

print '.' 

b = d<threshold 

print '.' 

durations=[] 
for i in xrange(len(b)): 
    if b[i] and (i==0 or not b[i-1]): 
     counter=1 
    if i>0 and b[i-1] and b[i]: 
     counter+=1 
    if (b[i-1] and not b[i]) or i==len(b)-1: 
     durations.append(counter) 

print '.' 

Answer 1

Anche se non numpy primitive , itertools funzioni sono spesso molto veloce, in modo da fare dare questa una prova (e misurare volte per varie soluzioni, tra cui questo, ovviamente):

def runs_of_ones(bits): 
    for bit, group in itertools.groupby(bits): 
    if bit: yield sum(group) 

Se sono necessari i valori in un elenco, è possibile utilizzare solo l'elenco (runs_of_ones (bit)), naturalmente; ma forse un elenco di comprensione potrebbe essere marginalmente ancora più veloce:

def runs_of_ones_list(bits): 
    return [sum(g) for b, g in itertools.groupby(bits) if b] 

Trasferirsi in possibilità "NumPy-native", che dire:

def runs_of_ones_array(bits): 
    # make sure all runs of ones are well-bounded 
    bounded = numpy.hstack(([0], bits, [0])) 
    # get 1 at run starts and -1 at run ends 
    difs = numpy.diff(bounded) 
    run_starts, = numpy.where(difs > 0) 
    run_ends, = numpy.where(difs < 0) 
    return run_ends - run_starts 

Anche in questo caso: essere sicuri di soluzioni di riferimento contro ogni altri in realistic- per-voi esempi!

Answer 2

durations = [] 
counter = 0 

for bool in b: 
    if bool: 
     counter += 1 
    elif counter > 0: 
     durations.append(counter) 
     counter = 0 

if counter > 0: 
    durations.append(counter) 

Answer 3

Solo nel caso qualcuno è curioso (e dal momento che lei ha citato MATLAB di passaggio), ecco un modo per risolverlo in MATLAB:

threshold = 7; 
d = 10*rand(1,100000); % Sample data 
b = diff([false (d < threshold) false]); 
durations = find(b == -1)-find(b == 1); 

Io non sono troppo familiarità con Python, ma forse questo potrebbe aiutare darti qualche idea. =)

Answer 4

Ecco una soluzione che utilizza solo matrici: accetta un array contenente una sequenza di bool e conta la lunghezza delle transizioni.

>>> from numpy import array, arange 
>>> b = array([0,0,0,1,1,1,0,0,0,1,1,1,1,0,0], dtype=bool) 
>>> sw = (b[:-1]^b[1:]); print sw 
[False False True False False True False False True False False False 
    True False] 
>>> isw = arange(len(sw))[sw]; print isw 
[ 2 5 8 12] 
>>> lens = isw[1::2] - isw[::2]; print lens 
[3 4] 

sw contiene una vera dove c'è un interruttore, isw li converte in indici. Gli articoli di isw vengono quindi sottratti a coppie in lens.

Si noti che se la sequenza inizia con un 1, si conterebbe la lunghezza delle sequenze dello 0: questo può essere risolto nell'indicizzazione per calcolare l'obiettivo. Inoltre, non ho testato casi d'angolo tali sequenze di lunghezza 1.

---

funzione completa che restituisce posizioni di partenza e le lunghezze di tutte le True -subarrays.

import numpy as np 

def count_adjacent_true(arr): 
    assert len(arr.shape) == 1 
    assert arr.dtype == np.bool 
    if arr.size == 0: 
     return np.empty(0, dtype=int), np.empty(0, dtype=int) 
    sw = np.insert(arr[1:]^arr[:-1], [0, arr.shape[0]-1], values=True) 
    swi = np.arange(sw.shape[0])[sw] 
    offset = 0 if arr[0] else 1 
    lengths = swi[offset+1::2] - swi[offset:-1:2] 
    return swi[offset:-1:2], lengths 

testato per diverse bool 1D-array (array vuoto; elementi singoli/multipli; pari/dispari lunghezze; iniziato con True/False; con solo True/False elementi).

Answer 5

RLE completamente numerico e vettoriale generico per qualsiasi array (funziona anche con stringhe, booleani ecc.).

Emette tuple di lunghezze di corsa, posizioni iniziali e valori.

import numpy as np 

def rle(inarray): 
     """ run length encoding. Partial credit to R rle function. 
      Multi datatype arrays catered for including non Numpy 
      returns: tuple (runlengths, startpositions, values) """ 
     ia = np.asarray(inarray)     # force numpy 
     n = len(ia) 
     if n == 0: 
      return (None, None, None) 
     else: 
      y = np.array(ia[1:] != ia[:-1])  # pairwise unequal (string safe) 
      i = np.append(np.where(y), n - 1) # must include last element posi 
      z = np.diff(np.append(-1, i))  # run lengths 
      p = np.cumsum(np.append(0, z))[:-1] # positions 
      return(z, p, ia[i]) 

Abbastanza veloce (i7):

xx = np.random.randint(0, 5, 1000000) 
%timeit yy = rle(xx) 
100 loops, best of 3: 18.6 ms per loop 

più tipi di dati:

rle([True, True, True, False, True, False, False]) 
Out[8]: 
(array([3, 1, 1, 2]), 
array([0, 3, 4, 5]), 
array([ True, False, True, False], dtype=bool)) 

rle(np.array([5, 4, 4, 4, 4, 0, 0])) 
Out[9]: (array([1, 4, 2]), array([0, 1, 5]), array([5, 4, 0])) 

rle(["hello", "hello", "my", "friend", "okay", "okay", "bye"]) 
Out[10]: 
(array([2, 1, 1, 2, 1]), 
array([0, 2, 3, 4, 6]), 
array(['hello', 'my', 'friend', 'okay', 'bye'], 
     dtype='|S6')) 

stessi risultati come Alex Martelli di cui sopra:

xx = np.random.randint(0, 2, 20) 

xx 
Out[60]: array([1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1]) 

am = runs_of_ones_array(xx) 

tb = rle(xx) 

am 
Out[63]: array([4, 5, 2, 5]) 

tb[0][tb[2] == 1] 
Out[64]: array([4, 5, 2, 5]) 

%timeit runs_of_ones_array(xx) 
10000 loops, best of 3: 28.5 µs per loop 

%timeit rle(xx) 
10000 loops, best of 3: 38.2 µs per loop 

Leggermente più lento di Alex (ma ancora molto veloce) e molto più flessibile.