How to compile this function

Question

I want to compile my function that calculates the energy of a particle system. This is what I've tried:

n = 512;
nl = Ceiling[Power[n, (3)^-1]];
r = 1.;
σ = 2 r;
ρ = 0.01/((2 r)^3);
l = Power[n/ρ, (3)^-1];
ε = 1;
εr = .01;
z = 1;

v = Flatten[Table[l ({i, j, k} - .5)/(nl), {i, nl}, {j, nl}, {k, nl}], 2];
c = RandomSample[Join[Table[-1, {n/2}], Table[1, {n/2}]]];
tc = Times[#[[All, 1]], #[[All, 2]]] &@Subsets[c, {2}];

UC = 
 Compile[{{vv, _Real, 2}}, 
  Total[(4 ε ((σ/#[[1]])^12 - (σ/#[[1]])^6) + 
       1/(4 π εr) (#[[2]] z)/#[[1]]) & /@ 
    Transpose[{Sqrt[Total[((#[[1]] - #[[2]]) - 
             l*Round[(#[[1]] - #[[2]])/l])^2]] & /@ Subsets[vv, {2}],
             tc
    }]
  ]
  , CompilationTarget -> "C", 
  RuntimeOptions -> {"Speed", "CatchMachineUnderflow" -> True}, 
  CompilationOptions -> {"InlineExternalDefinitions" -> True}]

UC[v]

However, the compilation of UC never finishes and after some time, the kernel becomes unresponsive and has to be killed. So I can't get the speed advantages of compilation. How can I persuade Mathematica to compile my function? I don't think the function is too complex to be compiled, there are only basic operations used.

Since you only try to compile some of your code, while referencing other uncompiled parts, this is probably not what you want, i.e. I don't think this will give you any benefits. Try to compile to WML instead of C, and use CompilePrint to look at the compiled code. It calls the main evaluator almost everywhere. — Marius Ladegård Meyer, Jun 02 '15 at 12:16
You can look here to see which built-in functions are compilable. Notice for instance that Subsets are not in the list. — Marius Ladegård Meyer, Jun 02 '15 at 12:19
Thanks, I missed the Subsets function. However, I can still put the Subsets part outside and the rest should compile. — shrx, Jun 02 '15 at 12:22
Evaluate this System`SetSystemOptions["CompileOptions" -> "CompileReportExternal" -> True] and then try to compile your function. — Sektor, Jun 02 '15 at 13:07

Marius Ladegård Meyer · Answer 1 · 2015-06-02T14:19:10.857

Here's how I would compile it. As Oleksandr R. commented, we want to inline everything except the huge lists.

UC = Compile[
  {{v, _Real, 3}, {tc, _Real, 1}},
  Block[{inner},
   inner =
    Sqrt[Total[((#[[1]] - #[[2]]) - 
           l*Round[(#[[1]] - #[[2]])/l])^2]] & /@ v;
   Total[4 ε ((σ/inner)^12 - (σ/
         inner)^6) + (tc z)/((4 π εr) inner)]
   ]
  , CompilationTarget -> "C", 
  RuntimeOptions -> {"Speed", "CatchMachineUnderflow" -> True}, 
  CompilationOptions -> {"InlineExternalDefinitions" -> True, 
    "InlineCompiledFunctions" -> True}
  ]
UC[Subsets[v,{2}],tc]

A crucial improvement that makes both compiled and uncompiled versions much faster is to apply the outer function to the full vectors inner and tc, instead of Maping to each pair of elements.

Note however that since the values of n, nl etc. are "hard coded" at Compile-time, one would need to recompile if one changes any of the parameters that in turn can change v or tc. A way around this is, of course, to simply make all the parameters arguments to the compiled function.

shrx · Accepted Answer · 2015-06-02T14:28:01.220

It seems the problem lies with the inclusion of tc, which in this case contains 130816 elements, in InlineExternalDefinitions. This is apparently too much for Mathematica to handle, but it never complains about it. I'm not aware of any documented limitations of InlineExternalDefinitions, so users should take the appropriate precautions when including large amounts of data in compiled functions.

Update

Thanks to the helpful suggestions from @MariusLadegårdMeyer and @Sektor, I managed to rewrite the function so that it can be fully compiled except for the Subsets part. The compiled function runs about twice as fast as the solution provided by @MariusLadegårdMeyer, however in this case all variables need to be passed to the compiled function at runtime (not an issue):

UC2 = Compile[{{vs, _Real, 3}, {tc, _Real, 1}, {l, _Real}, {ε, _Real},
               {σ, _Real}, {εr, _Real}, {z, _Real}},
  Total[(4 ε ((σ/#[[1]])^12 - (σ/#[[1]])^6) + 
       1/(4 π εr) (#[[2]] z)/#[[1]]) &@{Sqrt[
      Total[Transpose[(((#[[1]] - #[[2]]) - 
              l*Round[(#[[1]] - #[[2]])/l])^2) &@Transpose[vs]]]], 
     tc}], CompilationTarget -> "C", 
  RuntimeOptions -> {"Speed", "CatchMachineUnderflow" -> True}, 
  CompilationOptions -> {"InlineExternalDefinitions" -> True, 
    "InlineCompiledFunctions" -> True}
  ]

Speed comparison:

AbsoluteTiming[UC[Subsets[v, {2}], tc]] (*from @MariusLadegårdMeyer*)

{0.207897, -75.0581}

AbsoluteTiming[UC2[Subsets[v, {2}], tc, l, ε, σ, εr, z]] (*my solution*)

{0.096658, -75.0581}

It is probably the sheer size of the code to be compiled that causes a problem. You could try compiling without optimizations (i.e., common subexpression elimination) and see if that allows the process to finish. Alternatively, pass Subsets[c, {2}] and Subsets[vv, {2}] as parameters. — Oleksandr R., Jun 02 '15 at 13:27
@OleksandrR. Can you explain what you mean with "common subexpression elimination" and how to compile without it? — shrx, Jun 02 '15 at 13:46

How to compile this function

2 Answers2