python - combinazione di argsort con mascheramento per ottenere valori più vicini nella finestra mobile

Question

Ho un codice per calcolare i valori mancanti in un'immagine, in base ai valori adiacenti in una finestra circolare 2D. Utilizza anche i valori di una o più immagini adiacenti temporaneamente nelle stesse posizioni (vale a dire la stessa finestra 2D spostata nella terza dimensione).

Per ogni posizione mancante, devo calcolare il valore basato non necessariamente su tutti i valori disponibili nell'intera finestra, ma solo sulle n celle spazialmente più vicine che hanno valori (in entrambe le immagini/Z- posizioni degli assi), dove n è un valore inferiore al numero totale di celle nella finestra 2D.

Al minuto, è molto più veloce calcolare per ogni cosa nella finestra, perché il mio mezzo di ordinamento per ottenere le celle n più vicine con i dati è la parte più lenta della funzione in quanto deve essere ripetuta ogni volta anche se il le distanze in termini di coordinate della finestra non cambiano. Non sono sicuro che sia necessario e sento che devo essere in grado di ottenere le distanze ordinate una volta e quindi mascherarle nel processo di selezione delle sole celle disponibili.

Ecco il mio codice per la selezione dei dati da utilizzare all'interno di una finestra della posizione della cella gap:

# radius will in reality be ~100 
radius = 2 
y,x = np.ogrid[-radius:radius+1, -radius:radius+1] 
dist = np.sqrt(x**2 + y**2) 
circle_template = dist > radius 

# this will in reality be a very large 3 dimensional array 
# representing daily images with some gaps, indicated by 0s 
dataStack = np.zeros((2,5,5)) 
dataStack[1] = (np.random.random(25) * 100).reshape(dist.shape) 
dataStack[0] = (np.random.random(25) * 100).reshape(dist.shape) 

testdata = dataStack[1] 
alternatedata = dataStack[0] 
random_gap_locations = (np.random.random(25) * 30).reshape(dist.shape) > testdata 
testdata[random_gap_locations] = 0 
testdata[radius, radius] = 0 

# in reality we will go through every gap (zero) location in the data 
# for each image and for each gap use slicing to get a window of 
# size (radius*2+1, radius*2+1) around it from each image, with the 
# gap being at the centre i.e. 
# testgaplocation = [radius, radius] 
# and the variables testdata, alternatedata below will refer to these 
# slices 

locations_to_exclude = np.logical_or(circle_template, np.logical_or 
            (testdata==0, alternatedata==0)) 
# the places that are inside the circular mask and where both images 
# have data 
locations_to_include = ~locations_to_exclude 
number_available = np.count_nonzero(locations_to_include) 

# we only want to do the interpolation calculations from the nearest n 
# locations that have data available, n will be ~100 in reality 
number_required = 3 

available_distances = dist[locations_to_include] 
available_data = testdata[locations_to_include] 
available_alternates = alternatedata[locations_to_include] 

if number_available > number_required: 
    # In this case we need to find the closest number_required of elements, based 
    # on distances recorded in dist, from available_data and available_alternates 
    # Having to repeat this argsort for each gap cell calculation is slow and feels 
    # like it should be avoidable 
    sortedDistanceIndices = available_distances.argsort(kind = 'mergesort',axis=None) 
    requiredIndices = sortedDistanceIndices[0:number_required] 
    selected_data = np.take(available_data, requiredIndices) 
    selected_alternates = np.take(available_alternates , requiredIndices) 
else: 
    # we just use available_data and available_alternates as they are... 

# now do stuff with the selected data to calculate a value for the gap cell 

Questo funziona, ma più della metà del tempo totale della funzione è preso nella argsort della mascherato dati sulla distanza spaziale. (~ 900 μS di 1,4 mS totali - e questa funzione verrà eseguita decine di miliardi di volte, quindi questa è una differenza importante!)

Sono sicuro che devo essere in grado di fare questo argsort una volta fuori dal funzione, quando la finestra della distanza spaziale è originariamente impostata e quindi includere gli indici di ordinamento nel mascheramento, per ottenere il primo comeManyToCalcolare gli indici senza dover ripetere l'ordinamento. La risposta potrebbe consistere nel mettere i vari bit da cui stiamo estraendo, in un array di record - ma non riesco a capire come, in caso affermativo. Qualcuno può vedere come posso rendere questa parte del processo più efficiente?

Answer 1

Così si vuole fare l'esterno di ordinamento del ciclo:

sorted_dist_idcs = dist.argsort(kind='mergesort', axis=None) 

Quindi, utilizzando alcune variabili dal codice originale, questo è quello che ho potuto venire con, se ci si sente ancora come un importante round- viaggio ..

loc_to_incl_sorted = locations_to_include.take(sorted_dist_idcs) 
sorted_dist_idcs_to_incl = sorted_dist_idcs[loc_to_incl_sorted] 

required_idcs = sorted_dist_idcs_to_incl[:number_required] 
selected_data = testdata.take(required_idcs) 
selected_alternates = alternatedata.take(required_idcs) 

Annotare il required_idcs si riferiscono a posizioni nel testdata e non available_data come nel codice originale. E questo frammento ho usato take allo scopo di indicizzare opportunamente l'array appiattito.

Answer 2

@moarningsun - grazie per il commento e la risposta. Questi mi hanno messo sulla strada giusta, ma non funzionano abbastanza per me quando il gap è il raggio < dal bordo dei dati: in questo caso uso una finestra attorno alla cella gap che viene "tagliata" ai limiti dei dati. In questa situazione gli indici riflettono la finestra "piena" e quindi non possono essere utilizzati per selezionare le celle dalla finestra delimitata.

