Probieren Sie gegebene Punkte stochastisch in einem 3D-Raum mit minimalem Abstand zum nächsten Nachbarn und maximaler Dichte aus
Ich habe n
Punkte in einem 3D-Raum. Ich möchte eine Teilmenge von Punkten stochastisch abtasten, wobei alle Abstände zum nächsten Nachbarn größer als sind r
. Die Größe der Teilmenge m
ist unbekannt, aber ich möchte, dass die abgetasteten Punkte so dicht wie möglich sind.
Es gibt ähnliche Fragen, aber es geht nur darum, Punkte zu generieren, anstatt aus bestimmten Punkten zu stichproben.
Generieren Sie zufällige Punkte im 3D-Raum mit minimalem Abstand zum nächsten Nachbarn
3D-Zufallspunkte mit minimalem Abstand zwischen ihnen generieren?
Angenommen, ich habe 300 zufällige 3D-Punkte.
import numpy as np
n = 300
points = np.random.uniform(0, 10, size=(n, 3))
Was ist der schnellste Weg, um eine Teilmenge von m
Punkten mit minimalem Abstand zum nächsten Nachbarn zu erhalten und r = 1
gleichzeitig zu maximieren m
?
Antworten
Es gibt wahrscheinlich ein effizientes Bicriteria-Approximationsschema, aber warum sollte man sich die Mühe machen, wenn die Ganzzahlprogrammierung im Durchschnitt so schnell ist?
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()]))
Dies ist nicht optimal groß für eine Teilmenge, sollte aber nahe sein, ohne sehr lange zu dauern KDTree
, um die Entfernungsberechnungen zu optimieren:
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()
Das iterative Entfernen des am meisten verbundenen Knotens ist nicht optimal (behält ungefähr 200 der 300 Punkte bei), ist jedoch wahrscheinlich der schnellste Algorithmus, der nahezu optimal ist (das einzige, was schneller ist, ist das zufällige Entfernen). Sie könnten @njit
reduce_pairs
diesen Teil möglicherweise schneller machen (werde es versuchen, wenn ich später Zeit habe).
Testen der Antwort von @David Eisenstat mit 30 gegebenen Punkten:
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 ],])
Wenn der erwartete Mindestabstand 1 beträgt:
subset1 = subset_David_Eisenstat(points, r=1.)
print(len(subset1))
# Output: 18
Überprüfen Sie den Mindestabstand:
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
Ändern Sie den erwarteten Mindestabstand auf 2:
subset2 = subset_David_Eisenstat(points, r=2.)
print(len(subset2))
# Output: 10
Überprüfen Sie den Mindestabstand:
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