Cannot differentiate through `cuhre`. Explicit type inference

[mmikhasenko]

I love the library and I don't any other better way in Julia to do >1 dim integrals.

Just discovered an issed trying calculate gradient of the integrated function

using Cuba
using ForwardDiff
#
b(cosθ,ϕ,c) = c[1]*cosθ*sin(ϕ) + c[2]*cosθ^2*cos(ϕ)^2
f(c) = cuhre((x,f)->f[1]=b(2x[1]-1,π*(2x[2]-1),c),2,1)
ForwardDiff.gradient(c->f(c).integral[1], [1,1.1])

giving the error

MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{typeof(χ2),Float64},Float64,3})
Closest candidates are:
  Float64(::Real, !Matched::RoundingMode) where T<:AbstractFloat at rounding.jl:200
  Float64(::T) where T<:Number at boot.jl:718
  Float64(!Matched::Int8) at float.jl:60
  ...

most likely at dointegral function https://github.com/giordano/Cuba.jl/blob/3681ed83b8fa0fafba03c7d67d2689d0028e920c/src/Cuba.jl#L181

Interestingly a gradient of 1dim integral taken with QuadGK works (well, native Julia).

[mmikhasenko]

Found way out (not solution, unfortunately) checking references at your README (thanks)

using  ForwardDiff
using HCubature
b(cosθ,ϕ,c) = c[1]*cosθ*sin(ϕ) + c[2]*cosθ^2*cos(ϕ)^2
f(c) = hcubature(x->b(x[1],x[2],c), [-1.0,-π], [1.0,π])[1]
ForwardDiff.gradient(f, [1.1,1.1])

I am curious still if it is in principle possible to ForwardDiff wrappers with ccalls?

[mmikhasenko]

Ah, it seems that I still need Cuba.jl :) For my integrals, HCubature takes 15 times longer

34.546465 seconds (239.96 M allocations: 5.404 GiB, 3.64% gc time)

with Cuba.jl:

 1.424789 seconds (13.52 M allocations: 316.171 MiB, 4.57% gc time)

[giordano]

I am curious still if it is in principle possible to ForwardDiff wrappers with ccalls?

I'm not entirely sure, but I fear not :disappointed: I'd have suggested you to use a pure-Julia library like HCubature.jl or NIntegration.jl

[mmikhasenko]

If we are sure that there are no tricks to pull Dual through c++ implementation, the issue can be closed.

It seems to me that it is not possible, running the same code on different types is one of the main advantages of julia. It does not work so naturally with external calls.

[giordano]

Honestly I don't think there is anything we can do here in terms of automatic differentiation, so I'm going to close this issue.

PS: the Cuba library is written in C, not C++ :wink:

[mmikhasenko]

I copy here the reply of @simeonschaub on discourse for completeness. He provides an example of differentiation with Cuba.jl!

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