2
$\begingroup$

I am investigating Compile lately and I came across the following problem.

I define a 2-variable function

  f[x_, y_] := 
  If[x < 0 || y < 0, 0, 
  PDF[PoissonDistribution[3], x]*PDF[PoissonDistribution[2], y]
];

and then I have a new function, which I intend to compile. The non-compiled version of it is:

g[gh_, ga_] := 
  Sum[f[i, j], {i, 1 - gh - ga, 10}, {j, 0, i - 1 + gh - ga}];

So when I call for N[g[0, 0] I get 0.584997.

This is my attempt to compile it

gC = Compile[{{gh, _Integer}, {ga, _Integer}}, 
   Sum[f[i, j], {i, 1 - gh - ga, 10}, {j, 0, i - 1 + gh - ga}]
   ];

and when I call for N[gC[0, 0] I get the very same answer (0.584997) along with the following error.

CompiledFunction::cfex: Could not complete external evaluation at instruction 13; proceeding with uncompiled evaluation. >>

Can anyone understand what am I doing wrong?

$\endgroup$
  • 1
    $\begingroup$ Look at mathematica.stackexchange.com/questions/1096/… - you'll find PDF is not compilable. $\endgroup$ – dr.blochwave Mar 7 '16 at 12:15
  • 1
    $\begingroup$ You're against rule 1, 2, 3 here. $\endgroup$ – xzczd Mar 7 '16 at 12:18
  • 1
    $\begingroup$ Try f[x_, y_] = PDF[PoissonDistribution[3], x] PDF[PoissonDistribution[2], y] // FunctionExpand instead and report back. $\endgroup$ – J. M. will be back soon Mar 7 '16 at 12:21
  • $\begingroup$ @blochwavw @xzczd I followed xzczd's methos in my previous post here mathematica.stackexchange.com/questions/108121/… and created the following (more naive) problem, and compile worked fine. f1[x_]:=PDF[PoissonDistribution[l],x] cf1 = Compile[{l}, #,RuntimeOptions -> "EvaluateSymbolically" -> False] &@Sum[f1[i], {i, 0, 2.5}]; FindRoot[cf1@l == 0.4605263157894738, {l, 1}] $\endgroup$ – Tom Zinger Mar 7 '16 at 14:06
  • 1
    $\begingroup$ Guys, I don't think this question is simple or is easily found in the documentation… $\endgroup$ – xzczd Mar 8 '16 at 2:39
7
$\begingroup$

Your function is simple enough that you don't need to rely on pre-defined functions like PDF or PoissonDistribution, just code it all yourself

gC = Compile[{{gh, _Integer}, {ga, _Integer}}, 
   Sum[If[x >= 0 && y >= 0, 3^x/(E^3 x!) 2^y/(E^2 y!), 0], {x, 
     1 - gh - ga, 10}, {y, 0, x - 1 + gh - ga}]];

gC[0, 0]
(* 0.584997 *)

Edit

As pointed out by J.M., you are safer defining the function as

gC = Compile[{{gh, _Integer}, {ga, _Integer}}, 
   Sum[If[x >= 0 && y >= 0, 
     Exp[x Log[3] - 3 - LogGamma[x + 1]] Exp[
       y Log[2] - 2 - LogGamma[y + 1]], 0], {x, 1 - gh - ga, 10}, {y, 
     0, x - 1 + gh - ga}]];

which does not give numerical errors for any of the inputs I tested.

$\endgroup$
  • 4
    $\begingroup$ Factorial[] is not compilable tho. For safety against overflow, 2^y/(E^2 y!) should be replaced by Exp[y Log[2] - 2 - LogGamma[y + 1]], and similarly for the other Poisson PDF. $\endgroup$ – J. M. will be back soon Mar 7 '16 at 12:26
  • $\begingroup$ @J.M. this is interesting to me, because my first reaction to your comment was that it obviously does compile in this instance. But it appears to fail when the number is above 20: If you define fac = Compile[{{y, _Integer}}, y! ]; then fac[20] executes without fail, but fac[21] falls back to uncompiled evaluation. $\endgroup$ – Jason B. Mar 7 '16 at 13:19
  • $\begingroup$ But, if you instead define fac = Compile[{{y, _Real}}, y!]; then the threshold is higher, fac[170] evaluates, but fac[171] falls back to the uncompiled form. $\endgroup$ – Jason B. Mar 7 '16 at 13:21
  • $\begingroup$ @JasonB Yes you are right. My initial goal was to use it in a similar case, where f was recursive. I got the very same errors I get with my initial post's error. Then I thought that it is the already discussed (in this site) issue with compile and recursive functions, but to my surprise problem occurred to a far more simple function too. $\endgroup$ – Tom Zinger Mar 7 '16 at 15:01
1
$\begingroup$

Just wanted to provide an answer to the OP's question

Can anyone understand what am I doing wrong?

which I feel hasn't really been addressed. The original error stems from the compiler assuming f[i,j] returns an integer instead of a real (when i and j are positive), because the inputs of the compiled function gC are integers and the type of f[i,j] is unspecified. To see that this is the case, On["CompilerWarnings"] prior to compilation and check the problematic instruction 13 of

Needs["CompiledFunctionTools`"];
CompilePrint@gC

to see

13 I16 = MainEvaluate[ Hold[f][ I14, I15]]

To my understanding, which may be wrong because the output of CompilePrint isn't documented, CompilePrint prefixes integer registers with I, so this instruction implies an integer register is being assigned the result of f. When N[gC[0,0]] is evaluated, f[i,j] is called with positive i and j, returning a real, which instruction 13 attempts to assign to an integer register. Because this instruction cannot be executed, the kernel reverts to uncompiled evaluation, which works because it makes no such assumptions on the return type of f.

Also, the MainEvaluate in instruction 13 also requires a call to the Wolfram Language evaluator, which means f is effectively uncompiled. This is less relevant to the OP's question because it causes a performance penalty instead of an error, which the OP was concerned about. Also, just because f is uncompiled doesn't mean the compiled code will definitely not be sped up, since Sum is still compiled. However, the accepted answer provides a way to compute the same answer while avoiding MainEvaluate, and should be preferentially used in production runs over other solutions.

For purely academic interest though, we may avoid the error by making sure the compiler knows the return type of f is real. There are several ways to do this, but I will just describe the two most direct ones. The canonical way is using the 3rd argument of compile

gC = Compile[{{gh, _Integer}, {ga, _Integer}}, 
Sum[f[i, j], {i, 1 - gh - ga, 10}, {j, 0, i - 1 + gh - ga}], 
{{_f, _Real}}];

The less memory efficient way is to compile with real input,

gC = Compile[{{gh, _Real}, {ga, _Real}}, 
Sum[f[i, j], {i, 1 - gh - ga, 10}, {j, 0, i - 1 + gh - ga}]];

In this case, the compiler still warns it will assume the return type of f. However, it correctly assumes f returns a real, given gC takes real arguments.

I should stress that these methods do not avoid the MainEvaluate and will not be as efficient as the accepted answer.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.