Error in Compiled function at runtime in spatial competition model

Question

I just did a model of spatial competition between two species with different competitive ability & dispersal rates on a lattice with variables nrs of individuals per site. To speed it up I tried to use Compile - but I get an error at runtime. It works OK though if I run it uncompiled. Anybody any ideas what I might be doing wrong? My code is

gridsize = 30;
k = 1000;
zeroMatrix[n_, m_] := Table[0, {n}, {m}];
initSpace[n_, m_, k_] := Block[{s},
   s = { zeroMatrix[n, m], zeroMatrix[n, m] };
   s[[1, 1, 1]] = k;
   s[[2, n, m]] = k;
   Return[s];
   ];
current = initSpace[gridsize, gridsize, k];

This function does one time step iteration, using synchronous updating:

stepC = Compile[{{current, _Real, 3}, {b, _Real, 1}, {disp, _Real, 
     1}, {fitness, _Real, 1}},
    Block[{tmp, k, k2, n, m, curr = current},
    n = Dimensions[curr][[2]];
    m = Dimensions[curr][[3]];
    Do[
     curr[[species]] *= b[[species]];
     curr[[species]] += 
      Map[Round, 
       GaussianFilter[curr[[species]], disp[[species]]], {2}];
     , {species, 2}];
    tmp = (fitness[[1]]*curr[[1]] + fitness[[2]]*curr[[2]]);
    k = Max[curr];
    Do[
     k2 = If[tmp[[i, j]] < k, tmp[[i, j]], k];
     tmp[[i, j]] = If[ tmp[[i, j]] == 0, 1, tmp[[i, j]] ];
     curr[[1, i, j]] = k2*fitness[[1]]*curr[[1, i, j]] / tmp [[i, j]];
     curr[[2, i, j]] = k2*fitness[[2]]*curr[[2, i, j]] / tmp [[i, j]];
     , {j, m}, {i, n}];
    curr
    ],
   CompilationTarget -> "C", Parallelization -> True, 
   RuntimeOptions -> "Speed"];

I then define a function to display the densities of the two species:

getImage[current_] := Module[{img},
   img = Map[{#/k, 0, 0} &, current[[1]], {2} ];
   img += Map[{0, 0, #/k} &, current[[2]], {2} ];
   Return[Graphics[Raster[img]]];
   ];

And iterate and display an animation of the dynamics:

b = {2.0, 2.0};
d = {3.0, 1.0};
f = {1.0, 1.3};
current = initSpace[gridsize, gridsize, k];
ngens = 30;
steps = 5;
imgs = Table[0, {ngens}];
Timing[Do[
  imgs[[i]] = getImage[current];
  Do[current = stepC[current, b, d, f], {steps}];
  , {i, ngens}]]
Animate[
 imgs[[i]]
 , {i, 1, ngens, 1}, AnimationRepetitions -> 1, 
 AnimationRunning -> False]

During runtime I then get an error "CompiledFunction::cfne: Numerical error encountered; proceeding with uncompiled evaluation." Any thoughts what's wrong?

Answer 1

You can get some speed increase by writing your functions differently. It is almost always faster to act on whole lists in Mathematica than to explicitly loop over the individual elements.

A couple of notes:

• Round is listable and does not have to be mapped over the array.
• The definition of tmp is just a dot product.
• The alterations to tmp in your If statements can be done using Clip and Unitize, which are both listable.
• The update to curr can also be done in a listable fashion, removing the need to loop over i and j.

Here are alternative definitions for stepC and getImage which provide something like 5 to 10 times speed increase:

step2[current_, b_, d_, f_] := Module[{curr, tmp, k2},
  curr = MapThread[Times, {current, b}];
  curr += Round@MapThread[GaussianFilter, {curr, d}];
  tmp = f.curr;
  k2 = Clip[tmp, {-Infinity, Max[curr]}];
  tmp = tmp + 1 - Unitize[tmp];
  MapThread[#1 #2 k2/tmp &, {curr, f}]]

getImage2[current_] := 
 Graphics@Raster[
   Apply[{#1, 0, #2} &, Transpose[current, {3, 1, 2}], {2}]/k]

I'd also use Nest instead of a Do loop for the individual iterations, and Table to collect the results into imgs:

imgs = Table[
   getImage2[current = Nest[step2[#, b, d, f] &, current, steps]], {i, ngens}];

Answer 2

As already pointed out in the comments, compiling to "C" is useless, when you use e.g. functions like GaussianFilter wich cannot be compiled down and which require a call to the evaluating kernel.

Regarding your error message: Why don't you investigate in the parameters to stepC when the error occurs? You could do this by wrapping a Check around which stops when a message is thrown and the you insert debugging code at this place. You could for instance try to define a new step function

step[current_, b_, d_, f_] := Check[stepC[current, b, d, f],
   Throw[{current, b, d, f}]];

which you then call in your Do loop

res = Catch[Timing[Do[imgs[[i]] = getImage[current];
    Do[current = step[current, b, d, f], {steps}];, {i, ngens}]]]

Now res contains the erroneous input.