Lấy mẫu các điểm đã cho một cách ngẫu nhiên trong không gian 3D với khoảng cách hàng xóm gần nhất tối thiểu và mật độ tối đa

Jan 10 2021

Tôi có nđiểm trong không gian 3D. Tôi muốn lấy mẫu ngẫu nhiên một tập hợp con các điểm với tất cả các khoảng cách lân cận gần nhất đều lớn hơn r. Kích thước của tập hợp con mlà không xác định, nhưng tôi muốn các điểm được lấy mẫu càng dày đặc càng tốt.

Có những câu hỏi tương tự, nhưng tất cả đều là về việc tạo ra các điểm, thay vì lấy mẫu từ các điểm đã cho.
Tạo các điểm ngẫu nhiên trong không gian 3D với khoảng cách hàng xóm gần nhất tối thiểu

Tạo các điểm ngẫu nhiên 3-d với khoảng cách giữa mỗi điểm là nhỏ nhất?

Giả sử tôi có 300 điểm 3D ngẫu nhiên,

import numpy as np
n = 300
points = np.random.uniform(0, 10, size=(n, 3))

Cách nhanh nhất để có được một tập hợp con các mđiểm với khoảng cách nhỏ nhất-láng giềng gần nhất r = 1trong khi tối đa hóa là mgì?

Trả lời

2 DavidEisenstat Jan 14 2021 at 08:00

Có thể có một lược đồ xấp xỉ nhị phân hiệu quả, nhưng tại sao phải bận tâm khi lập trình số nguyên trung bình quá nhanh?

import numpy as np

n = 300
points = np.random.uniform(0, 10, size=(n, 3))

from ortools.linear_solver import pywraplp

solver = pywraplp.Solver.CreateSolver("SCIP")
variables = [solver.BoolVar("x[{}]".format(i)) for i in range(n)]
solver.Maximize(sum(variables))
for j, q in enumerate(points):
    for i, p in enumerate(points[:j]):
        if np.linalg.norm(p - q) <= 1:
            solver.Add(variables[i] + variables[j] <= 1)
solver.EnableOutput()
solver.Solve()
print(len([i for (i, variable) in enumerate(variables) if variable.SolutionValue()]))
1 DanielF Jan 18 2021 at 16:16

Đây không phải là lớn tối ưu của một tập hợp con, nhưng phải gần trong khi không mất nhiều thời gian, sử dụng KDTreeđể tối ưu hóa các tính toán khoảng cách:

from scipy.spatial import KDTree
import numpy as np

def space_sample(n = 300, low = 0, high = 10, dist = 1):
    points = np.random.uniform(low, high, size=(n, 3))
    k = KDTree(points)
    pairs = np.array(list(k.query_pairs(dist)))
    
    def reduce_pairs(pairs, remove = []):  #iteratively remove the most connected node
        p = pairs[~np.isin(pairs, remove).any(1)]
        if p.size >0:
            count = np.bincount(p.flatten(), minlength = n)
            r = remove + [count.argmax()]
            return reduce_pairs(p, r)
        else:
            return remove
    
    return np.array([p for i, p in enumerate(points) if not(i in reduce_pairs(pairs))])

subset = space_sample()

Loại bỏ lặp đi lặp lại nút được kết nối nhiều nhất không phải là tối ưu (giữ khoảng 200 trong số 300 điểm), nhưng có thể là thuật toán nhanh nhất gần với mức tối ưu (điều duy nhất nhanh hơn là chỉ xóa ngẫu nhiên). Bạn có thể @njit reduce_pairsthực hiện phần đó nhanh hơn (sẽ thử nếu tôi có thời gian sau).

ShaunHan Jan 19 2021 at 07:07

Kiểm tra câu trả lời của @David Eisenstat với 30 điểm cho trước:

from ortools.linear_solver import pywraplp
import numpy as np

def subset_David_Eisenstat(points, r):
    solver = pywraplp.Solver.CreateSolver("SCIP")
    variables = [solver.BoolVar("x[{}]".format(i)) for i in range(len(points))]
    solver.Maximize(sum(variables))
    for j, q in enumerate(points):
        for i, p in enumerate(points[:j]):
            if np.linalg.norm(p - q) <= r:
                solver.Add(variables[i] + variables[j] <= 1)
    solver.EnableOutput()
    solver.Solve()
    indices = [i for (i, variable) in enumerate(variables) if variable.SolutionValue()]
    return points[indices]

points = np.array(
[[ 7.32837882, 12.12765786, 15.01412241],
 [ 8.25164031, 11.14830379, 15.01412241],
 [ 8.21790113, 13.05647987, 13.05647987],
 [ 7.30031002, 13.08276009, 14.05452502],
 [ 9.18056467, 12.0800735 , 13.05183844],
 [ 9.17929647, 11.11270337, 14.03027534],
 [ 7.64737905, 11.48906945, 15.34274827],
 [ 7.01315886, 12.77870699, 14.70301668],
 [ 8.88132414, 10.81243313, 14.68685022],
 [ 7.60617372, 13.39792166, 13.39792166],
 [ 8.85967682, 12.72946394, 12.72946394],
 [ 9.50060705, 11.43361294, 13.37866092],
 [ 8.21790113, 12.08471494, 14.02824481],
 [ 7.32837882, 12.12765786, 16.98587759],
 [ 8.25164031, 11.14830379, 16.98587759],
 [ 7.30031002, 13.08276009, 17.94547498],
 [ 8.21790113, 13.05647987, 18.94352013],
 [ 9.17929647, 11.11270337, 17.96972466],
 [ 9.18056467, 12.0800735 , 18.94816156],
 [ 7.64737905, 11.48906945, 16.65725173],
 [ 7.01315886, 12.77870699, 17.29698332],
 [ 8.88132414, 10.81243313, 17.31314978],
 [ 7.60617372, 13.39792166, 18.60207834],
 [ 8.85967682, 12.72946394, 19.27053606],
 [ 9.50060705, 11.43361294, 18.62133908],
 [ 8.21790113, 12.08471494, 17.97175519],
 [ 7.32837882, 15.01412241, 12.12765786],
 [ 8.25164031, 15.01412241, 11.14830379],
 [ 7.30031002, 14.05452502, 13.08276009],
 [ 9.18056467, 13.05183844, 12.0800735 ],])

Khi khoảng cách tối thiểu dự kiến ​​là 1:

subset1 = subset_David_Eisenstat(points, r=1.)
print(len(subset1))
# Output: 18

Kiểm tra khoảng cách tối thiểu:

from scipy.spatial.distance import cdist
dist = cdist(subset1, subset1, metric='euclidean')
# Delete diagonal
res = dist[~np.eye(dist.shape[0],dtype=bool)].reshape(dist.shape[0],-1)
print(np.min(res))
# Output: 1.3285513450926985

Thay đổi khoảng cách tối thiểu dự kiến ​​thành 2:

subset2 = subset_David_Eisenstat(points, r=2.)
print(len(subset2))
# Output: 10

Kiểm tra khoảng cách tối thiểu:

from scipy.spatial.distance import cdist
dist = cdist(subset2, subset2, metric='euclidean')
# Delete diagonal
res = dist[~np.eye(dist.shape[0],dtype=bool)].reshape(dist.shape[0],-1)
print(np.min(res))
# Output: 2.0612041004376223