Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?
share|improve this question
    
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
    
Honestly, I don't have a clue what the problem might be here (i.e. why Aspectrum is not being inlined; everything else looks okay except that there's obviously something amiss with the CompilePrint output given the InputForm). Can you upload the complete code somewhere? –  Oleksandr R. Oct 28 '12 at 3:06
1  
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
show 4 more comments

1 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 *)
share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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