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

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?

-
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. –  rm -rf Nov 16 '12 at 23:46
thx for letting me know - I'm new here as you probably gathered :-) –  Tom Wenseleers Nov 16 '12 at 23:49
No problem; I hope you have fun on this site :) Just for reference, here's a list of editing tips. –  rm -rf Nov 17 '12 at 0:23

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}]
]

-
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? –  Tom Wenseleers Nov 17 '12 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}];
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 = InternalBag[]},
Do[InternalStuffBag[population,
space = step[b, mu, stepsize, y, r, e, space], 2], {nrsteps}];
Partition[#, gridsize] & /@
Partition[InternalBagPart[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]
`

-
Thanks so much - wow that was a lot faster indeed! –  Tom Wenseleers Nov 17 '12 at 21:41