Skip to main content
1 of 2
Michael E2
  • 258.7k
  • 21
  • 370
  • 830
  1. The problem arises in the compiled function created for some sort of singularity processing that NIntegrate[] decided was necessary.

  2. Is it a bug? Ask WRI. Here are two things to consider. (a) The singularity processing is not necessary. (b) The compiled function created has a return type Complex even though the function is Real. The function being compiled is derived from the integrand and is slightly more complicated (a linear combination of Re[] expressions). Compile[] seems to gives up on the analysis of this function and defaults to the safer type Complex for the result.

  3. The comments show some ways to avoid a complex output. A direct fix is to turn off the singularity handler:

NIntegrate[ Re[fg[{kx, ky}] Exp[2 \[Pi] I (kx + ky)]], {kx, -0.5, 0.5}, {ky, -0.5, 0.5}, MinRecursion -> 1, (* stops the GlobalAdaptive slow-con warning *) Method -> {"GlobalAdaptive", "SingularityHandler" -> None}] (* 0.014725 *) 

Here's how to analyze the integration. IntegrationMonitor returns information about the integral over each subregion.

NIntegrate[ Re[fg[{kx, ky}] Exp[2 \[Pi] I (kx + ky)]], {kx, -0.5, 0.5}, {ky, -0.5, 0.5}, IntegrationMonitor :> ((foo = #) &)] 

Pick out the subregions over which the integral is complex:

Position[Through[foo["Integral"]], _Complex] cplx = Extract[foo, %]; (* {{137}, {138}, {139}, {140}, {141}, {142}, {146}, {147}, {148}, \ {149}, {150}, {155}, {156}, {327}, {328}, {333}, {336}, {337}, {338}, \ {339}, {342}, {344}, {346}, {347}, {348}, {349}} *) 

Each subregion has an Experimental`NumericalFunction[] used to evaluate the integral, and the numerical function has a CompiledFunction (in the OP's case). Unsurprisingly, the ones related to Real integral values have Real integrand values, and those with Complex integral values have Complex integrand values.

nfR = foo[[1]]["NumericalFunction"]; cfR = nfR@"CompiledFunction"; nfC = cplx[[1]]["NumericalFunction"]; cfC = nfC@"CompiledFunction"; nfR[-0.375`, -0.375`] cfR[-0.375`, -0.375`] (* -0.125894 -0.125894 *) nfC[0.96875`, 1.] cfC[0.96875`, 1.] (* 0.00102451 + 0. I 0.00102451 + 0. I *) 

We can check that the return type agrees with the results (a redundant step).

Needs@"CompiledFunctionTools`" Cases[ToCompiledProcedure[cfR], _CompiledResult, Infinity] (* {CompiledResult[Register[Real, 2]]} *) Cases[ToCompiledProcedure[cfC], _CompiledResult, Infinity] (* {CompiledResult[Register[Complex, 17]]} *) 

One can examine the end of the output of CompilePrint@cfR and CompilePrint@cfC to see that both apply Re[] etc. The output is rather long, and I will omit it. One can also see whether the result is a real or complex register.

We can test the hypothesis that it is Compile that fails to parse the Real value of the function:

testC = Compile[{{kx, _Complex}, {ky, _Complex}}, Evaluate@nfC@"FunctionExpression"] testC[0.96875`, 1.] (* 0.00102451 + 0. I *) 

Here's the form of the function in nfC[]:

Shallow[nfC@"FunctionExpression", 5] (* kx (0.00236289 Re[<<1>>] + 0.00236289 Re[<<1>>]) *) 

Clearly, one wants the time Compile[] spends analyzing code to be very short, especially in auto-compiled code. It's hard to argue strongly that it is a bug. Normally, if the user is calling Compile[] there are workarounds. I don't know of any hooks in NIntegrate[] to control Compile[].

Michael E2
  • 258.7k
  • 21
  • 370
  • 830