Échantillonner des points donnés de manière stochastique dans un espace 3D avec une distance minimale du plus proche voisin et une densité maximale

Jan 10 2021

J'ai des npoints dans un espace 3D. Je veux échantillonner de manière stochastique un sous-ensemble de points avec toutes les distances du plus proche voisin supérieures à r. La taille du sous m- ensemble est inconnue, mais je veux que les points échantillonnés soient aussi denses que possible.

Il y a des questions similaires, mais elles concernent toutes la génération de points, plutôt que l'échantillonnage à partir de points donnés.
Générez des points aléatoires dans l'espace 3D avec une distance minimale du plus proche voisin

Générer des points aléatoires en 3D avec une distance minimale entre chacun d'eux?

Disons que j'ai 300 points 3D aléatoires,

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

Quel est le moyen le plus rapide d'obtenir un sous-ensemble de mpoints avec une distance minimale du voisin le plus proche r = 1tout en maximisant m?

Réponses

2 DavidEisenstat Jan 14 2021 at 08:00

Il existe probablement un schéma d'approximation bicritère efficace, mais pourquoi s'embêter quand la programmation en nombres entiers est si rapide en moyenne?

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

Ce n'est pas la taille optimale d'un sous-ensemble, mais devrait être proche sans prendre beaucoup de temps, en utilisant KDTreepour optimiser les calculs de distance:

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

La suppression itérative du nœud le plus connecté n'est pas optimale (conserve environ 200 des 300 points), mais c'est probablement l'algorithme le plus rapide qui est proche de l'optimum (la seule chose plus rapide étant simplement de supprimer au hasard). Vous pourriez éventuellement @njit reduce_pairsrendre cette partie plus rapide (j'essaierai si j'ai le temps plus tard).

ShaunHan Jan 19 2021 at 07:07

Tester la réponse de @David Eisenstat avec 30 points donnés:

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

Lorsque la distance minimale attendue est de 1:

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

Vérifiez la distance minimale:

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

Changez la distance minimale attendue sur 2:

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

Vérifiez la distance minimale:

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