5
$\begingroup$

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[90];
Points = RandomReal[1, {999, 2}];
VD = VoronoiDiagram[Points];

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

VertexConnection = Table[{}, {i, 1, Length@VD[[1]]}];
For[k = 1, k <= Length[VD[[2]]], 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[[2]]], 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[[1]]], UnboundedVertices]]];
pos = First@
   Position[BoundedVertices, 
    First@Nearest[BoundedVertices, {0.5, 0.5}, 1]];
Start = RandomSample[BC, 
  300];(*Take[Select[VD[[2]],Length@Intersection[pos,#[[2]]]>0&][[All,\
1]],2];*)

Manipulate[
 Refresh[
  If[updating || onestep,
   t++;
   onestep = False;
   ActiveCells = 
    Union[Select[ActiveCells, 
      MemberQ[{2, 3, 4}, 
        Length@Intersection[ActiveCells, NC[[#]]]] &], 
     Select[DeadCells, 
      Length@Intersection[ActiveCells, NC[[#]]] == 3 &]];
   DeadCells = Complement[BC, ActiveCells];
   ];
  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.

Steady state

$\endgroup$
  • 1
    $\begingroup$ Don't forget to add Needs["ComputationalGeometry`"]. $\endgroup$ – QuantumDot Jun 8 '18 at 15:36
  • $\begingroup$ 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. $\endgroup$ – MarcoB Jun 8 '18 at 15:39
  • $\begingroup$ Yes, sorry it requires 'ComputationalGeometry'. $\endgroup$ – fwgb Jun 8 '18 at 15:45
  • $\begingroup$ The problem is getting the neighbors for each cell faster and drawing so many polygons. $\endgroup$ – fwgb Jun 8 '18 at 15:46
  • 1
    $\begingroup$ Nice idea! Sounds as if this post was highly relevant. In particular Chip Hurst's answer. $\endgroup$ – Henrik Schumacher Jun 8 '18 at 15:46
5
$\begingroup$

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:

enter image description here

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}]
A = CellCellAdjacencyMatrix[
   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:

enter image description here

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"]];
   dead = Complement[Range[Length[A]], alive];
   Graphics[
    {EdgeForm[Thin], FaceForm[],
     GraphicsComplex[pts,
      {
       FaceForm[Darker@Green], Polygon[cells[[alive]]],
       FaceForm[Darker@Red], Polygon[cells[[dead]]]
       }
      ]},
    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[[1]], _?(Head[#] == Ray &), 1]];
BC = Select[VD[[2]], Length@Intersection[UnboundedVertices, #[[2]]] == 0 &][[All, 1]];
pts = DeleteCases[VD[[1]], _Ray];
cells = VD[[2, BC, 2]];

A = CellCellAdjacencyMatrix[Length[pts], cells];
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}]

enter image description here

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

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...

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.