Estendi la maschera numpy da n celle a destra per ogni valore errato, in modo efficiente

Question

Diciamo che ho un array di lunghezza 30 con 4 valori non validi. Voglio creare una maschera per quei valori sbagliati, ma poiché userò le funzioni della finestra mobile, mi piacerebbe anche un numero fisso di indici successivi dopo che ogni valore errato fosse contrassegnato come non valido. In seguito, n = 3:

mi piacerebbe fare questo nel modo più efficiente possibile, perché questa routine verrà eseguito molte volte su grandi serie di dati che contiene miliardi di datapoint. Quindi ho bisogno il più vicino possibile a una soluzione vettoriale stilizzata perché vorrei evitare i loop di Python.

A scanso di ribattitura, qui è la matrice:

import numpy as np 
a = np.array([4, 0, 8, 5, 10, 9, np.nan, 1, 4, 9, 9, np.nan, np.nan, 9,\ 
       9, 8, 0, 3, 7, 9, 2, 6, 7, 2, 9, 4, 1, 1, np.nan, 10]) 

Answer 1

Eppure un'altra risposta!
Basta prendere la maschera che hai già e applica le versioni logiche o spostate di se stesso. Ben vettorializzato e follemente veloce! : D

def repeat_or(a, n=4): 
    m = np.isnan(a) 
    k = m.copy() 

    # lenM and lenK say for each mask how many 
    # subsequent Trues there are at least 
    lenM, lenK = 1, 1 

    # we run until a combination of both masks will give us n or more 
    # subsequent Trues 
    while lenM+lenK < n: 
     # append what we have in k to the end of what we have in m 
     m[lenM:] |= k[:-lenM] 

     # swap so that m is again the small one 
     m, k = k, m 

     # update the lengths 
     lenM, lenK = lenK, lenM+lenK 

    # see how much m has to be shifted in order to append the missing Trues 
    k[n-lenM:] |= m[:-n+lenM] 

    return k 

Purtroppo non ho potuto ottenere m[i:] |= m[:-i] esecuzione ... probabilmente una cattiva idea sia per modificare e utilizzare la maschera di modificarsi. Funziona per m[:-i] |= m[i:], tuttavia questa è la direzione sbagliata.
In ogni caso, invece della crescita quadratica ora abbiamo una crescita simile a Fibonacci che è ancora migliore di quella lineare.
(ho mai pensato che avrei effettivamente scrivere un algoritmo che è in realtà legato alla sequenza di Fibonacci, senza essere qualche problema di matematica strano.)

test in condizioni "reali" con array di dimensioni e 1E6 Nans 1E5:

In [5]: a = np.random.random(size=1e6) 

In [6]: a[np.random.choice(np.arange(len(a), dtype=int), 1e5, replace=False)] = np.nan 

In [7]: %timeit reduceat(a) 
10 loops, best of 3: 65.2 ms per loop 

In [8]: %timeit index_expansion(a) 
100 loops, best of 3: 12 ms per loop 

In [9]: %timeit cumsum_trick(a) 
10 loops, best of 3: 17 ms per loop 

In [10]: %timeit repeat_or(a) 
1000 loops, best of 3: 1.9 ms per loop 

In [11]: %timeit agml_indexing(a) 
100 loops, best of 3: 6.91 ms per loop 

Lascerò ulteriori benchmark a Thomas.

Answer 2

qualcosa di simile?

maskleft = np.where(np.isnan(a))[0] 
maskright = maskleft + n 
mask = np.zeros(len(a),dtype=bool) 
for l,r in itertools.izip(maskleft,maskright): 
    mask[l:r] = True 

Oppure, dal momento che n è piccolo, potrebbe essere meglio per eseguire un ciclo su di esso invece:

maskleft = np.where(np.isnan(a))[0] 
mask = np.zeros(len(a),dtype=bool) 
for i in range(0,n): 
    thismask = maskleft+i 
    mask[thismask] = True 

Fatta eccezione per il ciclo su n, quanto sopra è completamente vettorizzati. Ma il ciclo è completamente parallelizzabile, quindi potresti essere in grado di ottenere un aumento del fattore n usando per es. multiprocessing o Cython, se sei disposto ad andare nei guai.

Modifica: per @askewchan la soluzione 2 può potenzialmente causare errori fuori intervallo. Ha anche problemi di indicizzazione nello range(0,n). possibilità di correzione:

