Repository: projectmesa/mesa
Pull request #2084 (2024-03-21 00:41:22)

Voronoi Tesselation based Discrete Space
## First version of Voronoi Tesselation based Discrete Space, 
described on issue #1895 . 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.
![image](https://github.com/projectmesa/mesa/assets/46361092/a966046b-a6b4-4680-b396-c22d2db2064d)


The library chosen for the triangulation algorithm implementation was [PyHull](https://github.com/materialsvirtuallab/pyhull), a wrapper of qhull C library, which deals with spatial calculations. The choice was made considering performance, size and usability (SciPy has a similar module but much heavier).   


Based on my discussion with @EwoutH on the issue thread, i thought it would be useful to inherit DiscreteSpace class. 
``` python
class VoronoiGrid(DiscreteSpace):
```

### Example
``` python
from mesa import Model
from mesa.experimental.cell_space import VoronoiGrid, CellAgent

class MyAgent(CellAgent):
    def __init__(self, unique_id, model):
        super().__init__(unique_id, model)
```
Users can input cells centroids
``` python
points = [
    [0, 0], [0, 1], [0, 2], 
    [1, 0], [1, 1], [1, 2],
    [2, 0], [2, 1], [2, 2]]
grid = VoronoiGrid(centroids_coordinates=points, capacity=1, random=model.random)
```
and cells neighborhoods are defined by Delaunay triangulation

``` python
agent = MyAgent(1, model)
cell = grid.select_random_empty_cell()
agent.move_to(cell)
print(cell)
print(cell.neighborhood())
```
```
Cell([1, 0], [<__main__.MyAgent object at 0x7f2177268450>])
CellCollection({Cell([0, 1], []): [], Cell([1, 1], []): [], Cell([0, 0], []): [], Cell([2, 1], []): [], Cell([2, 0], []): []})
```
### Next Steps
Definitely, next steps involve computing the boundaries of each cell using Voronoi Tesselation. This allow us to define capacity of each cell based on its area/volume.

EwoutH: Thanks, it's a great start! Using DelaunayTri is a smart approach.

I guess `move_to` currently move an agent to the cells centroid? While that is logical behavior, I would also like to see an method in which agents are moved to A) the closest point in the cell from their current position and B) a random spot in the cell. This could be a separate PR.

Another thing that it's needed in the long term is a method to check in which cell a point (and thus an agent) is.

@wang-boyu I would also appreciate a review from your mesa-geo expertise.

For the currently functionality, I would like some additional docstring and unittests. Then you can decide if you want to add further functionality to this PR, or merge this first and follow up with new PRs.
wang-boyu: Thanks for the PR. I have to say that I'm still a bit confused about why this is based on discrete space, not continuous space (as I mentioned in https://github.com/projectmesa/mesa/issues/1895#issuecomment-2016514829 and https://github.com/projectmesa/mesa/issues/1895#issuecomment-2016521840).

As a concrete example, consider this simple space:

| Cell 1 | Cell 2 | Cell 3 |
|:---:|:---:|:---:|
| O    | ?   | X  |

If it is a Voronoi space, should Cell 2 be split into halves, where left half belongs to agent of type O and right half belongs to agent of type X? How does this work in discrete space?
tpike3: Apologies, a little late to this discussion, this is really cool @vitorfrois. 

