# Catching and debugging numerical errors in compiled functions

I have a compiled function (that was created using the second method of this question). Basically it's using a compiled closure. The problem is that for some parameter values it does not evaluate. Here is a small example:

(*just some function that has a singularity at x=5*)
f = Compile[{{x, _Real, 0}, {a, _Real, 0}}, Sin[a*x*Pi/180]/(x - 5)
, RuntimeOptions -> {"Speed", "EvaluateSymbolically" -> False}];

(*The wrapper just increments the x values by diff. Aa is the
variable that is supposed to be inlined in the closure*)
Block[{Aa},
With[{wrapper =
Compile[{{x, _Real, 0}, {diff, _Real, 0}},
Evaluate@f[x + diff, Aa]]},
min = Compile[{{xs, _Real, 1}, {Aa, _Real, 0}, {diff, _Real, 0}},
First@Sort[wrapper[#, diff] & /@ xs]
, CompilationOptions -> {"InlineCompiledFunctions" -> True,
"InlineExternalDefinitions" -> True,
"ExpressionOptimization" -> True},
RuntimeOptions -> {"Speed"}](*compile*)
]
];

(*The warning is harmeless*)
(*CompiledFunction::cfsa:Argument diff+x at position 1 should be a machine-size real number. >>*)

(*The function is compiled*)
Needs["CompiledFunctionTools"]
CompilePrint[min]

3 arguments
9 Integer registers
11 Real registers
3 Tensor registers
Underflow checking off
Overflow checking off
Integer overflow checking off
RuntimeAttributes -> {}

T(R1)0 = A1
R0 = A2
R1 = A3
I7 = 0
I6 = -5
I3 = 1
R6 = 3.141592653589793
I5 = 180
Result = R7

1   I2 = Length[ T(R1)0]
2   I8 = I7
3   T(R1)1 = Table[ I2]
4   I4 = I7
5   goto 20
6   R3 = GetElement[ T(R1)0, I4]
7   R2 = R1
8   R5 = R2 + R3
9   R4 = R0
10  R8 = I5
11  R9 = Reciprocal[ R8]
12  R8 = R6 * R9
13  R9 = R4 * R5 * R8
14  R8 = Sin[ R9]
15  R9 = I6
16  R10 = R5 + R9
17  R9 = Reciprocal[ R10]
18  R8 = R8 * R9
19  Element[ T(R1)1, I8] = R8
20  if[ ++ I4 < I2] goto 6
21  T(R1)2 = Sort[ T(R1)1]]
22  R7 = Part[ T(R1)2, I3]
23  Return


But then here are the use cases:

(*everything works OK even though the singularity is in the arguments*)
min[{1,5,5},1,0]
(*=> -0.0043631*)


But for some arguments it does not:

(*But if all the arguments are singularites: BLAM!*)
min[{5,5,5},1,0]

(*=> CompiledFunction[{x,a},Sin[(a x \[Pi])/180]/(x-5),-CompiledCode-][5,Aa] *)

   CompiledFunction::cfne: Numerical error encountered; proceeding with uncompiled evaluation. >>
CompiledFunction::cfse: Compiled expression Aa should be a machine-size real number. >>
CompiledFunction::cfex: Could not complete external evaluation at instruction 2; proceeding with uncompiled evaluation. >>


Questions:

• Why is this working if some of the arguments give division by zero, but does not evaluate if all of the arguments are incorrect?
• How can one debug such problems? There is nothing obviously wrong at instruction 2 in the compiled output...
• How to guard against and/or catch such errors in compiled code?
• I can't see the code so it's hard to investigate. But in here, I get a result, with these 2 warnings: "CompiledFunction::cfse: Compiled expression blalbab should be a machine-size real number", and "CompiledFunction::cfex: Could not complete external evaluation at instruction 27; proceeding with uncompiled evaluation" – Rojo Oct 26 '12 at 14:34
• So perhaps something inside your code doesn't return what it should for those inputs? – Rojo Oct 26 '12 at 14:35
• @Rojo Thanks for confirming the error, I thought I was going crazy... :=) I don't see what could possibly be returning an unexpected output (but I'll triple check). By the way: This function is a compiled Nelder-Mead algorithm. – Ajasja Oct 26 '12 at 19:47
• The problem arises when you try to return an infinity; that this occurs only when the arguments are all the same is a side effect of your returning the smallest element (try Last@Sort[...] instead). Interestingly the VM appears to be able to deal with machine infinities during calculations, so it's really a question of how to convert these into ComplexInfinity without causing a numerical error. If you compile this to C you could manually edit the C code, but I don't really know how to do it just using Mathematica. – Oleksandr R. Oct 30 '12 at 13:45
• Er... wherever those parameters are used, could you switch param for Abs[param]? The N-M method doesn't readily admit constraints imposed by returning a large value (which is a method I often use with differential evolution) because the simplex can easily get stuck in a marginally feasible region with this approach. Instead I think you will have better luck if you can transform your objective function so that its domain is bounded. Be careful that it remains unimodal, otherwise you'll run into additional difficulties. – Oleksandr R. Oct 30 '12 at 15:49

Note: the following applies only to Mathematica 8 and above. In previous versions the numerical error occurs immediately when a NaN or inf value is produced, so that avoiding it will be much more difficult.

The reason this occurs is that, although the VM can work with IEEE NaN/inf internally without a problem, Mathematica does not actually support these values at a higher level, using Indeterminate and DirectedInfinity instead. But, these are symbolic values that cannot appear in a packed array, so that it is basically impossible to pass them in and out of compiled code. As a result, even though the compiled code arrives at the expected answer without a problem, the attempt to return this to the top level triggers a numerical error.

To avoid it, one can perhaps test the result of the calculation and send back not only the value itself, but also a flag indicating the existence of a non-numerical or infinite quantity if this appears. For instance:

cf = With[{$MaxMachineNumber =$MaxMachineNumber},
Compile[{{a, _Real, 0}, {b, _Real, 0}},
Block[{
inf = 2. \$MaxMachineNumber,
result = a/b, flag = 0
},
Which[
result == inf, result = 0; flag = 1,
result != result, result = 0; flag = 2,
True, Null
];
{result, flag}
]
]
];

cf[1, 2]
(* -> {0.5, 0.}, result is numeric *)

cf[1, 0]
(* -> {0., 1.}, result is infinity *)

cf[0, 0]
(* -> {0., 2.}, result is not-a-number *)
`