A problem with Compile and Map

Question

There is a weird behavior of Compile that I do not understand. Consider the following example:

(*initial parameters*)
lL = 64.;
pi = 3.14;

tab = Flatten[Table[{kx, ky},
{kx, 0., pi, pi/lL}, {ky, 0., pi, pi/lL}], 1];

enn = Compile[{{kx, _Real}, {ky, _Real}}, Cos[kx] Cos[ky], 
CompilationTarget -> "C"];

First approach

fF = Compile[{{En, _Real}}, 1./(1. + Exp[500. En]), 
   CompilationTarget -> "C"];

sumT = Compile[{{tab, _Real, 2}},
   Total[Map[fF[enn[#[[1]], #[[2]]]] &, tab, 1]],
   CompilationTarget -> "C"];

sumT[tab] // AbsoluteTiming
{2.389137, 131584.5}

Second approach

ffd = Compile[{{kx, _Real}, {ky, _Real}}, 
   1./(1. + Exp[500. enn[kx, ky]]), {{enn, _Real}}, 
    CompilationTarget -> "C"];

sumt = Compile[{{tab, _Real, 2}},
   Total[Map[ffd[#[[1]], #[[2]]] &, tab]],
   CompilationTarget -> "C"];

sumt[tab] // AbsoluteTiming
{0.061003, 131584.5}

The first approach is much slower and returns the following error:

"Argument {{0.,0.},(...),<<263119>>} 
at position 1 should be a "machine-size real number"

If the argument of fF is declared as _Real, shouldn't it work as well? I've tried to add {{enn,_Real}} as a declaration of subexpression in sumT, but it didn't change a thing.

Somehow, in the problem I try to tackle, I'd prefer to use the first approach (if it would work:-).

Answer 1

While I'm not sure why the error you see is generated, you can fix your sumT function by taking the fF call out of the Map form:

sumT = Compile[
  {{tab, _Real, 2}},
  Total[fF /@ Map[enn[#[[1]], #[[2]]] &, tab]],
  CompilationTarget -> "C"];

This worked fine for me in version 10.1 of Mathematica:

sumT[tab] // AbsoluteTiming
(* ==> {0.000667, 2112.49} *)

This said, you can get nearly identical performance (at least in my own tests) without using compilation at all, assuming you rewrite your code to use listable functions instead of Map:

sumT2[tab_] := With[
  {tr = Transpose@tab},
  Total[1.0 / (1.0 + Exp[500.0 * Cos[tr[[1]]]*Cos[tr[[2]]]])]];

Mean[Table[AbsoluteTiming[#[tab]][[1]], {1000}]]& /@ {sumT, sumT2}
(* ==> {0.000470342, 0.000475503} *)

Answer 2

It happens because Compile cannot infer the types returned by the subsidiary functions enn and fF. In the second approach, you partially solved this by specifying it manually for enn. In principle you could have done that in the first approach as well if you had specified it for fF too, but in practice getting the type inference to work out correctly is not always easy, and this seems to be one of those cases. However, it is not a good option anyway, because it results in a call out of compiled code, which comes with a significant performance penalty.

Actually, both of your approaches are slightly wrong in this respect. You can correct both problems by giving the option CompilationOptions -> "InlineExternalDefinitions" -> True. This way, Compile is able to determine the return types from the code.

I use your given initial parameters, except that lL should be 512, not 64, in order to give the same results as you show in the question.

Approach 1:

fF = Compile[{{En, _Real}}, 1./(1. + Exp[500. En])];

sumT = Compile[{{tab, _Real, 2}}, 
  Total[Map[fF[enn[#[[1]], #[[2]]]] &, tab, 1]], 
  CompilationOptions -> "InlineExternalDefinitions" -> True
 ];

sumT[tab] // Timing
(* -> {0.093750, 131584.} (with no messages) *)

Approach 2:

ffd = Compile[{{kx, _Real}, {ky, _Real}}, 
  1./(1. + Exp[500. enn[kx, ky]]), 
  CompilationOptions -> "InlineExternalDefinitions" -> True
 ];

sumt = Compile[{{tab, _Real, 2}}, 
  Total[Map[ffd[#[[1]], #[[2]]] &, tab]], 
  CompilationOptions -> {
   "InlineExternalDefinitions" -> True, 
   "InlineCompiledFunctions" -> True
  }
 ];

sumt[tab] // AbsoluteTiming
(* {0.093750, 131584.} (with no messages) *)

The code generated for each of these is basically equivalent, so they perform the same.

I should just add that, personally, I don't particularly like using CompilationOptions -> "InlineExternalDefinitions" -> True. I find it better to use With, because it works in older versions as well and gives you greater control over what to inline.

Finally, you're better off not to use compilation here anyway. The excellent suggestion by user21382 to use listable functions results in simpler, clearer, less error-prone, and better than twice as fast code (as tested on my computer with lL set to 512).