maskleft = np.where(np.isnan(a))[0] 
ghost_mask = np.zeros(len(a)+n,dtype=bool) 
for i in range(0, n+1): 
    thismask = maskleft + i 
    ghost_mask[thismask] = True 
mask = ghost_mask[:len(ghost_mask)-n] 

Answer 3

È possibile utilizzare np.ufunc.reduceat con np.bitwise_or:

import numpy as np 
a = np.array([4, 0, 8, 5, 10, 9, np.nan, 1, 4, 9, 9, np.nan, np.nan, 9, 
       9, 8, 0, 3, 7, 9, 2, 6, 7, 2, 9, 4, 1, 1, np.nan, 10]) 
m = np.isnan(a) 
n = 4 
i = np.arange(1, len(m)+1) 
ind = np.column_stack([i-n, i]) # may be a faster way to generate this 
ind.clip(0, len(m)-1, out=ind) 

np.bitwise_or.reduceat(m, ind.ravel())[::2] 

Sul dati:

print np.column_stack([m, reduced]) 
[[False False] 
[False False] 
[False False] 
[False False] 
[False False] 
[False False] 
[ True True] 
[False True] 
[False True] 
[False True] 
[False False] 
[ True True] 
[ True True] 
[False True] 
[False True] 
[False True] 
[False False] 
[False False] 
[False False] 
[False False] 
[False False] 
[False False] 
[False False] 
[False False] 
[False False] 
[False False] 
[False False] 
[False False] 
[ True True] 
[False True]] 

Answer 4

Questo potrebbe anche essere considerato un problema morphological dilation, usando qui il scipy.ndimage.binary_dilation:

