# How to use \`integration\_qawc\` with function wrapped with local variables?

**URL:** <https://discourse.julialang.org/t/how-to-use-integration-qawc-with-function-wrapped-with-local-variables/113269>\
**Category:** General Usage\
**Created:** [April 20, 2024, 11:14am UTC](https://discourse.julialang.org/t/how-to-use-integration-qawc-with-function-wrapped-with-local-variables/113269 "2024-04-20T11:14:59Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![zxct](https://avatars.discourse-cdn.com/v4/letter/z/5daacb/32.png) [@zxct](https://discourse.julialang.org/u/zxct)\
**Post date:** [April 20, 2024, 11:14am UTC](https://discourse.julialang.org/t/how-to-use-integration-qawc-with-function-wrapped-with-local-variables/113269/1 "2024-04-20T11:14:59Z")

</div>

Hi, I want to calculate the cauchy principle value integration of a three arguments function inside a function with two arguments decided at runtime with local variables. But it seemed that `integration_qawc` will not find the local variables. Minimal code to repeate the situation is given as follows,

```julia
using GSL

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

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

    if α == 1
        f = @gsl_function(x -> (1 - fdist(x, μ[1], kbt)))
    else
        f = @gsl_function(x -> (1 - fdist(x, μ[2], kbt)))
    end

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

    integration_qawc(f, -1000, 1000, ϵ[j], 1e-9, 1e-9, 100, ws, result, abserr)

    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)

```

with `ERROR: UndefVarError: μ not defined`

Any suggestions for making the code to work?

---

<div class="post-metadata">

**Author:** ![lxvm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lxvm/32/50010_2.png) [@lxvm](https://discourse.julialang.org/u/lxvm)\
**Post date:** [April 21, 2024, 5:59am UTC](https://discourse.julialang.org/t/how-to-use-integration-qawc-with-function-wrapped-with-local-variables/113269/2 "2024-04-21T05:59:47Z")

</div>

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:

> [@Passing a local to @cfunction fails](https://discourse.julialang.org/t/passing-a-local-to-cfunction-fails/40815):
>
> I want to pass a local function via @cfunction to a library, but can’t do it. (This function is a callback, and should store some data as it’s called). Minimal non-working example of what I would like: const badlocaldata = Int[] function badbar(i::Int) global badlocaldata push!(badlocaldata,i) end function foo() localdata = Int[] function bar(i::Int) push!(localdata,i) end # this fails with error message "UndefVarError: bar not defined" ccall(:libfunction,…

You may also find [this manual section](https://docs.julialang.org/en/v1/manual/calling-c-and-fortran-code/) useful.  
One solution for passing a local is to construct the `gsl_function` callback yourself like this:

```julia
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](https://juliamath.github.io/QuadGK.jl/stable/quadgk-examples/#Cauchy-principal-values). If you also notice that `1-fdist` causes [catastrophic cancellation](https://en.wikipedia.org/wiki/Catastrophic_cancellation) I would rewrite the problem like this:

```julia
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)

```

---

<div class="post-metadata">

**Author:** ![zxct](https://avatars.discourse-cdn.com/v4/letter/z/5daacb/32.png) [@zxct](https://discourse.julialang.org/u/zxct)\
**Post date:** [April 24, 2024, 2:05am UTC](https://discourse.julialang.org/t/how-to-use-integration-qawc-with-function-wrapped-with-local-variables/113269/3 "2024-04-24T02:05:09Z")

</div>

Thanks a lot for your reply. Now I know it’s the limitation of the C interface. I will adopt the QuadGK method.