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}]
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:
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}]
Needs["ComputationalGeometry`"]. $\endgroup$