How to use `integration_qawc` with function wrapped with local variables?

General Usage
2

Hi, to know how a @gsl_function will work, you also need to know how a @cfunction will work, and the issue you report is the same as:

You may also find this manual section useful.
One solution for passing a local is to construct the gsl_function callback yourself like this:

using GSL

fdist(ε::Float64, μ::Float64, kbt::Float64) = 1.0 / (exp((ε - μ)/kbt) + 1.0)

function integrand(x::Cdouble, p::Tuple{Cdouble,Cdouble})::Cdouble
    1.0-fdist(x, p...)
end

function genΛ2h(i, j, α, ϵ, t, μ, kbt; ws_size = 100)
    ws = integration_workspace_alloc(ws_size)

    result = Cdouble[0]
    abserr = Cdouble[0]

    p = Ref((α == 1 ? μ[1] : μ[2], kbt))
    GC.@preserve p begin
        ptr = pointer_from_objref(p)
        f = gsl_function(@cfunction(integrand, Cdouble, (Cdouble, Ref{Tuple{Cdouble,Cdouble}})), ptr)
        integration_qawc(f, -1000, 1000, ϵ[j], 1e-9, 1e-9, ws_size, ws, result, abserr)
    end
    integration_workspace_free(ws)
    return - t[i,α] * t[j,α]' * result[1] / (2*π)
end

genΛ2h(1, 1, 1, [1., 1.], [[1e-2 1e-2];[1e-2 1e-2]], [1.5, -1.5], 0.6)

Another option is that the development branch of GSL.jl defines a cconvert, unsafe_convert pair for gsl_function that makes it possible to pass an anonymous function with a closure as the integrand directly to the routine, which is far easier to use.

It may also be more convenient to use QuadGK.jl for this calculation, whose manual gives a method for principal value integrals. If you also notice that 1-fdist causes catastrophic cancellation I would rewrite the problem like this:

using QuadGK

function cauchy_quadgk(g, a, b; kws...)
    a < 0 < b || throw(ArgumentError("domain must include 0"))
    g₀ = g(0)
    g₀int = b == -a ? zero(g₀) : g₀ * log(abs(b/a)) / (b - a)
    return quadgk_count(x -> (g(x)-g₀)/x + g₀int, a, 0, b; kws...)
end
function genΛ2h_gk(i, j, α, ϵ, t, μ, kbt)
    p = α == 1 ? μ[1] : μ[2]
    c = ϵ[j]
    result, err, numevals = cauchy_quadgk(-1000-c, 1000-c; atol=1e-9, rtol=1e-9) do x
        e = exp((x - p + c)/kbt)
        (e > one(e) ? one(e)-one(e)/(e+one(e)) : e/(e+one(e)))
    end
    return - t[i,α] * t[j,α]' * result[1] / (2*π)
end
genΛ2h_gk(1, 1, 1, [1., 1.], [[1e-2 1e-2];[1e-2 1e-2]], [1.5, -1.5], 0.6)
1 Like