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
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 m
là 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 = 1
trong khi tối đa hóa là m
gì?
Trả lời
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()]))
Đâ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_pairs
thực hiện phần đó nhanh hơn (sẽ thử nếu tôi có thời gian sau).
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