Muestra los puntos dados estocásticamente en un espacio 3D con una distancia mínima al vecino más cercano y una densidad máxima

Question

Tengo npuntos en un espacio 3D. Quiero muestrear estocásticamente un subconjunto de puntos con todas las distancias de los vecinos más cercanos mayores que r. Se mdesconoce el tamaño del subconjunto , pero quiero que los puntos muestreados sean lo más densos posible.

Hay preguntas similares, pero todas tienen que ver con generar puntos, en lugar de tomar muestras de puntos dados.
Genere puntos aleatorios en el espacio 3D con una distancia mínima al vecino más cercano

¿Generar puntos aleatorios tridimensionales con una distancia mínima entre cada uno de ellos?

Digamos que tengo 300 puntos 3D aleatorios,

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

¿Cuál es la forma más rápida de obtener un subconjunto de mpuntos con la distancia mínima al vecino más cercano r = 1mientras se maximiza m?

Answer 1

Probablemente exista un esquema de aproximación de bicriterio eficiente, pero ¿por qué molestarse cuando la programación de enteros es tan rápida en promedio?

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

Answer 2

Esto no es un subconjunto demasiado grande de manera óptima, pero debe estar cerca sin tomar mucho tiempo, y se usa KDTreepara optimizar los cálculos de distancia:

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

Eliminar iterativamente el nodo más conectado no es óptimo (mantiene alrededor de 200 de los 300 puntos), pero es probable que sea el algoritmo más rápido cercano al óptimo (lo único más rápido es simplemente eliminar al azar). Posiblemente podría @njit reduce_pairshacer esa parte más rápido (lo intentaré si tengo tiempo más tarde).

Answer 3

Probando la respuesta de @David Eisenstat con 30 puntos dados:

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

Cuando la distancia mínima esperada es 1:

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

Compruebe la distancia mínima:

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

Cambie la distancia mínima esperada a 2:

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

Compruebe la distancia mínima:

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