added chebyshev ball, updated tests and added examples

docs/source/examples/index.rst

@@ -38,3 +38,4 @@ A few generic examples
     Find a half-plane representation of a point cloud <hspace_from_vert>
     Find a vertex representation of a set of half-planes <vert_from_hspace>
     Polytope manipulation example: Minkowski sum and intersection <manipulating_polytopes>
+    Polytope manipulation example: Inner approximation <manipulating_polytopes_approx>

docs/source/examples/manipulating_polytopes.rst

@@ -5,7 +5,7 @@ This example shows how to manipulate polytopes using the ``pycapacity`` package.
 In particular, it shows how to perform Minkowski sum and intersection of polytopes using the ``Polytope`` object's operators ``+`` and ``&`` respectively.
 
 .. code-block:: python
-    
+
     from pycapacity.objects import Polytope # import polytope object
     import numpy as np
 

docs/source/examples/manipulating_polytopes_approx.rst

+Polytope manipulation example: Inner sphere approximation
+=============================================================
+
+This example shows how to manipulate polytopes using the ``pycapacity`` package. 
+In particular, it shows how to calculate the inner sphere approximation of a polytope, using the Chebyshev ball algorithm.
+
+.. code-block:: python
+
+    from pycapacity.objects import Polytope # import polytope object
+    import numpy as np
+
+    # this seed is used to generate the same image 
+    # as in the examples in the docs 
+    np.random.seed(12345)
+
+    N = 10 # ten vertices
+    m = 3   # space dimension
+
+    points = np.array(np.random.rand(m,N))*2-1 # points
+
+    # create a polytope object
+    p = Polytope()
+    p.find_from_point_cloud(points)
+
+    # calculate the inner sphere approximation 
+    s = p.chebyshev_ball()
+
+    # plotting the polytope
+    import matplotlib.pyplot as plt
+    from pycapacity.visual import * # pycapacity visualisation tools
+
+    fig = plt.figure(4)
+    # draw polytopes
+    plot_polytope(plot=fig, 
+        polytope=p, 
+        label='polytope', 
+        face_color="blue", 
+        edge_color='black',  
+        alpha=0.2)
+    plot_ellipsoid(plot=fig, ellipsoid=s, label='sphere')
+    plt.legend()
+    plt.show()
+
+The output of the code is a visualisation of the the polytope and the inner sphere approximation.
+
+.. image:: ../images/polyman_ball.png
+

pycapacity/algorithms.py

@@ -528,18 +528,9 @@ def chebyshev_center(A,b):
     Returns:
         center(array): returns a chebyshev center of the polytope
     """
-    # calculate the vertices
-    Ab_mat = np.hstack((np.array(A),-np.array(b)))
-
-    # calculating chebyshev center
-    norm_vector = np.reshape(np.linalg.norm(Ab_mat[:, :-1], axis=1), (A.shape[0], 1))
-    c = np.zeros((Ab_mat.shape[1],))
-    c[-1] = -1
-    G = matrix(np.hstack((Ab_mat[:, :-1], norm_vector)))
-    h = matrix(- Ab_mat[:, -1:])
-    solvers_opt={'tm_lim': 100000, 'msg_lev': 'GLP_MSG_OFF', 'it_lim':10000}
-    res = cvxopt.glpk.lp(c=c,  G=G, h=h, options=solvers_opt)
-    return np.array(res[1][:-1]).reshape((-1,))
+    # calculate the chebyshev ball
+    c, r = chebyshev_ball(A,b)
+    return c
 
 def chebyshev_ball(A,b):
     """
@@ -618,7 +609,7 @@ def vertex_to_hspace(vertex):
         d(list): vector of half-space representation `Hx<d`
     """
     hull = ConvexHull(vertex.T, qhull_options='QJ')
-    return  hull.equations[:,:-1], -hull.equations[:,-1]
+    return  hull.equations[:,:-1], -hull.equations[:,-1].reshape((-1,1))
 
 
 

pycapacity/objects.py

@@ -92,8 +92,12 @@ class Polytope:
         self.vertices = vertices
         self.faces = faces
         self.face_indices = face_indices
-        self.H = H
-        self.d = d
+        if H is not None and d is not None:
+            self.H = H
+            self.d = d.reshape(-1,1)
+        else:
+            self.H = None
+            self.d = None
     
     def find_vertices(self):
         """
@@ -120,14 +124,15 @@ class Polytope:
         """
         A function calculating the face representation of the polytope from the vertex representation
         """
