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]]}];
   (* 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}];
     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?

  • $\begingroup$ If you indent your code blocks by 4 spaces (or select it and press the {} button inside the text box on this site), then it gets formatted as code with pretty syntax highlighting. $\endgroup$
    – rm -rf
    Commented Nov 16, 2012 at 23:46
  • $\begingroup$ thx for letting me know - I'm new here as you probably gathered :-) $\endgroup$ Commented Nov 16, 2012 at 23:49
  • $\begingroup$ No problem; I hope you have fun on this site :) Just for reference, here's a list of editing tips. $\endgroup$
    – rm -rf
    Commented Nov 17, 2012 at 0:23

2 Answers 2


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];
   i = Part[is, p]; j = Part[js, p]; q = Part[qs, p];
   val = space[[i, j]];
    (*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}]
  • $\begingroup$ Many thx for this - that goes quite a bit faster! Is there also a possibility to Compile the steps function though? When I tried it strangely enough ends up going more slowly than the uncompiled version. Any ideas why that migtht be the case? $\endgroup$ Commented Nov 17, 2012 at 8:18

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}];
   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[]},
    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}*)

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

enter image description here

  • $\begingroup$ Thanks so much - wow that was a lot faster indeed! $\endgroup$ Commented Nov 17, 2012 at 21:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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