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;
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]];
     curr[[species]] *= b[[species]];
     curr[[species]] += 
       GaussianFilter[curr[[species]], disp[[species]]], {2}];
     , {species, 2}];
    tmp = (fitness[[1]]*curr[[1]] + fitness[[2]]*curr[[2]]);
    k = Max[curr];
     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}];
   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} ];

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}];
  imgs[[i]] = getImage[current];
  Do[current = stepC[current, b, d, f], {steps}];
  , {i, ngens}]]
 , {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?

  • $\begingroup$ You're aware that GaussianFilter[] isn't compilable? $\endgroup$ Commented Nov 19, 2012 at 12:41
  • $\begingroup$ Ha sorry didn't know that! What would be the workaround then? Define that function yourself in another compiled function or something? $\endgroup$ Commented Nov 19, 2012 at 12:43
  • 1
    $\begingroup$ @TomWenseleers that would be one approach. If you go down that route, here's a bit of advice: don't even try to implement the FFT in Mathematica's compilable subset. Better to use LibraryLink instead, since LibraryFunctions can be called efficiently from within the VM. $\endgroup$ Commented Nov 19, 2012 at 13:30
  • $\begingroup$ Do you know where the time is being taken ? Could you compute the steps of the evolution of the species first and hold them as data, and then apply a ListAnimate to the precomputed data, or does this have to be a real time animation ? $\endgroup$ Commented Nov 19, 2012 at 14:06
  • $\begingroup$ @image_doctor: no it doesn't have to be a real time animation - computing the steps of the evolution first and displaying it using ListAnimate would be fine too! But seems like it would be quite a lot of hassle then to produce a compiled version of my function, so maybe I'll leave it uncompiled after all... :-( $\endgroup$ Commented Nov 19, 2012 at 18:58

2 Answers 2


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_] := 
   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}];

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.


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.