Catching and debugging numerical errors in compiled functions

Question

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?

Answer 1

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 *)