chebyshev ball added

pycapacity/algorithms.py

@@ -541,6 +541,34 @@ def chebyshev_center(A,b):
     res = cvxopt.glpk.lp(c=c,  G=G, h=h, options=solvers_opt)
     return np.array(res[1][:-1]).reshape((-1,))
 
+def chebyshev_ball(A,b):
+    """
+    Calculating chebyshev ball of a polytope given the half-space representation
+
+    https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.HalfspaceIntersection.html#r9b902253b317-1
+
+    Args:
+        A(list):  
+            matrix of half-space representation `Ax<b`
+        b(list): 
+            vector of half-space representation `Ax<b`
+    Returns:
+        center(array): returns a chebyshev center of the polytope
+        radius(float): returns a chebyshev radius 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,)), np.array(res[1][-1]).reshape((-1,))
+
 def hspace_to_vertex(H,d):
     """
     From half-space representation to the vertex representation

pycapacity/objects.py

@@ -166,6 +166,23 @@ class Polytope:
         d_int = np.vstack((self.d,p.d))
         return Polytope(H=H_int,d=d_int)
 
+    # implement chebyshev ball
+    def chebyshev_ball(self):
+        """
+        Function calculating the Chebyshev ball of the polytope, the largest ball that can be inscribed in the polytope. 
+        The function calculates the center and radius of the ball and returns an ellipsoid object representing the Chebyshev ball of the polytope
+
+        Th function uses one Linear Programming problem to find the Chebyshev ball of the polytope.
+        See also: https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.HalfspaceIntersection.html#r9b902253b317-1
+
+        Returns:
+            Ellipsoid(Ellipsoid): an ellipsoid object representing the Chebyshev ball of the polytope
+        
+        """
+        if self.H is None or self.d is None:
+            self.find_halfplanes()
+        center, radius = chebyshev_ball(self.H,self.d.reshape(-1,1))
+        return Ellipsoid(radii=np.ones(self.H.shape[1],)*radius,rotation=np.eye(self.H.shape[1]), center=center)
 
 
 class Ellipsoid:
@@ -184,6 +201,7 @@ class Ellipsoid:
     :ivar rotation: rotation of the ellipsoid SO3 matrix
     """
 
-    def __init__(self, radii=None, rotation=None):
+    def __init__(self, radii=None, rotation=None, center=None):
         self.radii = radii
-        self.rotation = rotation
\ No newline at end of file
+        self.rotation = rotation
+        self.center = center
\ No newline at end of file

pycapacity/visual.py

@@ -359,6 +359,7 @@ def plot_ellipsoid(radii=None, rotation=None, ellipsoid=None, center=None, plot=
     if ellipsoid is not None:
         radii = ellipsoid.radii
         rotation = ellipsoid.rotation
+        center = ellipsoid.center
 
     if ellipsoid is None and radii is None:
         print("no data to plot")

tests/testing_objects.py

@@ -71,4 +71,10 @@ 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)
\ No newline at end of file
+    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)
+
+# test chebychev ball of a polytope using a cube
+def chebyshev_ball():
+    p = Polytope(H = np.vstack((np.eye(3),-np.eye(3))), d = np.vstack((np.zeros((3,1)),np.ones((3,1)))))
+    e = p.chebyshev_ball()
+    assert np.all(e.radii == 0.5) and np.all(e.rotation == np.eye(3)) and np.all(e.center == 0.5*np.ones((3,1)))
\ No newline at end of file