Using ParallelMap to run multiple Game of Life simulations simultaniously

Question

My Assignment and Code

As an exercise in Parallel computing I've been tasked with implementing Conway's GoL in a way such that at every iteration only a fraction of cells are updated. I'm then to run several simulations starting with the same seed, in order to study effect of non synchronous updating.

I've implemented this sequentially as follows:

Initial variables:

gridSize = 10;
updateFraction = 0.9;
sampleSize = 8;
seed = RandomInteger[{0, 1}, {gridSize, gridSize}];

Normal GoL rules, and updated rules for randomly updating a fraction of the cells:

rulestest[{neighbours_, cell_}] := 
 Module[{}, 
  Piecewise[{{1, neighbours == 3}, {0, 
     neighbours <= 1 || neighbours > 3}}, cell]]

rulestestrandom[{neighbours_, cell_, randomseed_}] := 
 Module[{}, 
  Piecewise[{{cell, randomseed == 0}}, 
   Piecewise[{{1, neighbours == 3}, {0, 
      neighbours <= 1 || neighbours > 3}}, cell]]]

Applying the rules to the grid:

step[matrix_] := 
 Map[{Total@Delete[Flatten[#], 5], #[[2, 2]]} &, 
   Partition[ArrayPad[matrix, 1, "Periodic"], {3, 3}, {1, 1}], {2}]
     /.{p_, q_} -> rulestest[{p, q}]

steprandom[matrix_] := 
 Map[Flatten@{Total@Delete[Flatten[#], 5], #[[2, 2]], 
      RandomChoice[{updateFraction, 1 - updateFraction} -> {1, 0},1]}&, 
       Partition[ArrayPad[matrix, 1, "Periodic"], {3, 3}, {1, 1}], {2}]
     /. {p_, q_, r_} -> rulestestrandom[{p, q, r}]

Doing a bunch of iterations and plotting it:

controlMatrix = seed;
controlAnimateSet = {controlMatrix}; 
Do[controlMatrixtemp = step[controlMatrix]; 
 AppendTo[controlAnimateSet, controlMatrixtemp]; 
 controlMatrix = controlMatrixtemp;, {100}]

Map[Function[testMatrix[#] = seed], Range[sampleSize]];
Map[Function[testMatrixtemp[#] = {}], Range[sampleSize]];
Map[Function[testMatrixAnimateSet[#] = {testMatrix[#]}], 
  Range[sampleSize]];
Do[Map[(testMatrixtemp[#] = steprandom[testMatrix[#]]; 
     AppendTo[testMatrixAnimateSet[#], testMatrixtemp[#]]; 
     testMatrix[#] = testMatrixtemp[#];) &, Range[sampleSize]];, {100}]

Manipulate[{Table[
   ArrayPlot[testMatrixAnimateSet[m][[n]], Mesh -> True, 
    ImageSize -> Medium], {m, 1, 8}], 
  ArrayPlot[controlAnimateSet[[n]], Mesh -> True, 
   ImageSize -> Medium]}, {n, 1, 100, 1}]

The problem

The code runs fine just like this, sequentially. However the goal is to make it run in parallel. I figured there are two things that could potentially be parallelized.

Firstly updating the grid could be done in parallel:

step[matrix_] := 
 ParallelMap[{Total@Delete[Flatten[#], 5], #[[2, 2]]} &, 
   Partition[ArrayPad[matrix, 1, "Periodic"], {3, 3}, {1, 1}], {2}]
     /.{p_, q_} -> rulestest[{p, q}]

Now this works fine, but for grid sizes up to 100x100 (Which I've been instructed to use) this is less efficient when it comes to computing time, resulting in speedups of around 0.8 on my quadcore machine.

The second option, which in my mind seems more obvious, is to parallelize the iterations on the group of test matrices. I tried to do it as follows:

Do[ParallelMap[(testMatrixtemp[#] = steprandom[testMatrix[#]]; 
     AppendTo[testMatrixAnimateSet[#], testMatrixtemp[#]]; 
     testMatrix[#] = testMatrixtemp[#];) &, Range[sampleSize]];, {100}]

However this doesn't work, and I don't understand why. The code runs without errors and gives a very nice speedup of around 3.0, however the testMatrixAnimateSet never gets appended with future iterations of the grid. Is this an error in my code, or just Mathematica being quirky? Is there a way to work around this?

Answer 1

My final solution

After working on the code a little bit more, I've come up with a much better way to tackle the problem that I will leave here for posterity

First I define the size of the simulation, create the seed and launch kernels if they aren't already running.

(*Define some variables, generate the seed and launch the kernels.*)
gridSize = 50;
sampleSize = 100;
fractionAmount = 10;
seed = RandomInteger[{0, 1}, {gridSize, gridSize}];
If[Kernels[] == {}, LaunchKernels[]];

Then I define the rules for the GoL, and I define a function to update a grid according to these rules.

(*Define the rules for game of life asynchronously*)
rulestestrandom[{neighbours_, cell_, randomseed_}] := 
 Module[{}, 
  Piecewise[{{cell, randomseed == 0}}, 
   Piecewise[{{1, neighbours == 3}, {0, 
      neighbours <= 1 || neighbours > 3}}, cell]]]

(*Define an iteration of the grid using the rule set*)
steprandom[matrix_, updateFraction_] := 
 Map[Flatten@{Total@Delete[Flatten[#], 5], #[[2, 2]], 
      RandomChoice[{updateFraction, 1 - updateFraction} -> {1, 0}, 
       1]} &, Partition[
    ArrayPad[matrix, 1, "Periodic"], {3, 3}, {1, 1}], {2}] /. {p_, q_,
     r_} -> rulestestrandom[{p, q, r}]

I then define a function that allows me to update a single grid over a certain amount of time steps, while updating only a specific fraction of the cells in the grid.

(*Define the advancing of a single grid, updating a certain fraction \
of cells, over gridSize/2 timesteps*)
evaluateGrid[n_] := Module[
  {testMatrix = seed,
   testMatrixtemp,
   testHammingSet = {},
   testDensitySet = {},
   updateFraction = n/fractionAmount},
  Do[
   testMatrixtemp = steprandom[testMatrix, updateFraction];
   AppendTo[testDensitySet, Total[Flatten@testMatrix]/gridSize^2 // N];
   AppendTo[testHammingSet, 
    HammingDistance[Flatten@testMatrixtemp, Flatten@testMatrix]];
   testMatrix = testMatrixtemp;
   , {gridSize/2}];
  {testDensitySet, testHammingSet}]

Then, in parallel, I apply this function to the seed a certain amount of times to generate a nice sample of grids. I do this for different fractions. As output I get a lists of densities and hamming distances that allows the behaviour of the system to be studied.

solutions = 
  Table[ParallelTable[evaluateGrid[n], {sampleSize}], {n, 1, 
    fractionAmount, 1}];

To visualize the data I use the following code:

dEnvelope[n_] := {Max@#, Min@#} & /@ (Transpose@solutions[[n, All, 1]])
hEnvelope[n_] := {Max@#, Min@#} & /@ (Transpose@solutions[[n, All, 2]])
Manipulate[
 ListPlot[{dEnvelope[n][[All, 1]], dEnvelope[n][[All, 2]]}, 
  AxesLabel -> {"Time", "Density"}, PlotLabel -> "Density over time", 
  Joined -> True, Filling -> True, 
  PlotRange -> {0.1, 0.6}], {{n, 1, "Fraction of cells updated"}, 1, 
  fractionAmount, 1}]
Manipulate[
 ListPlot[{hEnvelope[n][[All, 1]], hEnvelope[n][[All, 2]]}, 
  AxesLabel -> {"Time", "Hamming distance"}, 
  PlotLabel -> "Hamming distance over time", Joined -> True, 
  Filling -> True, 
  PlotRange -> {0, gridSize^2/2}], {{n, 1, 
   "Fraction of cells updated"}, 1, fractionAmount, 1}]
Manipulate[
 BoxWhiskerChart[Transpose@solutions[[n, All, 1]], "Outliers", 
  PlotLabel -> "Distribution of Density over time", 
  PlotRange -> {0.1, 0.6}], {{n, 1, "Fraction of cells updated"}, 1, 
  fractionAmount, 1}]
Manipulate[
 BoxWhiskerChart[Transpose@solutions[[n, All, 2]], "Outliers", 
  PlotLabel -> "Distribution of Hamming distance over time", 
  PlotRange -> {0, gridSize^2/2}], {{n, 1, 
   "Fraction of cells updated"}, 1, fractionAmount, 1}]

Unlike the previous method, which achieved negative speedup because the worker kernels were constantly waiting on eachother to acces the same data, this method achieves a high speedup.