From c19f53e11e6d14a447e315fc8883967f91082e51 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?V=C3=ADtor=20Fr=C3=B3is?= <46361092+vitorfrois@users.noreply.github.com> Date: Fri, 30 Aug 2024 11:45:37 -0300 Subject: [PATCH] Voronoi Tesselation based Discrete Space (#2084) This feature allows the user to build a discrete space based on a random sample of points, where neighbors are defined by Delaunay Triangulation. More specifically, Delaunay Triangulation is a dual-graph representation of the Voronoi Tesselation. Using this algorithm, we can easily find nearest neighbors without delimiting cells edges. --- mesa/experimental/cell_space/__init__.py | 2 + mesa/experimental/cell_space/voronoi.py | 264 ++++++++++++++++++++ mesa/visualization/components/matplotlib.py | 53 ++++ tests/test_cell_space.py | 23 ++ tests/test_solara_viz.py | 11 + 5 files changed, 353 insertions(+) create mode 100644 mesa/experimental/cell_space/voronoi.py diff --git a/mesa/experimental/cell_space/__init__.py b/mesa/experimental/cell_space/__init__.py index dce296aebce..8db71025616 100644 --- a/mesa/experimental/cell_space/__init__.py +++ b/mesa/experimental/cell_space/__init__.py @@ -9,6 +9,7 @@ OrthogonalVonNeumannGrid, ) from mesa.experimental.cell_space.network import Network +from mesa.experimental.cell_space.voronoi import VoronoiGrid __all__ = [ "CellCollection", @@ -20,4 +21,5 @@ "OrthogonalMooreGrid", "OrthogonalVonNeumannGrid", "Network", + "VoronoiGrid", ] diff --git a/mesa/experimental/cell_space/voronoi.py b/mesa/experimental/cell_space/voronoi.py new file mode 100644 index 00000000000..4395ca4ead8 --- /dev/null +++ b/mesa/experimental/cell_space/voronoi.py @@ -0,0 +1,264 @@ +from collections.abc import Sequence +from itertools import combinations +from random import Random + +import numpy as np + +from mesa.experimental.cell_space.cell import Cell +from mesa.experimental.cell_space.discrete_space import DiscreteSpace + + +class Delaunay: + """ + Class to compute a Delaunay triangulation in 2D + ref: http://github.com/jmespadero/pyDelaunay2D + """ + + def __init__(self, center: tuple = (0, 0), radius: int = 9999) -> None: + """ + Init and create a new frame to contain the triangulation + center: Optional position for the center of the frame. Default (0,0) + radius: Optional distance from corners to the center. + """ + center = np.asarray(center) + # Create coordinates for the corners of the frame + self.coords = [ + center + radius * np.array((-1, -1)), + center + radius * np.array((+1, -1)), + center + radius * np.array((+1, +1)), + center + radius * np.array((-1, +1)), + ] + + # Create two dicts to store triangle neighbours and circumcircles. + self.triangles = {} + self.circles = {} + + # Create two CCW triangles for the frame + triangle1 = (0, 1, 3) + triangle2 = (2, 3, 1) + self.triangles[triangle1] = [triangle2, None, None] + self.triangles[triangle2] = [triangle1, None, None] + + # Compute circumcenters and circumradius for each triangle + for t in self.triangles: + self.circles[t] = self._circumcenter(t) + + def _circumcenter(self, triangle: list) -> tuple: + """ + Compute circumcenter and circumradius of a triangle in 2D. + """ + points = np.asarray([self.coords[v] for v in triangle]) + points2 = np.dot(points, points.T) + a = np.bmat([[2 * points2, [[1], [1], [1]]], [[[1, 1, 1, 0]]]]) + + b = np.hstack((np.sum(points * points, axis=1), [1])) + x = np.linalg.solve(a, b) + bary_coords = x[:-1] + center = np.dot(bary_coords, points) + + radius = np.sum(np.square(points[0] - center)) # squared distance + return (center, radius) + + def _in_circle(self, triangle: list, point: list) -> bool: + """ + Check if point p is inside of precomputed circumcircle of triangle. + """ + center, radius = self.circles[triangle] + return np.sum(np.square(center - point)) <= radius + + def add_point(self, point: Sequence) -> None: + """ + Add a point to the current DT, and refine it using Bowyer-Watson. + """ + point_index = len(self.coords) + self.coords.append(np.asarray(point)) + + bad_triangles = [] + for triangle in self.triangles: + if self._in_circle(triangle, point): + bad_triangles.append(triangle) + + boundary = [] + triangle = bad_triangles[0] + edge = 0 + + while True: + opposite_triangle = self.triangles[triangle][edge] + if opposite_triangle not in bad_triangles: + boundary.append( + ( + triangle[(edge + 1) % 3], + triangle[(edge - 1) % 3], + opposite_triangle, + ) + ) + edge = (edge + 1) % 3 + if boundary[0][0] == boundary[-1][1]: + break + else: + edge = (self.triangles[opposite_triangle].index(triangle) + 1) % 3 + triangle = opposite_triangle + + for triangle in bad_triangles: + del self.triangles[triangle] + del self.circles[triangle] + + new_triangles = [] + for e0, e1, opposite_triangle in boundary: + triangle = (point_index, e0, e1) + self.circles[triangle] = self._circumcenter(triangle) + self.triangles[triangle] = [opposite_triangle, None, None] + if opposite_triangle: + for i, neighbor in enumerate(self.triangles[opposite_triangle]): + if neighbor and e1 in neighbor and e0 in neighbor: + self.triangles[opposite_triangle][i] = triangle + + new_triangles.append(triangle) + + n = len(new_triangles) + for i, triangle in enumerate(new_triangles): + self.triangles[triangle][1] = new_triangles[(i + 1) % n] # next + self.triangles[triangle][2] = new_triangles[(i - 1) % n] # previous + + def export_triangles(self) -> list: + """ + Export the current list of Delaunay triangles + """ + triangles_list = [ + (a - 4, b - 4, c - 4) + for (a, b, c) in self.triangles + if a > 3 and b > 3 and c > 3 + ] + return triangles_list + + def export_voronoi_regions(self): + """ + Export coordinates and regions of Voronoi diagram as indexed data. + """ + use_vertex = {i: [] for i in range(len(self.coords))} + vor_coors = [] + index = {} + for triangle_index, (a, b, c) in enumerate(sorted(self.triangles)): + vor_coors.append(self.circles[(a, b, c)][0]) + use_vertex[a] += [(b, c, a)] + use_vertex[b] += [(c, a, b)] + use_vertex[c] += [(a, b, c)] + + index[(a, b, c)] = triangle_index + index[(c, a, b)] = triangle_index + index[(b, c, a)] = triangle_index + + regions = {} + for i in range(4, len(self.coords)): + vertex = use_vertex[i][0][0] + region = [] + for _ in range(len(use_vertex[i])): + triangle = next( + triangle for triangle in use_vertex[i] if triangle[0] == vertex + ) + region.append(index[triangle]) + vertex = triangle[1] + regions[i - 4] = region + + return vor_coors, regions + + +def round_float(x: float) -> int: + return int(x * 500) + + +class VoronoiGrid(DiscreteSpace): + triangulation: Delaunay + voronoi_coordinates: list + regions: list + + def __init__( + self, + centroids_coordinates: Sequence[Sequence[float]], + capacity: float | None = None, + random: Random | None = None, + cell_klass: type[Cell] = Cell, + capacity_function: callable = round_float, + cell_coloring_property: str | None = None, + ) -> None: + """ + A Voronoi Tessellation Grid. + + Given a set of points, this class creates a grid where a cell is centered in each point, + its neighbors are given by Voronoi Tessellation cells neighbors + and the capacity by the polygon area. + + Args: + centroids_coordinates: coordinates of centroids to build the tessellation space + capacity (int) : capacity of the cells in the discrete space + random (Random): random number generator + CellKlass (type[Cell]): type of cell class + capacity_function (Callable): function to compute (int) capacity according to (float) area + cell_coloring_property (str): voronoi visualization polygon fill property + """ + super().__init__(capacity=capacity, random=random, cell_klass=cell_klass) + self.centroids_coordinates = centroids_coordinates + self._validate_parameters() + + self._cells = { + i: cell_klass(self.centroids_coordinates[i], capacity, random=self.random) + for i in range(len(self.centroids_coordinates)) + } + + self.regions = None + self.triangulation = None + self.voronoi_coordinates = None + self.capacity_function = capacity_function + self.cell_coloring_property = cell_coloring_property + + self._connect_cells() + self._build_cell_polygons() + + def _connect_cells(self) -> None: + """ + Connect cells to neighbors based on given centroids and using Delaunay Triangulation + """ + self.triangulation = Delaunay() + for centroid in self.centroids_coordinates: + self.triangulation.add_point(centroid) + + for point in self.triangulation.export_triangles(): + for i, j in combinations(point, 2): + self._cells[i].connect(self._cells[j]) + self._cells[j].connect(self._cells[i]) + + def _validate_parameters(self) -> None: + if self.capacity is not None and not isinstance(self.capacity, float | int): + raise ValueError("Capacity must be a number or None.") + if not isinstance(self.centroids_coordinates, Sequence) or not isinstance( + self.centroids_coordinates[0], Sequence + ): + raise ValueError("Centroids should be a list of lists") + dimension_1 = len(self.centroids_coordinates[0]) + for coordinate in self.centroids_coordinates: + if dimension_1 != len(coordinate): + raise ValueError("Centroid coordinates should be a homogeneous array") + + def _get_voronoi_regions(self) -> tuple: + if self.voronoi_coordinates is None or self.regions is None: + self.voronoi_coordinates, self.regions = ( + self.triangulation.export_voronoi_regions() + ) + return self.voronoi_coordinates, self.regions + + @staticmethod + def _compute_polygon_area(polygon: list) -> float: + polygon = np.array(polygon) + x = polygon[:, 0] + y = polygon[:, 1] + return 0.5 * np.abs(np.dot(x, np.roll(y, 1)) - np.dot(y, np.roll(x, 1))) + + def _build_cell_polygons(self): + coordinates, regions = self._get_voronoi_regions() + for region in regions: + polygon = [coordinates[i] for i in regions[region]] + self._cells[region].properties["polygon"] = polygon + polygon_area = self._compute_polygon_area(polygon) + self._cells[region].properties["area"] = polygon_area + self._cells[region].capacity = self.capacity_function(polygon_area) + self._cells[region].properties[self.cell_coloring_property] = 0 diff --git a/mesa/visualization/components/matplotlib.py b/mesa/visualization/components/matplotlib.py index 83d0e3d8eaf..c35fe97bc5e 100644 --- a/mesa/visualization/components/matplotlib.py +++ b/mesa/visualization/components/matplotlib.py @@ -6,6 +6,7 @@ from matplotlib.ticker import MaxNLocator import mesa +from mesa.experimental.cell_space import VoronoiGrid @solara.component @@ -20,6 +21,8 @@ def SpaceMatplotlib(model, agent_portrayal, dependencies: list[any] | None = Non _draw_network_grid(space, space_ax, agent_portrayal) elif isinstance(space, mesa.space.ContinuousSpace): _draw_continuous_space(space, space_ax, agent_portrayal) + elif isinstance(space, VoronoiGrid): + _draw_voronoi(space, space_ax, agent_portrayal) else: _draw_grid(space, space_ax, agent_portrayal) solara.FigureMatplotlib(space_fig, format="png", dependencies=dependencies) @@ -150,6 +153,56 @@ def portray(space): _split_and_scatter(portray(space), space_ax) +def _draw_voronoi(space, space_ax, agent_portrayal): + def portray(g): + x = [] + y = [] + s = [] # size + c = [] # color + + for cell in g.all_cells: + for agent in cell.agents: + data = agent_portrayal(agent) + x.append(cell.coordinate[0]) + y.append(cell.coordinate[1]) + if "size" in data: + s.append(data["size"]) + if "color" in data: + c.append(data["color"]) + out = {"x": x, "y": y} + # This is the default value for the marker size, which auto-scales + # according to the grid area. + out["s"] = s + if len(c) > 0: + out["c"] = c + + return out + + x_list = [i[0] for i in space.centroids_coordinates] + y_list = [i[1] for i in space.centroids_coordinates] + x_max = max(x_list) + x_min = min(x_list) + y_max = max(y_list) + y_min = min(y_list) + + width = x_max - x_min + x_padding = width / 20 + height = y_max - y_min + y_padding = height / 20 + space_ax.set_xlim(x_min - x_padding, x_max + x_padding) + space_ax.set_ylim(y_min - y_padding, y_max + y_padding) + space_ax.scatter(**portray(space)) + + for cell in space.all_cells: + polygon = cell.properties["polygon"] + space_ax.fill( + *zip(*polygon), + alpha=min(1, cell.properties[space.cell_coloring_property]), + c="red", + ) # Plot filled polygon + space_ax.plot(*zip(*polygon), color="black") # Plot polygon edges in red + + @solara.component def PlotMatplotlib(model, measure, dependencies: list[any] | None = None): fig = Figure() diff --git a/tests/test_cell_space.py b/tests/test_cell_space.py index 050e39a09a6..a918f59e53b 100644 --- a/tests/test_cell_space.py +++ b/tests/test_cell_space.py @@ -11,6 +11,7 @@ Network, OrthogonalMooreGrid, OrthogonalVonNeumannGrid, + VoronoiGrid, ) @@ -372,6 +373,28 @@ def test_networkgrid(): assert connection.coordinate in G.neighbors(i) +def test_voronoigrid(): + points = [[0, 1], [1, 3], [1.1, 1], [1, 1]] + + grid = VoronoiGrid(points) + + assert len(grid._cells) == len(points) + + # Check cell neighborhood + assert len(grid._cells[0]._connections) == 2 + for connection in grid._cells[0]._connections: + assert connection.coordinate in [[1, 1], [1, 3]] + + with pytest.raises(ValueError): + VoronoiGrid(points, capacity="str") + + with pytest.raises(ValueError): + VoronoiGrid((1, 1)) + + with pytest.raises(ValueError): + VoronoiGrid([[0, 1], [0, 1, 1]]) + + def test_empties_space(): import networkx as nx diff --git a/tests/test_solara_viz.py b/tests/test_solara_viz.py index 5fef7e0b214..660277b3d7e 100644 --- a/tests/test_solara_viz.py +++ b/tests/test_solara_viz.py @@ -140,6 +140,17 @@ def test_call_space_drawer(mocker): model, agent_portrayal, dependencies=dependencies ) + # check voronoi space drawer + voronoi_model = mesa.Model() + voronoi_model.grid = mesa.experimental.cell_space.VoronoiGrid( + centroids_coordinates=[(0, 1), (0, 0), (1, 0)], + ) + solara.render( + SolaraViz( + model_class=voronoi_model, model_params={}, agent_portrayal=agent_portrayal + ) + ) + def test_slider(): slider_float = Slider("Agent density", 0.8, 0.1, 1.0, 0.1)