# Improve "Game of Life" on a mosaic

I had the idea of generalizing Conways Game of Life from a regular lattice to act on a Poisson Voronoi mosaic (tessellation). The code below works but is maybe not the best idea in terms of speed.

I'm using Mathematica 8, can you help improving the code? Thanks in advance.

SeedRandom;
Points = RandomReal[1, {999, 2}];
VD = VoronoiDiagram[Points];

UnboundedVertices =
Flatten[Position[VD[], _?(Head[#] == Ray &), 1]];
BC = Select[VD[],
Length@Intersection[UnboundedVertices, #[]] == 0 &][[All,
1]]; (* all bounded cells *)

VertexConnection = Table[{}, {i, 1, Length@VD[]}];
For[k = 1, k <= Length[VD[]], k++,
CellVertices = VD[[2, k, 2]];
For[j = 1, j <= Length[CellVertices], j++,
VertexConnection[[CellVertices[[j]]]] =
Append[VertexConnection[[CellVertices[[j]]]], k]
]
];
NC = {};(* neighbor cells for each cell *)
For[k = 1, k <= Length[VD[]], k++,
AppendTo[NC,
Complement[Union@Flatten@VertexConnection[[VD[[2, k, 2]]]], {k}]]
]

(* find two cells in the middle *)
BoundedVertices =
VD[[1, Complement[Range[Length@VD[]], UnboundedVertices]]];
pos = First@
Position[BoundedVertices,
First@Nearest[BoundedVertices, {0.5, 0.5}, 1]];
Start = RandomSample[BC,
300];(*Take[Select[VD[],Length@Intersection[pos,#[]]>0&][[All,\
1]],2];*)

Manipulate[
Refresh[
If[updating || onestep,
t++;
onestep = False;
ActiveCells =
Union[Select[ActiveCells,
MemberQ[{2, 3, 4},
Length@Intersection[ActiveCells, NC[[#]]]] &],
Length@Intersection[ActiveCells, NC[[#]]] == 3 &]];
];
Graphics[
Table[{Black, Polygon[VD[[1, VD[[2, ActiveCells[[i]], 2]]]]]}, {i,
1, Length@ActiveCells}], PlotRange -> {{0, 1}, {0, 1}},
Frame -> True, FrameTicks -> False, ImageSize -> Large],
UpdateInterval -> If[updating, 0, Infinity]
],
{{updating, False, "run simulation"}, {True, False}},
Button["update one step", onestep = True, ImageSize -> Medium],
Button["reset", updating = False; ActiveCells = Start;
DeadCells = Complement[BC, ActiveCells]; t = 0, ImageSize -> Small],
{{t, 0}, ControlType -> None},
{{ActiveCells, Start}, ControlType -> None},
{{DeadCells, Complement[BC, ActiveCells]}, ControlType -> None},
{{onestep, False}, ControlType -> None},
ControlPlacement -> Left
]


I tried different rules, but all patterns are either oscillating or converging to a steady state, no gliders. Maybe worth some research. • Don't forget to add Needs["ComputationalGeometry"]. Jun 8, 2018 at 15:36
• Your problem is interesting and perhaps somebody can help, but it would be much more beneficial if you could at least go through your (long) piece of code and identify bottlenecks for people to work on. Otherwise you are simply asking us to do all your work for you, and that may not be well received. Jun 8, 2018 at 15:39
• Yes, sorry it requires 'ComputationalGeometry'.
– JHT
Jun 8, 2018 at 15:45
• The problem is getting the neighbors for each cell faster and drawing so many polygons.
– JHT
Jun 8, 2018 at 15:46
• Nice idea! Sounds as if this post was highly relevant. In particular Chip Hurst's answer. Jun 8, 2018 at 15:46

The following employs VoronoiMesh, so it won't work with Mathematica 8.0. But it might help you to get an idea how to speed up the computional part.

We will assign an integer to each cell. 1 stands for a living cell, 0 for a dead one.

The following is a CompiledFunction that computes the the state (dead or alive) from the current state, the number of living neighbor cells, and the total number of neighbors (called degree). Note that the function has the RuntimeAttribute Listable so that it will thread over lists. You may your own rules of course and you will quite likely have to do that because generic Voronoi cells will have less than the 8 neighbors in the classical Game of Life.

computeStates =
Compile[{{state, _Integer}, {neighsalive, _Integer}, {degree,  _Integer}},
If[state == 1,
Switch[neighsalive,
2, 1,
3, 1,
4, 1,
_, 0
],
Switch[neighsalive,
3, 1,
_, 0
]
],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True,
RuntimeOptions -> "Speed"
];


Next is a function that computes the adjacency matrix A for the cells. That means, that A[[i,j]] will be one if cells i and j are neighbors to each other and 0 else. (A cell is never its own neighbor.) This matrix will be crucial in counting the number of living neighbors. If state is the vector of ones and zeros describing the states of the cells, the i-th entry in A.state will contain the number of living neighbors. As this involves a matrix-vector product and since this operation is highly optimized, there is hardly a more efficient way to count the number of living neigbors in a general mesh (there may be more performant ways on a rectangular grid, though).

CellCellAdjacencyMatrix[vertexcount_, cells_] :=
Module[{m, B, lens, nn, rp},
m = Length[cells];
B = If[m > 0,
lens = Length /@ cells;
rp = Join[{0}, Accumulate[lens]];
nn = rp[[-1]];
SparseArray @@ {Automatic, {m, vertexcount}, 0, {1, {
rp,
Partition[Flatten[cells], 1]
},
ConstantArray[1, nn]}}
,
{}
];
(# - DiagonalMatrix[Diagonal[#]]) &@Unitize[B.Transpose[B]]
];


# Example 1: Classical Game of Life

This produces a classical Game of Life on a regular $50 \times 50$ grid.

m = 50;
pts = Tuples[Range[m], 2];
R = VoronoiMesh[pts, {{0.5, m + 0.5}, {0.5, m + 0.5}}];


Computing the adjacencymatrix, creating a start vector and computing the degrees of each cell.

A = CellCellAdjacencyMatrix[
MeshCellCount[R, 0],
Join @@ MeshCells[R, 2, "Multicells" -> True][[All, 1]]
];
start = RandomInteger[{0, 1}, Length[A]];
degrees = Total[A];


This runs 5000 generations within one second on my laptop.

history = NestList[computeStates[#, A.#, degrees] &, start, 5000]; // AbsoluteTiming // First


0.983295

Here is a visualization. I am not an expert in interactive computations, so this just shows you the precomputed states.

Manipulate[
Dynamic@HighlightMesh[R,
{2, Flatten[SparseArray[history[[t]]]["NonzeroPositions"]]},
ImageSize -> Medium
],
{{t, 1}, 1, Length[history], 1},
TrackedSymbols :> t
]


Here is just the last frame: # Example 2: Random Voronoi diagram

xrange = {-1, 1};
yrange = {0, 1};
n = 999;
pts = Transpose[{RandomReal[xrange, n], RandomReal[yrange, n]}];
R = VoronoiMesh[pts, {xrange, yrange}]
MeshCellCount[R, 0],
Join @@ MeshCells[R, 2, "Multicells" -> True][[All, 1]]
];
start = RandomInteger[{0, 1}, Length[A]];
degrees = Total[A];

history = NestList[computeStates[#, A.#, degrees] &, start, 12000]; //
AbsoluteTiming // First


Again, the last frame: # Modification for legacy versions of Mathematica

This should also work with older versions.

A plotting routine:

plot[t_] := Module[{alive, dead},
alive = Flatten[SparseArray[history[[t]]]["NonzeroPositions"]];
Graphics[
{EdgeForm[Thin], FaceForm[],
GraphicsComplex[pts,
{
FaceForm[Darker@Green], Polygon[cells[[alive]]],
}
]},
PlotRange -> {xrange, yrange}
]
];


Some preparations...

xrange = {-1, 1};
yrange = {0, 1};
n = 999;
pts = Transpose[{RandomReal[xrange, n], RandomReal[yrange, n]}];
Needs["ComputationalGeometry"];
VD = VoronoiDiagram[pts];

UnboundedVertices = Flatten[Position[VD[], _?(Head[#] == Ray &), 1]];
BC = Select[VD[], Length@Intersection[UnboundedVertices, #[]] == 0 &][[All, 1]];
pts = DeleteCases[VD[], _Ray];
cells = VD[[2, BC, 2]];

start = RandomInteger[{0, 1}, Length[A]];
degrees = Total[A];


... some modified rules that will make the game more lively on random Voronoi diagrams

computeStates = Compile[{{state, _Integer}, {neighsalive, _Integer}, {degree, _Integer}},
If[state == 1,
If[2 <= neighsalive <= Quotient[degree, 2], 1, 0],
If[Quotient[degree, 4] <= neighsalive <= Quotient[degree, 2], 1,
0]
],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True,
RuntimeOptions -> "Speed"
];


... the actual simulation...

history = NestList[computeStates[#, A.#, degrees] &, start, 1000]; // AbsoluteTiming // First


0.085566

... and a visualization

Manipulate[
Dynamic[plot[t]],
{{t, 1}, 1, Length[history], 1}] • Looks pretty good, I'm about testing...
– JHT
Jun 9, 2018 at 17:19
• CellCellAdjacencyMatrix is not much faster getting the neighbors, but the code in history speeds up the thing massively, thanks!
– JHT
Jun 9, 2018 at 17:48
• You're welcome! Jun 9, 2018 at 17:51
• To slow it down Manipulate[ Dynamic[plot[t]], {t, 1, Length[history], 1, AnimationRate -> 2, RefreshRate -> 60}]
– JHT
Jun 9, 2018 at 18:06

I tried a couple of rules and its behavior.

{1,2},{2} means cells with 1 or 2 active neighbors stay active, dead cells with 2 active neighbors get active

{1}, {2} (* divergent *)
{2}, {2} (* divergent, interesting *)
{1, 2}, {2} (* divergent, partial convergent *)
{1, 4}, {2} (* divergent *)
{2, 3}, {2} (* divergent *)
{2, 4}, {2} (* divergent *)
{1, 3}, {2} (* divergent *)
{1, 3, 4}, {2} (* divergent *)
{2, 3}, {2, 3} (* divergent *)
{1, 2}, {2, 3} (* divergent *)
{1, 3}, {2, 3} (* divergent *)
{1, 4}, {2, 3} (* divergent *)
{1, 3, 4}, {2, 3} (* divergent *)
{1, 2, 4}, {2, 3} (* divergent *)

{1}, {3} (* convergent *)
{2}, {3} (* convergent *)
{3}, {3} (* convergent *)
{4}, {3} (* convergent *)
{3, 4}, {3} (* convergent *)
{2, 4}, {3} (* convergent *)
{1, 2}, {3} (* convergent *)
{2, 3}, {3} (* convergent *)
{1, 3}, {3} (* convergent *)
{1, 4}, {3} (* convergent *)
{1, 2, 3}, {3} (* convergent *)
{1, 2, 4}, {3} (* convergent *)
{2, 3, 4}, {3} (* convergent *)
{1, 3, 4}, {3} (* convergent *)
{3}, {2} (* convergent *)
{3, 4}, {2} (* convergent *)
{1, 2, 3}, {2} (* convergent *)
{2, 3, 4}, {2} (* convergent *)
{1, 2, 3, 4}, {2} (* convergent *)
{1, 2, 3}, {2, 3} (* convergent *)
{2, 3, 4}, {2, 3} (* convergent *)
{1, 2, 3, 4}, {2, 3} (* convergent *)


Still found no gliders, but hard to see...

• Maybe you have more luck on a more regular grid: pts = Flatten[ Table[N@{3 i + 3 ((-1)^j + 1)/4, Sqrt/2 j}, {i, 20}, {j, 40}], 1];. Jun 9, 2018 at 18:01
• Good idea, so one can study hexagonal grids and others too.
– JHT
Jun 9, 2018 at 18:15