Minimum en yakın komşu mesafesi ve maksimum yoğunluk ile 3 boyutlu uzayda stokastik olarak örnek verilen noktalar
n
3B alanda puanlarım var . Tüm en yakın komşu uzaklıkları şundan daha büyük olan bir nokta alt kümesini stokastik olarak örneklemek istiyorum r
. Alt kümenin boyutu m
bilinmiyor, ancak örneklenen noktaların olabildiğince yoğun olmasını istiyorum.
Benzer sorular var, ancak hepsi belirli noktalardan örnekleme yapmaktan ziyade noktalar oluşturmakla ilgili.
Minimum en yakın komşu mesafesi ile 3B alanda rastgele noktalar oluşturun
Her biri arasında minimum mesafe olacak şekilde 3 boyutlu rastgele noktalar üretilsin mi?
300 rastgele 3B noktam olduğunu varsayalım,
import numpy as np
n = 300
points = np.random.uniform(0, 10, size=(n, 3))
En üst düzeye çıkarırken m
minimum en yakın komşu mesafesi ile bir nokta alt kümesi elde etmenin en hızlı yolu nedir ?r = 1
m
Yanıtlar
Muhtemelen verimli bir iki kriterli yaklaşım şeması vardır, ancak tamsayı programlama ortalamada bu kadar hızlıyken neden uğraşasınız ki?
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()]))
Bu, bir alt küme için optimal olarak büyük değildir, ancak KDTree
mesafe hesaplamalarını optimize etmek için kullanarak çok uzun sürmemekle birlikte yakın olmalıdır :
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()
Yinelemeli olarak en çok bağlanan düğümü kaldırmak optimal değildir (300 noktanın yaklaşık 200'ünü tutar), ancak muhtemelen optimuma yakın olan en hızlı algoritmadır (daha hızlı olan tek şey rastgele kaldırmaktır). Bu @njit
reduce_pairs
bölümü daha hızlı yapabilirsiniz (daha sonra zamanım olursa denerim).
@ David Eisenstat'ın cevabını verilen 30 puanla test etmek:
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 ],])
Beklenen minimum mesafe 1 olduğunda:
subset1 = subset_David_Eisenstat(points, r=1.)
print(len(subset1))
# Output: 18
Minimum mesafeyi kontrol edin:
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
Beklenen minimum mesafeyi 2 olarak değiştirin:
subset2 = subset_David_Eisenstat(points, r=2.)
print(len(subset2))
# Output: 10
Minimum mesafeyi kontrol edin:
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