python - łączenie argsort z maskowaniem w celu uzyskania najbliższych wartości w ruchu okna

Question

Mam kod do obliczania brakujących wartości w obrazie, w oparciu o sąsiadujące wartości w okólniku 2D. Wykorzystuje również wartości z jednego lub więcej sąsiadujących ze sobą obrazów w tych samych lokalizacjach (to jest to samo okno 2D przesunięte w trzecim wymiarze).

Dla każdej pozycji, której brakuje, muszę obliczyć wartość na podstawie niekoniecznie na wszystkich wartościach dostępnych w całym oknie, ale tylko na najbliżej położonych n komórkach, które mają wartości (na obu obrazach/Z- pozycje osi), gdzie n jest wartością mniejszą niż całkowita liczba komórek w oknie 2D.

Z minuty na minutę, wszystko jest znacznie szybsze do obliczenia w oknie, ponieważ mój sposób sortowania w celu uzyskania najbliższych n komórek z danymi jest najwolniejszą częścią funkcji, ponieważ musi być powtarzany za każdym razem, mimo że odległości w zakresie współrzędnych okna nie zmieniają się. Nie jestem pewien, czy jest to konieczne i czuję, że muszę być w stanie uzyskać posortowane odległości jeden raz, a następnie zamaskować je w procesie tylko wybierania dostępnych komórek.

Oto mój kod do wybierania danych do wykorzystania w oknie lokalizacji komórek 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 

To działa, ale ponad połowę całkowitego czasu funkcji jest zrobione w argsort z zamaskowanym dane przestrzennej odległości. (~ 900uS o łącznej wartości 1,4 ms - i ta funkcja będzie działała dziesiątki miliardów razy, więc to jest ważna różnica!)

Jestem pewna, że ​​muszę być w stanie zrobić to po argsortie na zewnątrz funkcja, gdy okno odległości przestrzennej jest pierwotnie skonfigurowane, a następnie włącza te wskaźniki sortowania w maskowaniu, aby uzyskać pierwsze informacje o tym, jak obliczaćManyToCalculate bez konieczności ponownego sortowania. Odpowiedź może polegać na umieszczeniu różnych bitów, z których wydobywamy, na tablicę rekordów - ale nie wiem, jak to zrobić. Czy ktokolwiek może zobaczyć, jak sprawić, by ta część procesu była bardziej wydajna?

Answer 1

więc chcesz zrobić sortowania zewnątrz pętli:

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

Następnie za pomocą kilku zmiennych z oryginalnego kodu, to co mogłam wymyślić, choć nadal czuje się jak główny okrągłodennej podróż ..

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) 

Zanotuj required_idcs odnoszą się do miejsc, w testdata i nie available_data jak w oryginalnym kodzie. I ten fragment użyłem take w celu wygodnego indeksowania spłaszczonej macierzy.

Answer 2

@moarningsun - dzięki za komentarz i odpowiedź. To mnie poprawiło, ale nie działa dla mnie, gdy różnica wynosi < promień od krawędzi danych: w tym przypadku używam okna wokół komórki szczeliny, która jest "przycięta" do granic danych. W tej sytuacji indeksy odzwierciedlają "pełne" okno i dlatego nie można ich użyć do wybrania komórek z ograniczonego okna.

Niestety, zredagowałem tę część mojego kodu, kiedy wyjaśniłem pierwotne pytanie, ale okazało się to istotne.

Zdałem sobie sprawę, że jeśli ponownie użyjesz argsort na wyjściu z argsort, otrzymasz szereg; tj. pozycja, którą każdy element miałby po posortowaniu całej tablicy. Możemy je bezpiecznie zamaskować, a następnie pobrać najmniejszy z nich (i zrobić to na tablicy strukturalnej, aby uzyskać odpowiednie dane w tym samym czasie).

Oznacza to inny rodzaj w pętli, ale w rzeczywistości możemy użyć partycjonowanie zamiast pełnej rodzaju, ponieważ wszystko, czego potrzebujemy jest najmniejszym num_required przedmiotów. Jeśli num_required jest znacznie mniejsza niż liczba elementów danych, jest to znacznie szybsze niż wykonanie argsort.

Na przykład num_required = 80, a = 15000 num_available pełne argsort odbywa ~ 900μs, natomiast argpartition następnie przez indeks i plaster aby pierwsza trwa 80 ~ 110μs. Nadal musimy wykonać argsort, aby uzyskać stopnie od samego początku (a nie tylko partycjonowanie na podstawie odległości), aby uzyskać stabilność połączenia, a tym samym uzyskać "właściwą", gdy odległość nie jest unikalna.

Mój kod, jak pokazano poniżej, teraz działa w ~ 610uS na rzeczywistych danych, w tym rzeczywistych obliczeń, które nie są tutaj pokazane. Cieszę się z tego teraz, ale wydaje się, że istnieje kilka innych pozornie drobnych czynników, które mogą mieć wpływ na środowisko uruchomieniowe, co jest trudne do zrozumienia.

Na przykład oddanie circle_template w uporządkowany tablicy wraz z dist, ranks, a inna dziedzina nie pokazano tutaj, podwaja czas pracy ogólnej funkcji (nawet jeśli nie mamy dostępu circle_template w pętli!). Co gorsza, użycie struktury np.partition na tablicy strukturalnej z order=['ranks'] zwiększa ogólny czas działania funkcji o prawie o dwa rzędy wielkości w porównaniu z użyciem np.argpartition, jak pokazano poniżej!

# 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