def dilation(a, n): 
    m = np.isnan(a) 
    s = np.full(n, True, bool) 
    return ndimage.binary_dilation(m, structure=s, origin=-(n//2)) 

Nota sulla origin: questo argomento assicura la structure (lo chiamerei un kernel) inizia un po 'a sinistra del input (la maschera m). Normalmente, il valore a out[i] corrisponde alla dilatazione con il centro di structure (che sarebbe structure[n//2]) a in[i], ma si desidera che lo sia allo in[i].

Si può anche fare questo con un kernel che è imbottita sulla sinistra con False s, che è quello che sarebbe necessario se si è utilizzato il binary_dilation da scikit-image:

def dilation_skimage(a, n): 
    m = np.isnan(a) 
    s = np.zeros(2*n - n%2, bool) 
    s[-n:] = True 
    return skimage.morphology.binary_dilation(m, selem=s) 

Timing non sembra cambiare troppo tra i due:

dilation_scipy 
small: 10 loops, best of 3: 47.9 ms per loop 
large: 10000 loops, best of 3: 88.9 µs per loop 

dilation_skimage 
small: 10 loops, best of 3: 47.0 ms per loop 
large: 10000 loops, best of 3: 91.1 µs per loop 

Answer 5

È possibile utilizzare lo stesso trucco cumSum come si farebbe per un filtro medio in esecuzione:

def cumsum_trick(a, n): 
    mask = np.isnan(a) 
    cs = np.cumsum(mask) 
    cs[n:] -= cs[:-n].copy() 
    return cs > 0 

Sfortunatamente l'ulteriore .copy() è necessario, a causa di un po 'di buffering che va avanti internamente a dall'ordine delle operazioni. E 'possibile convincere NumPy per applicare la sottrazione in senso inverso, ma per quella di lavorare la matrice cs deve avere un passo negativo:

def cumsum_trick_nocopy(a, n): 
    mask = np.isnan(a) 
    cs = np.cumsum(mask, out=np.empty_like(a, int)[::-1]) 
    cs[n:] -= cs[:-n] 
    out = cs > 0 
    return out 

ma questo sembra fragile e non è poi così molto più veloce in ogni caso.

mi chiedo se c'è una compilato funzione di elaborazione del segnale da qualche parte che fa questa operazione esatto ..

---

Per le maschere iniziali sparse e piccole n questo è anche abbastanza veloce:

def index_expansion(a, n): 
    mask = np.isnan(a) 
    idx = np.flatnonzero(mask) 
    expanded_idx = idx[:,None] + np.arange(1, n) 
    np.put(mask, expanded_idx, True, 'clip') 
    return mask 

Answer 6

OP qui con i risultati del benchmark. Ho incluso il mio ("op") con il quale avevo iniziato, che scorre sugli indici non validi e aggiunge 1 ... n a loro quindi accetta gli uniques per trovare gli indici delle maschere. Puoi vederlo nel codice qui sotto con tutte le altre risposte.

In ogni caso, ecco i risultati. Le sfaccettature sono la dimensione della matrice lungo x (10 thru 10e7) e la dimensione della finestra lungo y (5, 50, 500, 5000). Quindi è coder in ogni aspetto, con un punteggio di log-10 perché stiamo parlando di microsecondi attraverso i minuti.

@swenzel sembra essere il vincitore con la sua seconda risposta, spostando prima risposta di @ moarningsun (seconda risposta di moarningsun stava crollando la macchina attraverso l'uso di memoria di massa, ma questo è probabilmente perché non è stato progettato per n grande o non sparse a).

Il grafico non rende giustizia al più veloce di questi contributi a causa della scala del registro (necessaria). Sono dozzine, centinaia di volte più veloci delle soluzioni di looping decenti. swenzel1 è 1000 volte più veloce di op nel caso più grande, e op sta già facendo uso di numpy.

Si noti che ho utilizzato una versione numpy compilata rispetto alle librerie Intel MKL ottimizzate che fanno pieno uso delle istruzioni AVX presenti dal 2012. In alcuni casi di utilizzo di vettori questo aumenterà la velocità di i7/Xeon di un fattore di 5 Alcuni dei contributi potrebbero giovare più di altri.

Ecco il codice completo per eseguire tutte le risposte inviate finora, inclusa la mia. La funzione allagree() assicura che i risultati siano corretti, mentre timeall() ti darà un Dataframe di panda di lunga durata con tutti i risultati in secondi.

È possibile eseguirlo abbastanza facilmente con un nuovo codice o modificare le mie ipotesi. Si prega di tenere presente che non ho preso in considerazione altri fattori come l'utilizzo della memoria. Inoltre, ho fatto ricorso a R ggplot2 per la grafica in quanto non conosco Seaborn/Matplotlib abbastanza bene da farlo fare ciò che voglio.

Per completezza, tutti i risultati sono d'accordo:

In [4]: allagree(n = 7, asize = 777) 
Out[4]: 
      AGML0 AGML1 askewchan0 askewchan1 askewchan2 moarningsun0 \ 
AGML0   True True  True  True  True   True 
AGML1   True True  True  True  True   True 
askewchan0 True True  True  True  True   True 
askewchan1 True True  True  True  True   True 
askewchan2 True True  True  True  True   True 
moarningsun0 True True  True  True  True   True 
swenzel0  True True  True  True  True   True 
swenzel1  True True  True  True  True   True 
op   True True  True  True  True   True 

      swenzel0 swenzel1 op 
AGML0   True  True True 
AGML1   True  True True 
askewchan0  True  True True 
askewchan1  True  True True 
askewchan2  True  True True 
moarningsun0  True  True True 
swenzel0   True  True True 
swenzel1   True  True True 
op    True  True True 

Grazie a tutti coloro che ha presentato!

Codice in materia di grafica dopo l'esportazione di uscita del timeall() utilizzando pd.to_csv e read.csv in R:

ww <- read.csv("ww.csv")  
ggplot(ww, aes(x=coder, y=value, col = coder)) + geom_point(size = 3) + scale_y_continuous(trans="log10")+ facet_grid(nsize ~ asize) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + ggtitle("Fastest by coder") + ylab("time (seconds)") 

codice per il test:

# test Stack Overflow 32706135 nan shift routines 

import numpy as np 
import pandas as pd 
import seaborn as sns 
import matplotlib.pyplot as plt 
from timeit import Timer 
from scipy import ndimage 
from skimage import morphology 
import itertools 
import pdb 
np.random.seed(8472) 


def AGML0(a, n):        # loop itertools 
    maskleft = np.where(np.isnan(a))[0] 
    maskright = maskleft + n 
    mask = np.zeros(len(a),dtype=bool) 
    for l,r in itertools.izip(maskleft,maskright): 
     mask[l:r] = True 
    return mask 


def AGML1(a, n):        # loop n 
    nn = n - 1 
    maskleft = np.where(np.isnan(a))[0] 
    ghost_mask = np.zeros(len(a)+nn,dtype=bool) 
    for i in range(0, nn+1): 
     thismask = maskleft + i 
     ghost_mask[thismask] = True 
    mask = ghost_mask[:len(ghost_mask)-nn] 
    return mask 


def askewchan0(a, n): 
    m = np.isnan(a) 
    i = np.arange(1, len(m)+1) 
    ind = np.column_stack([i-n, i]) # may be a faster way to generate this 
    ind.clip(0, len(m)-1, out=ind) 
    return np.bitwise_or.reduceat(m, ind.ravel())[::2] 


def askewchan1(a, n): 
    m = np.isnan(a) 
    s = np.full(n, True, bool) 
    return ndimage.binary_dilation(m, structure=s, origin=-(n//2)) 


def askewchan2(a, n): 
    m = np.isnan(a) 
    s = np.zeros(2*n - n%2, bool) 
    s[-n:] = True 
    return morphology.binary_dilation(m, selem=s) 


def moarningsun0(a, n): 
    mask = np.isnan(a) 
    cs = np.cumsum(mask) 
    cs[n:] -= cs[:-n].copy() 
    return cs > 0 


def moarningsun1(a, n): 
    mask = np.isnan(a) 
    idx = np.flatnonzero(mask) 
    expanded_idx = idx[:,None] + np.arange(1, n) 
    np.put(mask, expanded_idx, True, 'clip') 
    return mask 


def swenzel0(a, n): 
    m = np.isnan(a) 
    k = m.copy() 
    for i in range(1, n): 
     k[i:] |= m[:-i] 
    return k 


def swenzel1(a, n=4): 
    m = np.isnan(a) 
    k = m.copy() 

    # lenM and lenK say for each mask how many 
    # subsequent Trues there are at least 
    lenM, lenK = 1, 1 

    # we run until a combination of both masks will give us n or more 
    # subsequent Trues 
    while lenM+lenK < n: 
     # append what we have in k to the end of what we have in m 
     m[lenM:] |= k[:-lenM] 

     # swap so that m is again the small one 
     m, k = k, m 

     # update the lengths 
     lenM, lenK = lenK, lenM+lenK 

    # see how much m has to be shifted in order to append the missing Trues 
    k[n-lenM:] |= m[:-n+lenM] 
    return k 


def op(a, n): 
    m = np.isnan(a) 
    for x in range(1, n): 
     m = np.logical_or(m, np.r_[False, m][:-1]) 
    return m 


# all the functions in a list. NB these are the actual functions, not their names 
funcs = [AGML0, AGML1, askewchan0, askewchan1, askewchan2, moarningsun0, swenzel0, swenzel1, op] 

def allagree(fns = funcs, n = 10, asize = 100): 
    """ make sure result is the same from all functions """ 
    fnames = [f.__name__ for f in fns] 
    a = np.random.rand(asize) 
    a[np.random.randint(0, asize, int(asize/10))] = np.nan 
    results = dict([(f.__name__, f(a, n)) for f in fns]) 
    isgood = [[np.array_equal(results[f1], results[f2]) for f1 in fnames] for f2 in fnames] 
    pdgood = pd.DataFrame(isgood, columns = fnames, index = fnames) 
    if not all([all(x) for x in isgood]): 
     print "not all results identical" 
     pdb.set_trace() 
    return pdgood 


def timeone(f): 
    """ time one of the functions across the full range of a nd n """ 
    print "Timing", f.__name__ 
    Ns = np.array([10**x for x in range(0, 4)]) * 5 # 5 to 5000 window size 
    As = [np.random.rand(10 ** x) for x in range(1, 8)] # up to 10 million data data points 
    for i in range(len(As)): # 10% of points are always bad 
     As[i][np.random.randint(0, len(As[i]), len(As[i])/10)] = np.nan 
    results = np.array([[Timer(lambda: f(a, n)).timeit(number = 1) if n < len(a) \ 
         else np.nan for n in Ns] for a in As]) 
    pdresults = pd.DataFrame(results, index = [len(x) for x in As], columns = Ns) 
    return pdresults 


def timeall(fns = funcs): 
    """ run timeone for all known funcs """ 
    testd = dict([(x.__name__, timeone(x)) for x in fns]) 
    testdf = pd.concat(testd.values(), axis = 0, keys = testd.keys()) 
    testdf.index.names = ["coder", "asize"] 
    testdf.columns.names = ["nsize"] 
    testdf.reset_index(inplace = True) 
    testdf = pd.melt(testdf, id_vars = ["coder", "asize"]) 
    return testdf