-        if self.vertices is not None:
-            if self.face_indices is not None:
-                self.faces = face_index_to_vertex(self.vertices,self.face_indices)
-            else:
-                self.face_indices = vertex_to_faces(self.vertices)
-                self.faces = face_index_to_vertex(self.vertices,self.face_indices)
+        if self.vertices is None:
+            print("No vertex representation of the polytope is available - calculating it")
+            self.find_vertices()
+            
+        if self.face_indices is not None:
+            self.faces = face_index_to_vertex(self.vertices,self.face_indices)
         else:
-            print("No vertex representation of the polytope is available")
+            self.face_indices = vertex_to_faces(self.vertices)
+            self.faces = face_index_to_vertex(self.vertices,self.face_indices)
 
     def find_from_point_cloud(self, points):
         """

pycapacity/visual.py

@@ -19,7 +19,7 @@ from matplotlib.patches import Ellipse
 import numpy as np
 from pycapacity.objects import *
 
-def plot_polytope(polytope, plot=None, face_color=None, edge_color=None, vertex_color='black', alpha=None,label=None, center=None, scale=1.0, show_vertices=True, wireframe=False):
+def plot_polytope(polytope, plot=None, color=None, vertex_color='black', face_color=None, edge_color=None, alpha=None,label=None, center=None, scale=1.0, show_vertices=True, wireframe=False):
     """
     A polytope plotting function in 2d and 3d. It plots the polytope faces and vertices from the polytope object.
 
@@ -44,6 +44,18 @@ def plot_polytope(polytope, plot=None, face_color=None, edge_color=None, vertex_
         print("no polytope provided")
         return plot
     
+    # if color is one value, use it for all
+    if color is not None and (isinstance(color, str) or len(color) == 1):
+        face_color=color, 
+        edge_color=color, 
+        vertex_color=color
+    # if color is a list of 3 values, use it for face, edge and vertex
+    elif color is not None and len(color) == 3:
+        face_color, edge_color, vertex_color = color
+
+    if vertex_color is None:
+        show_vertices = False
+
     if label is None:
         label = ''
     if show_vertices:

tests/testing_objects.py

@@ -18,7 +18,7 @@ def test_polytope_from_point_cloud():
     p = Polytope()
     p.find_from_point_cloud(points)
     print(points, p.vertices, np.isin(p.vertices, points))
-    assert p.vertices.shape[1] <= points.shape[1] and np.all(p.H@points - p.d[:,None] < 0.001) and np.all(np.isin(np.round(p.vertices,2), points))
+    assert p.vertices.shape[1] <= points.shape[1] and np.all(p.H@points - p.d.reshape(-1,1)< 0.001) and np.all(np.isin(np.round(p.vertices,2), points))
 
 # write a test for polytope half-plane representation from a set of vertices
 def test_polytope_from_halfplane():
@@ -26,7 +26,7 @@ def test_polytope_from_halfplane():
     p = Polytope()
     p.find_from_point_cloud(points)
     p.find_halfplanes()
-    assert np.all(p.H@points - p.d[:,None] < 0.001)
+    assert np.all(p.H@points - p.d.reshape(-1,1)< 0.001)
 
 # write a test for polytope face representation from a set of vertices
 def test_polytope_from_face():
@@ -60,7 +60,7 @@ def test_polytope_minkowski_sum():
     p1 = Polytope(H = np.vstack((np.eye(3),-np.eye(3))), d = np.vstack((np.ones((3,1)),np.ones((3,1)))))
     p2 = Polytope(H = np.vstack((np.eye(3),-np.eye(3))), d = np.vstack((2*np.ones((3,1)),2*np.ones((3,1)))))
     p_sum = p1 + p2
-    assert np.all(p_sum.H@p1.vertices - p_sum.d[:,None]  < 1e-5) and  np.all(p_sum.H@p2.vertices - p_sum.d[:,None]  < 1e-5)
+    assert np.all(p_sum.H@p1.vertices - p_sum.d.reshape(-1,1)  < 1e-5) and  np.all(p_sum.H@p2.vertices - p_sum.d.reshape(-1,1) < 1e-5)
 
 # test intersection of two polytopes
 # testing for a cube
@@ -71,7 +71,7 @@ def test_polytope_intersection():
     p1.find_halfplanes()
     p2.find_halfplanes()
     p_int.find_vertices()
-    assert np.all(p1.H@p_int.vertices - p1.d[:,None]  < 1e-5) and  np.all(p2.H@p_int.vertices - p2.d[:,None]  < 1e-5)
+    assert np.all(p1.H@p_int.vertices - p1.d.reshape(-1,1) < 1e-5) and  np.all(p2.H@p_int.vertices - p2.d.reshape(-1,1) < 1e-5)
 
 # test chebychev ball of a polytope using a cube
 def chebyshev_ball():