Próbkuj dane punkty stochastycznie w przestrzeni 3D z minimalną odległością do najbliższego sąsiada i maksymalną gęstością

Jan 10 2021

Mam npunkty w przestrzeni 3D. Chcę stochastycznie próbkować podzbiór punktów ze wszystkimi odległościami najbliższego sąsiada większymi niż r. Rozmiar podzbioru mjest nieznany, ale chcę, aby próbkowane punkty były jak najbardziej gęste.

Są podobne pytania, ale wszystkie dotyczą generowania punktów, a nie pobierania próbek z danych punktów.
Generuj losowe punkty w przestrzeni 3D z minimalną odległością od najbliższego sąsiada

Wygenerować 3-d losowe punkty z minimalną odległością między nimi?

Powiedzmy, że mam 300 losowych punktów 3D,

import numpy as np
n = 300
points = np.random.uniform(0, 10, size=(n, 3))

Jaki jest najszybszy sposób na uzyskanie podzbioru mpunktów z minimalną odległością do najbliższego sąsiada r = 1przy maksymalizacji m?

Odpowiedzi

2 DavidEisenstat Jan 14 2021 at 08:00

Prawdopodobnie istnieje skuteczny schemat przybliżenia dwukryteriów, ale po co zawracać sobie głowę, skoro programowanie liczb całkowitych jest średnio tak szybkie?

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

Nie jest to optymalnie duży podzbiór, ale powinien być blisko, ale nie zajmować zbyt wiele czasu, używając KDTreedo optymalizacji obliczeń odległości:

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()

Iteracyjne usuwanie najbardziej połączonego węzła nie jest optymalne (zachowuje około 200 z 300 punktów), ale jest prawdopodobnie najszybszym algorytmem, który jest bliski optymalnego (jedyną szybszą rzeczą jest usuwanie losowe). Możesz ewentualnie @njit reduce_pairsprzyspieszyć tę część (spróbuję, jeśli będę miał czas później).

ShaunHan Jan 19 2021 at 07:07

Testowanie odpowiedzi @David Eisenstat z 30 punktami:

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 ],])

Gdy oczekiwana minimalna odległość wynosi 1:

subset1 = subset_David_Eisenstat(points, r=1.)
print(len(subset1))
# Output: 18

Sprawdź minimalną odległość:

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

Zmień oczekiwaną minimalną odległość na 2:

subset2 = subset_David_Eisenstat(points, r=2.)
print(len(subset2))
# Output: 10

Sprawdź minimalną odległość:

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