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 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?

share|improve this question
    
You're aware that GaussianFilter[] isn't compilable? –  J. M. Nov 19 '12 at 12:41
    
Ha sorry didn't know that! What would be the workaround then? Define that function yourself in another compiled function or something? –  Tom Wenseleers Nov 19 '12 at 12:43
1  
@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. –  Oleksandr R. Nov 19 '12 at 13:30
    
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 ? –  image_doctor Nov 19 '12 at 14:06
    
@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... :-( –  Tom Wenseleers Nov 19 '12 at 18:58
add comment

2 Answers 2

up vote 2 down vote accepted

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

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.

share|improve this answer
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.