From 94322006701c3941b0ddab3a6ee0ae3b6086a0f5 Mon Sep 17 00:00:00 2001 From: Tom Schierenbeck Date: Fri, 6 Sep 2024 11:12:10 +0200 Subject: [PATCH] Docs for polytopes. --- docs/_toc.yml | 1 + docs/polytope.md | 68 +++++++++++++++++++++++++++++++++++ requirements-dev.txt | 3 +- src/random_events/polytope.py | 15 ++++++-- test/test_polytope.py | 10 +++--- 5 files changed, 89 insertions(+), 8 deletions(-) create mode 100644 docs/polytope.md diff --git a/docs/_toc.yml b/docs/_toc.yml index e318bc2..4ac5968 100644 --- a/docs/_toc.yml +++ b/docs/_toc.yml @@ -8,4 +8,5 @@ chapters: - file: conceptual_guide - file: technical_guide - file: advanced_use +- file: polytope - file: autoapi/index \ No newline at end of file diff --git a/docs/polytope.md b/docs/polytope.md new file mode 100644 index 0000000..c4958c1 --- /dev/null +++ b/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() +``` diff --git a/requirements-dev.txt b/requirements-dev.txt index 16fe816..46f9ae3 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -4,4 +4,5 @@ jupyter pillow polytope ortools -scipy \ No newline at end of file +scipy +cvxopt \ No newline at end of file diff --git a/src/random_events/polytope.py b/src/random_events/polytope.py index 7c80b21..fd9a03c 100644 --- a/src/random_events/polytope.py +++ b/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 diff --git a/test/test_polytope.py b/test/test_polytope.py index efdd251..b2b585a 100644 --- a/test/test_polytope.py +++ b/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__':