Amostra pontos dados estocasticamente em um espaço 3D com distância mínima do vizinho mais próximo e densidade máxima
Tenho n
pontos em um espaço 3D. Quero amostrar estocasticamente um subconjunto de pontos com todas as distâncias dos vizinhos mais próximos maiores que r
. O tamanho do subconjunto m
é desconhecido, mas quero que os pontos amostrados sejam os mais densos possível.
Existem questões semelhantes, mas todas tratam de gerar pontos, em vez de amostrar a partir de determinados pontos.
Gere pontos aleatórios no espaço 3D com distância mínima do vizinho mais próximo
Gerar pontos aleatórios 3-d com distância mínima entre cada um deles?
Digamos que eu tenha 300 pontos 3D aleatórios,
import numpy as np
n = 300
points = np.random.uniform(0, 10, size=(n, 3))
Qual é a maneira mais rápida de obter um subconjunto de m
pontos com a distância mínima do vizinho mais próximo r = 1
ao maximizar m
?
Respostas
Provavelmente existe um esquema de aproximação bicritério eficiente, mas por que se preocupar quando a programação de inteiros é tão rápida em média?
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()]))
Isso não é idealmente grande em um subconjunto, mas deve estar próximo, embora não demore muito, usando KDTree
para otimizar os cálculos de distância:
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()
Remover iterativamente o nó mais conectado não é ideal (mantém cerca de 200 dos 300 pontos), mas é provavelmente o algoritmo mais rápido que está próximo do ideal (a única coisa mais rápida é apenas remover aleatoriamente). Você poderia @njit
reduce_pairs
tornar essa parte mais rápida (vou tentar se eu tiver tempo mais tarde).
Testando a resposta de @David Eisenstat com 30 pontos dados:
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 ],])
Quando a distância mínima esperada é 1:
subset1 = subset_David_Eisenstat(points, r=1.)
print(len(subset1))
# Output: 18
Verifique a distância mínima:
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
Altere a distância mínima esperada para 2:
subset2 = subset_David_Eisenstat(points, r=2.)
print(len(subset2))
# Output: 10
Verifique a distância mínima:
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