Purtroppo ho modificato quella parte del mio codice quando ho chiarito la domanda originale ma si è rivelato rilevante.

Mi sono reso conto ora che se si utilizza argsort di nuovo sull'output di argsort si ottengono ranghi; vale a dire la posizione che ogni articolo avrebbe quando l'array complessivo è stato ordinato. Possiamo mascherarli in modo sicuro e quindi prendere il più piccolo (e farlo su un array strutturato per ottenere i dati corrispondenti allo stesso tempo).

Ciò implica un altro ordinamento all'interno del ciclo, ma in realtà è possibile utilizzare il partizionamento piuttosto che un ordinamento completo, poiché tutto ciò di cui abbiamo bisogno è il più piccolo num_required elementi. Se num_required è sostanzialmente inferiore al numero di elementi di dati, questo è molto più veloce rispetto a argsort.

Ad esempio con num_required = 80 e num_available = 15000 pieno argsort prende ~ 900μs mentre argpartition seguite da indice e fetta di ottenere il primo 80 prende ~ 110μs. Dobbiamo ancora fare il argsort per ottenere i ranghi all'inizio (piuttosto che il solo partizionamento basato sulla distanza) al fine di ottenere la stabilità del mergesort, e quindi ottenere il "giusto" quando la distanza non è unica.

Il mio codice come mostrato di seguito ora viene eseguito in ~ 610uS su dati reali, inclusi i calcoli effettivi che non sono mostrati qui. Ne sono felice ora, ma sembrano esserci molti altri fattori apparentemente minori che possono influenzare il tempo di esecuzione che è difficile da capire.

Ad esempio mettendo il circle_template nella matrice strutturata insieme dist, ranks, e un altro campo non illustrata, raddoppia il tempo di esecuzione della funzione globale (anche se non accede circle_template nel loop!). Ancora peggio, utilizzando np.partition nell'array strutturato con order=['ranks'] aumenta il runtime della funzione complessiva di quasi due ordini di grandezza rispetto all'utilizzo di np.argpartition come illustrato di seguito!

# radius will in reality be ~100 
radius = 2 
y,x = np.ogrid[-radius:radius+1, -radius:radius+1] 
dist = np.sqrt(x**2 + y**2) 
circle_template = dist > radius 
ranks = dist.argsort(axis=None,kind='mergesort').argsort().reshape(dist.shape) 
diam = radius * 2 + 1 

# putting circle_template in this array too doubles overall function runtime! 
fullWindowArray = np.zeros((diam,diam),dtype=[('ranks',ranks.dtype.str), 
             ('thisdata',dayDataStack.dtype.str), 
             ('alternatedata',dayDataStack.dtype.str), 
             ('dist',spatialDist.dtype.str)]) 
fullWindowArray['ranks'] = ranks 
fullWindowArray['dist'] = dist 

# this will in reality be a very large 3 dimensional array 
# representing daily images with some gaps, indicated by 0s 
dataStack = np.zeros((2,5,5)) 
dataStack[1] = (np.random.random(25) * 100).reshape(dist.shape) 
dataStack[0] = (np.random.random(25) * 100).reshape(dist.shape) 

testdata = dataStack[1] 
alternatedata = dataStack[0] 
random_gap_locations = (np.random.random(25) * 30).reshape(dist.shape) > testdata 
testdata[random_gap_locations] = 0 
testdata[radius, radius] = 0 

# in reality we will loop here to go through every gap (zero) location in the data 
# for each image 
gapz, gapy, gapx = 1, radius, radius 
desLeft, desRight = gapx - radius, gapx + radius+1 
desTop, desBottom = gapy - radius, gapy + radius+1 
extentB, extentR = dataStack.shape[1:] 

# handle the case where the gap is < search radius from the edge of 
# the data. If this is the case, we can't use the full 
# diam * diam window 
dataL = max(0, desLeft) 
maskL = 0 if desLeft >= 0 else abs(dataL - desLeft) 
dataT = max(0, desTop) 
maskT = 0 if desTop >= 0 else abs(dataT - desTop) 
dataR = min(desRight, extentR) 
maskR = diam if desRight <= extentR else diam - (desRight - extentR) 
dataB = min(desBottom,extentB) 
maskB = diam if desBottom <= extentB else diam - (desBottom - extentB) 

# get the slice that we will be working within 
# ranks, dist and circle are already populated 
boundedWindowArray = fullWindowArray[maskT:maskB,maskL:maskR] 
boundedWindowArray['alternatedata'] = alternatedata[dataT:dataB, dataL:dataR] 
boundedWindowArray['thisdata'] = testdata[dataT:dataB, dataL:dataR] 

locations_to_exclude = np.logical_or(boundedWindowArray['circle_template'], 
            np.logical_or 
            (boundedWindowArray['thisdata']==0, 
             boundedWindowArray['alternatedata']==0)) 
# the places that are inside the circular mask and where both images 
# have data 
locations_to_include = ~locations_to_exclude 
number_available = np.count_nonzero(locations_to_include) 

# we only want to do the interpolation calculations from the nearest n 
# locations that have data available, n will be ~100 in reality 
number_required = 3 

data_to_use = boundedWindowArray[locations_to_include] 

if number_available > number_required: 
    # argpartition seems to be v fast when number_required is 
    # substantially < data_to_use.size 
    # But partition on the structured array itself with order=['ranks'] 
    # is almost 2 orders of magnitude slower! 

    reqIndices = np.argpartition(data_to_use['ranks'],number_required)[:number_required] 
    data_to_use = np.take(data_to_use,reqIndices) 
else: 
    # we just use available_data and available_alternates as they are... 
    pass 
# now do stuff with the selected data to calculate a value for the gap cell