Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

share|improve this question
    
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
add comment

2 Answers

up vote 8 down vote accepted

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}]
  ]
share|improve this answer
    
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
add comment

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]

enter image description here

share|improve this answer
    
Thanks so much - wow that was a lot faster indeed! –  Tom Wenseleers Nov 17 '12 at 21:41
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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