As dependencies and keeping Mesa lightweight is a constant struggle, but form a very superficial look is there a reason we wont use networkx  (https://networkx.org/documentation/stable/reference/algorithms/voronoi.html) and this has a corresponding example that may work for mesa_geo (https://networkx.org/documentation/stable/auto_examples/geospatial/plot_delaunay.html)? 

Let me know what you think/ what I am missing? 
vitorfrois: @wang-boyu 
> I have to say that I'm still a bit confused about why this is based on discrete space, not continuous space. How does this work in discrete space?

I'm considering the discrete space for Voronoi Tesselation as a **static/non changing graph**, where neighbors are obtained by **Delaunay Triangulation**. As @EwoutH said, it is only a special case for HexGrid where cell are not equally spaced. 

I think its more about a choice of implementation where I considered mainly the Cholera example. As I suggested yesterday, 
> since the result of Delaunay Triangulation is a graph, we can use NetworkGrid to represent it, and in the continuous case, implement in the Mesa Geo using Shapely.

Does that make sense?




wang-boyu: > I'm considering the discrete space for Voronoi Tesselation as a **static/non changing graph**, where neighbors are obtained by **Delaunay Triangulation**.

Thanks for clarifying. Appreciate it.

Question: Is the **Delaunay Triangulation** done in continuous space?
vitorfrois: > Question: Is the **Delaunay Triangulation** done in continuous space?

Yes, it is
vitorfrois: **Updates**
- Implemented Voronoi from scratch (basically copied a simple already implemented algorithm)
- Now using the regions resulting from Delaunay/Voronoi Triangulation to draw Voronoi cells in JupyterViz
- Altough, could not cover `visualization/components/matplotlib.py` with tests :(
vitorfrois: Omg, sorry I did not notice the license. I copied like 70% of it's
implementation, changing variables, function names and other minor things.

On Sun, Jul 21, 2024, 08:44 rht ***@***.***> wrote:

> ***@***.**** commented on this pull request.
> ------------------------------
>
> In mesa/experimental/cell_space/voronoi.py
> <https://github.com/projectmesa/mesa/pull/2084#discussion_r1685719329>:
>
> > @@ -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
>
> How much percentage of the code is based on this GitHub repo? Their
> license is GPLv3, and we actually can't copy their code because there will
> be a licensing issue.
>
> —
> Reply to this email directly, view it on GitHub
> <https://github.com/projectmesa/mesa/pull/2084#pullrequestreview-2190309373>,
> or unsubscribe
> <https://github.com/notifications/unsubscribe-auth/ALBWUBC63DYQCEJC5TPSGOLZNONJHAVCNFSM6AAAAABFAPRYMCVHI2DSMVQWIX3LMV43YUDVNRWFEZLROVSXG5CSMV3GSZLXHMZDCOJQGMYDSMZXGM>
> .
> You are receiving this because you were mentioned.Message ID:
> ***@***.***>
>

rht: There are 2 ways to proceed:
1. Use ChatGPT to do code laundry, i.e. ask it to do Delaunay triangulation code. It probably is going to be based on that GPLv3 repo, but you can blame ClosedAI for that instead.
2. Use SciPy, but have SciPy to be an optional dependency, only imported when the user needs `VoronoiGrid`

I'm more in favor of option 2, because option 1 is rather shady.
EwoutH: Yeah the license issue is a big one. Another (far fetched) option is to contact the maintainer to see if he’s willing to publish on another, more permissive license.
EwoutH: @vitorfrois I would like to try to contact the maintainers. Do you want to write them a message or would you like me to do it?
EwoutH: I can do it tomorrow, then hopefully we can move this forward next week!

Patch:
diff --git a/mesa/experimental/cell_space/__init__.py b/mesa/experimental/cell_space/__init__.py
index dce296ae..8db71025 100644
--- a/mesa/experimental/cell_space/__init__.py
+++ b/mesa/experimental/cell_space/__init__.py
@@ -9,6 +9,7 @@ from mesa.experimental.cell_space.grid import (
     OrthogonalVonNeumannGrid,
 )
 from mesa.experimental.cell_space.network import Network
+from mesa.experimental.cell_space.voronoi import VoronoiGrid
 
 __all__ = [
     "CellCollection",
@@ -20,4 +21,5 @@ __all__ = [
     "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 00000000..4395ca4e
--- /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 83d0e3d8..c35fe97b 100644
--- a/mesa/visualization/components/matplotlib.py
+++ b/mesa/visualization/components/matplotlib.py
@@ -6,6 +6,7 @@ from matplotlib.figure import Figure
 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 _draw_continuous_space(space, space_ax, agent_portrayal):
     _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()


Test patch:
diff --git a/tests/test_cell_space.py b/tests/test_cell_space.py
index 050e39a0..a918f59e 100644
--- a/tests/test_cell_space.py
+++ b/tests/test_cell_space.py
@@ -11,6 +11,7 @@ from mesa.experimental.cell_space import (
     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 5fef7e0b..660277b3 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)