Вершины выпуклой оболочки n-мерного точечного множества
У меня есть заданный набор точек в размерности n.из них я хочу найти те, которые являются вершинами (углами) выпуклой оболочки. Я хочу решить эту проблему с помощью Python (но могу вызвать и другие программы).
Правка: все координаты являются натуральными числами. В качестве вывода я ищу индексы вершин.
Гугл обычно давал задачу в 2D, или просил перечислить грани, что вычислительно намного сложнее.
Мои собственные попытки так далеко
-
scipy.spatial.ConvexHull()
: выдает ошибку для моего текущего примера. И я думаю, я читал, что это не работает для измерения выше 10. Мой начальник тоже не советовал этого делать. -
Normaliz
(как часть polymake): работает, но слишком медленно. Но, может быть, я сделал что-то не так.import PyNormaliz def find_column_indices(A,B): return [i for i in range(A.shape[1]) if list(A[:,i]) in B.T.tolist()] def convex_hull(A): poly = PyNormaliz.Cone(polytope = A.T.tolist()) hull_cone = poly.IntegerHull() hull_vertices = np.array([entry[:-1] for entry in hull_cone.VerticesOfPolyhedron()]).T hull_indices = find_column_indices(A, hull_vertices) return hull_indices
-
Решение с линейными программами: работает, но полностью не оптимизировано, поэтому я думаю, что должно быть лучшее решение.
import subprocess from multiprocessing import Pool, cpu_count from scipy.optimize import linprog import numpy as np def is_in_convex_hull(arg): A,v = arg m = A.shape[1] A_ub = -np.eye(m,dtype = np.int) b_ub = np.zeros(m) res = linprog([1 for _ in range(m)],A_ub,b_ub,A,v) return res['success'] def convex_hull2(A): pool = Pool(processes = cpu_count()) res = pool.map(is_in_convex_hull,[(np.delete(A,i,axis=1),A[:,i]) for i in range(A.shape[1])]) return [i for i in range(A.shape[1]) if not res[i]]
Пример:
A = array([[ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
...: [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 2, 2, 4, 6, 6, 6, 8, 1, 2, 2, 2, 2, 3, 2, 1, 2, 1, 2, 1, 1, 1, 1, 2, 2, 2, 3, 1, 2, 2, 1, 2, 1, 1, 1, 2, 1, 2],
...: [ 0, 0, 0, 0, 0, 2, 2, 4, 6, 0, 0, 2, 4, 0, 0, 2, 2, 2, 1, 2, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 2, 0, 1, 0, 1],
...: [ 0, 0, 2, 4, 4, 2, 2, 0, 0, 0, 6, 2, 0, 4, 0, 2, 4, 0, 1, 1, 2, 2, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 2, 1, 2, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1],
...: [ 0, 6, 0, 2, 4, 0, 6, 4, 2, 2, 0, 0, 8, 4, 8, 4, 0, 2, 4, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 2, 2, 2, 2, 2, 4, 2, 2, 1, 1, 1, 2, 3, 2, 2, 2, 2, 1, 2],
...: [ 0, 2, 14, 0, 4, 6, 0, 0, 4, 0, 2, 0, 4, 4, 4, 0, 0, 2, 2, 2, 2, 2, 2, 1, 2, 4, 1, 3, 2, 1, 1, 1, 1, 2, 1, 4, 1, 1, 0, 1, 1, 1, 2, 3, 1, 1, 1, 1],
...: [ 0, 0, 0, 2, 6, 0, 4, 6, 0, 0, 6, 2, 2, 0, 0, 2, 2, 0, 1, 1, 2, 1, 2, 1, 1, 1, 1, 1, 1, 1, 2, 2, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 0, 1],
...: [ 0, 2, 2, 12, 2, 0, 0, 2, 0, 8, 2, 4, 0, 4, 0, 4, 0, 0, 2, 1, 2, 1, 1, 2, 1, 1, 4, 2, 1, 2, 3, 1, 3, 2, 2, 2, 1, 1, 2, 1, 1, 1, 1, 1, 3, 1, 1, 2],
...: [ 0, 8, 2, 0, 0, 2, 2, 4, 4, 4, 2, 4, 0, 0, 2, 0, 0, 6, 2, 2, 1, 1, 1, 3, 2, 2, 1, 2, 2, 1, 2, 1, 2, 1, 3, 1, 2, 1, 1, 1, 1, 3, 1, 3, 2, 2, 0, 1]])
Доходы время
In [44]: %timeit convex_hull(A)
1.79 s ± 16.8 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
In [45]: %timeit convex_hull2(A)
337 ms ± 3.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
Для чуть большего примера соотношение было хуже, так что оно также не может быть объяснено параллелизацией.
Любая помощь или улучшение приветствуются.
1 ответ:
Вы можете изменить свой метод
is_in_convex_hull
следующим образом:def convex_hull(A): vertices = A.T.tolist() vertices = [ i + [ 1 ] for i in vertices ] poly = PyNormaliz.Cone(vertices = vertices) hull_vertices = np.array([entry[:-1] for entry in poly.VerticesOfPolyhedron()]).T hull_indices = find_column_indices(A, hull_vertices) return hull_indices
Нормализованная версия алгоритма будет работать намного быстрее.