Contoh titik yang diberikan secara stokastik dalam ruang 3D dengan jarak tetangga terdekat minimum dan kepadatan maksimum

Jan 10 2021

Saya memiliki npoin dalam ruang 3D. Saya ingin mengambil sampel secara stokastik subset titik dengan semua jarak tetangga terdekat lebih besar dari r. Ukuran subset mtidak diketahui, tetapi saya ingin titik sampelnya sepadat mungkin.

Ada pertanyaan serupa, tetapi semuanya tentang menghasilkan poin, bukan pengambilan sampel dari poin yang diberikan.
Hasilkan titik acak dalam ruang 3D dengan jarak tetangga terdekat minimum

Hasilkan titik acak 3-d dengan jarak minimum di antara masing-masing titik?

Katakanlah saya memiliki 300 titik 3D acak,

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

Apa cara tercepat untuk mendapatkan subset mpoin dengan jarak tetangga terdekat minimum r = 1sambil memaksimalkan m?

Jawaban

2 DavidEisenstat Jan 14 2021 at 08:00

Mungkin ada skema perkiraan bikriteria yang efisien, tetapi mengapa repot-repot ketika pemrograman bilangan bulat rata-rata begitu cepat?

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

Ini bukan bagian yang besar secara optimal, tetapi harus dekat sementara tidak membutuhkan waktu yang lama, gunakan KDTreeuntuk mengoptimalkan penghitungan jarak:

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

Menghapus node yang paling terhubung secara berulang tidak optimal (menyimpan sekitar 200 dari 300 poin), tetapi kemungkinan merupakan algoritme tercepat yang mendekati optimal (satu-satunya hal yang lebih cepat adalah menghapus secara acak). Anda mungkin bisa @njit reduce_pairsmembuat bagian itu lebih cepat (akan coba jika saya punya waktu nanti).

ShaunHan Jan 19 2021 at 07:07

Menguji jawaban @David Eisenstat dengan 30 poin yang diberikan:

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

Jika jarak minimum yang diharapkan adalah 1:

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

Periksa jarak minimum:

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

Ubah jarak minimum yang diharapkan menjadi 2:

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

Periksa jarak minimum:

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