Compiling a spatial predator-prey model and optimising it for speed

Question

I just made a spatial predator-prey model but it runs a little slow and I am seeking some advice on how to compile it to optimise it for speed.

I first define matrix with 0 for an empty site, 1 for a predator and 2 for a prey present:

e = 0;r = 1;y = 2;
initSpace[n_] := RandomInteger[{0, 2} , {n, n}];
space = initSpace[10];

I then defined a getCell function which looks for a neighbour at cell position $(i,j)$ in a matrix sp with a certain value x, and returns the coordinates of that cell; if it doesn't find a cell with that value it returns coordinates $(i,j)$.

getCell = 
 Compile[{{sp, _Integer, 
    2}, {i, _Integer}, {j, _Integer}, {x, _Integer}}, 
  Block[{n, m, k2, l2, cells}, {n, m} = Dimensions[sp];
   cells = {{i, j}};
   Do[k2 = Mod[i + k, n, 1]; (* This is the neighborhood *)
    l2 = Mod[j + l, m, 1];
    If[(k2 != i || l2 != j) && sp[[k2, l2]] == x, 
     AppendTo[cells, {k2, l2}]], {l, -1, 1}, {k, -1, 1}];
   If[Length[cells] == 1, {i, j}, RandomChoice[Rest[cells]]]], 
  CompilationTarget -> "C", Parallelization -> True, 
  RuntimeOptions -> "Speed"]

I then define a module step which calculates the space matrix in the next time step based on a given birth rate b and death rate mu after one event at a random location in the space.

step[b_, mu_] := Block[{i, j, k, l, dim},
  dim = Dimensions[space];
  i = Random[Integer, {1, dim[[1]]}];
  j = Random[Integer, {1, dim[[2]]}];
  Which[
   (* Predator *)
   space[[i, j]] == r,
   If[Random[Real] < mu,
    (* Predator dying *)
    space[[i, j]] = e,
    (* Predator eating *)
    {k, l} = getCell[space, i, j, y];
    space[[k, l]] = r;
    ],
   (* Prey *)
   space[[i, j]] == y,
   If[Random[Real] < b,
    (* Prey reproducing *)
    {k, l} = getCell[space, i, j, e];
    space[[k, l]] = y
    ]
   ]
  ]

Finally, I do many iterations

 gridsize = 200;
space = initSpace[gridsize];
nrsteps = 100;
stepsize = 5000;
b = 0.9; mu = 0.4;
population = Table[0, {nrsteps}];
Do[Do[step[b, mu];, {stepsize}];
  population[[i]] = space;
  , {i, 1, nrsteps}];
    Animate[
     ArrayPlot[population[[i]], ColorRules -> {e -> Black, r -> Red, y -> Blue}],
     {i, 1, nrsteps/stepsize, 1}, AnimationRepetitions -> 1, 
     AnimationRunning -> False]

Ideally I would like to have the latter two modules compiled as well, but when I tried to put Compile[] around them I got a bunch of errors. Does anyone know what I might be doing wrong? Or have any other advice on how I should optimise my code for speed? Or should I try to put everything in a single compiled Module?

Answer 1

Because random number generator requires initialization, multiple calls to RandomInteger result in much slower execution than a single call to generate many samples:

In[18]:= Table[RandomInteger[{1, 10}], {10^7}]; // AbsoluteTiming

Out[18]= {1.300000, Null}

In[19]:= RandomInteger[{1, 10}, 10^7]; // AbsoluteTiming

Out[19]= {0.210000, Null}

This is what happens in your code, since each of the numerous calls to function step generate 3 calls to random number generator. You code can be sped-up by a factor of 2 by replacing Do[ step[b,mu], {stepsize}] with a call to steps[b, mu, stepsize], where steps is as follows:

steps[b_, mu_, rep_] := Block[{i, j, k, l, dim, q, is, js, qs, val},
  dim = Dimensions[space];
  is = RandomInteger[{1, dim[[1]]}, rep];
  js = RandomInteger[{1, dim[[2]]}, rep];
  qs = RandomReal[1, rep];
  Do[
   i = Part[is, p]; j = Part[js, p]; q = Part[qs, p];
   val = space[[i, j]];
   Which[
    (*Predator*)val == r,
    If[q < mu,
     (*Predator dying*)space[[i, j]] = e,
     (*Predator eating*){k, l} = getCell[space, i, j, y];
     space[[k, l]] = r;
     ],
    (*Prey*)val == y,
    If[q < b,(*Prey reproducing*)
     {k, l} = getCell[space, i, j, e];
     space[[k, l]] = y]], {p, 1, rep}]
  ]

Answer 2

I think a big part of your problem are the many global variables. Avoiding this and using @Sasha's suggestions, you can compile the step function. Note the use of CompilationOptions that inlines your getCell function.

step = Compile[{{b, _Real}, {mu, _Real}, {rep, _Integer}, {y,  _Integer}, 
   {r, _Integer}, {e, _Integer}, {space, _Integer, 2}}, 
   Block[{i, j, k, l, dim, q, is, js, qs, val, spaceCopy = space}, 
    dim = Dimensions[space];
    is = RandomInteger[{1, dim[[1]]}, rep];
    js = RandomInteger[{1, dim[[2]]}, rep];
    qs = RandomReal[1, rep];
    Do[i = Part[is, p]; j = Part[js, p]; q = Part[qs, p];
     val = spaceCopy[[i, j]];
     Which[(*Predator*)val == r, 
      If[q < mu,(*Predator dying*)
       spaceCopy[[i, j]] = e,(*Predator eating*){k, l} = 
        getCell[spaceCopy, i, j, y];
       spaceCopy[[k, l]] = r;],(*Prey*)val == y, 
      If[q < b,(*Prey reproducing*){k, l} = 
        getCell[spaceCopy, i, j, e];
       spaceCopy[[k, l]] = y]], {p, 1, rep}];
    spaceCopy], 
   CompilationOptions -> "InlineExternalDefinitions" -> True];

You can then re-write your simulation as follows.

sim[y_: 2, r_: 1, e_: 0, gridsize_: 200, nrsteps_: 100, 
  stepsize_: 5000, b_: 0.9, mu_: 0.4] :=
 Block[{space = initSpace[gridsize], population = Internal`Bag[]},
  Do[Internal`StuffBag[population, 
    space = step[b, mu, stepsize, y, r, e, space], 2], {nrsteps}];
  Partition[#, gridsize] & /@ 
   Partition[Internal`BagPart[population, All], gridsize*gridsize]
  ]

This is now fairly fast (though I'm sure additional improvements could be made).

AbsoluteTiming[population = sim[];]

(*{0.6708012, Null}*)

Animate[ArrayPlot[population[[i]], 
  ColorRules -> {e -> Black, r -> Red, y -> Blue}], {i, 1, 
  stepsize/nrsteps, 1}, AnimationRepetitions -> 2, 
 AnimationRunning -> False]