python - kết hợp argsort với mặt nạ để lấy các giá trị gần nhất trong cửa sổ di chuyển

Question

Tôi có một số mã để tính giá trị bị thiếu trong một hình ảnh, dựa trên các giá trị lân cận trong cửa sổ tròn 2D. Nó cũng sử dụng các giá trị từ một hoặc nhiều hình ảnh liền kề tạm thời tại cùng một vị trí (tức là cùng một cửa sổ 2D được dịch chuyển trong thứ nguyên thứ 3).

Đối với mỗi vị trí bị thiếu, tôi cần tính giá trị dựa trên không nhất thiết trên tất cả các giá trị có sẵn trong toàn bộ cửa sổ, nhưng chỉ trên các ô không gian gần nhất có giá trị (trong cả hai hình ảnh/Z- vị trí trục), trong đó n là một số giá trị nhỏ hơn tổng số ô trong cửa sổ 2D.

Vào phút, việc tính toán mọi thứ trong cửa sổ nhanh hơn, vì phương tiện sắp xếp để lấy các ô gần nhất với dữ liệu là phần chậm nhất của hàm vì nó phải được lặp lại mỗi lần mặc dù khoảng cách về tọa độ cửa sổ không thay đổi. Tôi không chắc chắn điều này là cần thiết và cảm thấy tôi phải có khả năng để có được khoảng cách được sắp xếp một lần, và sau đó mặt nạ những người trong quá trình chỉ chọn các tế bào có sẵn.

Dưới đây là mã của tôi để lựa chọn dữ liệu để sử dụng trong một cửa sổ của địa điểm di động khoảng cách:

# 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 

này hoạt động, nhưng hơn một nửa tổng thời gian của hàm được thực hiện trong argsort của đeo mặt nạ dữ liệu khoảng cách không gian. (~ 900uS trong tổng số 1,4mS - và chức năng này sẽ chạy hàng chục tỷ lần, vì vậy đây là một sự khác biệt quan trọng!)

Tôi chắc chắn rằng tôi phải có khả năng thực hiện điều này một lần bên ngoài chức năng, khi cửa sổ khoảng cách không gian ban đầu được thiết lập và sau đó bao gồm các chỉ mục sắp xếp đó trong mặt nạ, để có chỉ mục howManyToCalculate đầu tiên mà không phải thực hiện lại việc sắp xếp. Câu trả lời có thể liên quan đến việc đưa các bit khác nhau mà chúng ta đang trích xuất từ, vào một mảng bản ghi - nhưng tôi không thể tìm ra cách, nếu có. Bất cứ ai có thể xem làm thế nào tôi có thể làm cho một phần của quá trình này hiệu quả hơn?

Answer 1

Vì vậy, bạn muốn làm bên ngoài sắp xếp vòng lặp:

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

Sau đó, sử dụng một số biến từ mã gốc, đây là những gì tôi có thể đưa ra, mặc dù nó vẫn cảm thấy giống như một round- lớn chuyến đi ..

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) 

Lưu ý required_idcs tham khảo đến các địa điểm trong testdata và không available_data như trong mã gốc. Và đoạn mã này tôi đã sử dụng take với mục đích thuận tiện lập chỉ mục mảng được làm phẳng.

Answer 2

@moarningsun - cảm ơn nhận xét và câu trả lời. Điều này khiến tôi đi đúng hướng, nhưng không hoàn toàn phù hợp với tôi khi khoảng cách là bán kính < từ cạnh của dữ liệu: trong trường hợp này tôi sử dụng một cửa sổ xung quanh ô khoảng cách được "cắt" thành các giới hạn dữ liệu. Trong tình huống này các chỉ số phản ánh cửa sổ "đầy đủ" và do đó không thể được sử dụng để chọn các ô từ cửa sổ bị chặn.

Thật không may là tôi đã chỉnh sửa một phần mã của mình khi tôi làm rõ câu hỏi gốc nhưng hóa ra lại có liên quan.

Tôi đã nhận ra ngay bây giờ nếu bạn sử dụng lại argsort trên đầu ra của argsort thì bạn sẽ nhận được thứ hạng; tức là vị trí mà mỗi mục sẽ có khi mảng tổng thể được sắp xếp. Chúng ta có thể an toàn che giấu chúng và sau đó lấy number_required của chúng nhỏ nhất (và thực hiện điều này trên một mảng có cấu trúc để nhận dữ liệu tương ứng cùng một lúc).

Điều này ngụ ý một loại khác trong vòng lặp, nhưng trên thực tế chúng tôi có thể sử dụng phân vùng thay vì sắp xếp đầy đủ, vì tất cả những gì chúng tôi cần là các mục num_required nhỏ nhất. Nếu num_required là ít hơn đáng kể so với số lượng các mục dữ liệu thì điều này nhanh hơn nhiều so với thực hiện argsort.

Ví dụ với num_required = 80 và num_available = 15000 đầy đủ argsort mất ~ 900μs trong khi argpartition tiếp theo chỉ số và lát để có được 80 đầu tiên diễn ~ 110μs. Chúng ta vẫn cần phải làm argsort để có được thứ hạng ngay từ đầu (thay vì chỉ phân vùng dựa trên khoảng cách) để có được sự ổn định của mergesort, và do đó có được "đúng" khi khoảng cách không phải là duy nhất.

Mã của tôi như được hiển thị bên dưới hiện chạy trong ~ 610uS trên dữ liệu thực, bao gồm các tính toán thực tế không được hiển thị ở đây. Tôi hạnh phúc với điều đó bây giờ, nhưng dường như có một số yếu tố khác dường như nhỏ có thể có ảnh hưởng đến thời gian chạy khó hiểu.

Ví dụ đặt circle_template trong mảng cấu trúc cùng với dist, ranks, và một lĩnh vực không được hiển thị ở đây, tăng gấp đôi thời gian chạy của hàm tổng thể (thậm chí nếu chúng ta không truy cập circle_template trong vòng lặp!). Thậm chí tệ hơn, sử dụng np.partition trên mảng có cấu trúc với order=['ranks'] tăng thời gian chạy chức năng tổng thể lên gần như hai đơn đặt hàng có độ lớn và sử dụng np.argpartition như được hiển thị bên dưới!

# 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