Multidimentional integrals in julia: speed and automatic differentiation

Specific Domains Numerics
6

With Zygote, you could solve it like this, by differentiating through the integrand:

using Cuba, Zygote

function model(cosθ,ϕ; αp1 = 2.3)
    α1 = α2 = αp1-1.0
    v = (exp(-α1*(cosθ+1)) + exp(-α2*(1-cosθ)))*(1+sin(ϕ))^(1.2)
    return abs2(v)
end

integrand(x, αp1) = model(2x[1]-1, π*(2x[2]-1);αp1=αp1)
I_cuba(αp1) = cuhre((x, f) -> f[1]=integrand(x, αp1),2,1).integral[1]*(4π)
Zygote.@adjoint I_cuba(αp1) = I_cuba(αp1), Δ -> (Δ * cuhre((x, f) -> f[1]=Zygote.pullback(a -> integrand(x, a), αp1)[2](1)[1], 2, 1).integral[1] * (4π),)

And then get the gradient:

julia> Zygote.gradient(x -> I_cuba(x[1]), [2.3])
([-13.368478696191737],)

It should be fairly easy to replace αp1 with a vector if you have multiple parameters.

2 Likes