Вершины выпуклой оболочки 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 2

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

Нормализованная версия алгоритма будет работать намного быстрее.