Compile function - error but correct calculation

Question

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?

Answer 1

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.

Answer 2

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.