Minimum en yakın komşu mesafesi ve maksimum yoğunluk ile 3 boyutlu uzayda stokastik olarak örnek verilen noktalar

Jan 10 2021

n3B 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 mbilinmiyor, 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 mminimum en yakın komşu mesafesi ile bir nokta alt kümesi elde etmenin en hızlı yolu nedir ?r = 1m

Yanıtlar

2 DavidEisenstat Jan 14 2021 at 08:00

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()]))
1 DanielF Jan 18 2021 at 16:16

Bu, bir alt küme için optimal olarak büyük değildir, ancak KDTreemesafe 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_pairsbölümü daha hızlı yapabilirsiniz (daha sonra zamanım olursa denerim).

ShaunHan Jan 19 2021 at 07:07

@ 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