Docs for polytopes.

docs/_toc.yml

@@ -8,4 +8,5 @@ chapters:
 - file: conceptual_guide
 - file: technical_guide
 - file: advanced_use
+- file: polytope
 - file: autoapi/index

docs/polytope.md

@@ -0,0 +1,68 @@
+---
+jupytext:
+  cell_metadata_filter: -all
+  formats: md:myst
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.11.5
+kernelspec:
+  display_name: Python 3
+  language: python
+  name: python3
+---
+
+
+# Polytopes
+
+The tutorial walks through the conversion from polytopes to events from the product algebra.
+
+First, let's define a polytope in 2D.
+We start by generating some random points and rotating them. Then we create a polytope over the set of points.
+
+```{code-cell} ipython3
+import numpy as np
+from scipy.spatial import ConvexHull
+from random_events.polytope import Polytope
+import plotly
+plotly.offline.init_notebook_mode()
+import plotly.graph_objects as go
+
+def rotation_matrix(theta: float):
+    return np.array([[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]])
+
+np.random.seed(69)
+points = np.random.uniform(-1, 1, (100, 2))
+points = points @ rotation_matrix(0.3)
+
+polytope = Polytope.from_2d_points(points)
+```
+
+Let's have a look at the generated data and polytope.
+
+```{code-cell} ipython3
+fig = go.Figure()
+fig.add_trace(go.Scatter(x=points[:, 0], y=points[:, 1], mode='markers', name='Points'))
+convex_hull = ConvexHull(points)
+hull_points = np.vstack([points[convex_hull.vertices], points[convex_hull.vertices[0]]])
+fig.add_trace(go.Scatter(x=hull_points[:, 0], y=hull_points[:, 1], mode='lines', name='Convex Hull'))
+fig.show()
+```
+
+Next, let's create the inner box approximation of the polytope.
+
+```{code-cell} ipython3
+inner_event = polytope.inner_box_approximation()
+fig.add_traces(inner_event.plot("aquamarine"))
+fig.show()
+```
+
+Let's create the outer box approximation of the polytope.
+
+```{code-cell} ipython3
+outer_event = polytope.outer_box_approximation()
+outer_event = outer_event - inner_event
+fig.add_traces(outer_event.plot("aqua"))
+fig.show()
+```

requirements-dev.txt

@@ -4,4 +4,5 @@ jupyter
 pillow
 polytope
 ortools
-scipy
+scipy
+cvxopt

src/random_events/polytope.py

@@ -11,6 +11,12 @@
 
 
 class Polytope(polytope.Polytope):
+    """
+    Extension of the polytope class from the polytope library.
+
+    This class enables conversion to simple events and provides the inner box and outer box approximation
+    from https://cse.lab.imtlucca.it/~bemporad/publications/papers/compgeom-boxes.pdf.
+    """
 
     @classmethod
     def from_polytope(cls, polytope_: polytope.Polytope) -> Self:
@@ -79,7 +85,7 @@ def inner_box_approximation(self, minimum_volume: float = 0.1) -> Event:
 
     def as_box_polytope(self) -> Self:
         """
-        Convert the polytope to a box polytope.
+        :return: The polytope as box polytope.
         """
         lower, upper = self.bounding_box
         return self.from_box([(a[0], b[0]) for a, b in zip(lower, upper)])
@@ -93,6 +99,8 @@ def split_on_axis_value(self, axis: int, value: np.ndarray) -> Tuple[Self, Self]
 
         :param axis: The axis to split on.
         :param value: The value to split on.
+
+        :return: The left and right split of the polytope.
         """
         a_vector = np.zeros((1, self.A.shape[1]))
         a_vector[0, axis] = 1.
@@ -113,7 +121,10 @@ def split_on_axis_value(self, axis: int, value: np.ndarray) -> Tuple[Self, Self]
     def outer_box_approximation(self, minimum_volume: float = 0.1) -> Event:
         """
         Compute an outer box approximation of the polytope.
+        This implements Algorithm 6 in https://cse.lab.imtlucca.it/~bemporad/publications/papers/compgeom-boxes.pdf
+
         :param minimum_volume: The minimum volume (epsilon) for the approximation.
+
         :return: The outer box approximation of the polytope as a random event.
         """
         polytopes_to_split = deque([self])
@@ -125,7 +136,7 @@ def outer_box_approximation(self, minimum_volume: float = 0.1) -> Event:
 
             # if the box is too small, skip
             volume = bounding_box_of_current_polytope.volume
-            if volume < minimum_volume:
+            if volume < minimum_volume or bounding_box_of_current_polytope <= current_polytope:
                 resulting_boxes.append(current_polytope)
                 continue
 

test/test_polytope.py

@@ -42,11 +42,11 @@ def test_outer_box_approximation(self):
         polytope = Polytope.from_2d_points(self.points)
         result = polytope.outer_box_approximation(0.1)
         self.assertTrue(result.is_disjoint())
-        fig = go.Figure()
-        fig.add_trace(go.Scatter(x=self.points[:, 0], y=self.points[:, 1], mode='markers', name='points'))
-        fig.add_traces(result.plot())
-        fig.update_layout(result.plotly_layout())
-        fig.show()
+        # fig = go.Figure()
+        # fig.add_trace(go.Scatter(x=self.points[:, 0], y=self.points[:, 1], mode='markers', name='points'))
+        # fig.add_traces(result.plot())
+        # fig.update_layout(result.plotly_layout())
+        # fig.show()
 
 
 if __name__ == '__main__':