# [Help] Nested ForwardDiff.jacobian of an Integrals.IntegralProblem

I have an integration problem (actually vector-valued, but problem here applies to scalar-valued ones as well) as a part of my model which it might be helpful to calculate its Hessian. I follow the documentation of Integrals.jl, however the nested Jacobian doesn’t seem to work for some solver algorithms. I suspect it’s due to the problem with propagating Dual types, but not sure, and hence also don’t know a priori which solvers are safe, and if there are work-around for the problematic solvers for this case. I would be happy if someone can help me to understand the issue here. Thank you!

Here is MWE, where `HCubatureJL` works fine but `CubatureJLh` doesn’t:

``````using Integrals, ForwardDiff
using IntegralsCubature
my_parameters = [1.0, 2.0]
my_function(x, p) = x^2 + p[1] * x + p[2]
function my_integration(p)
my_problem = IntegralProblem(my_function, -1.0, 1.0, p)
# return solve(my_problem, HCubatureJL(), reltol=1e-3, abstol=1e-3)  # Works
return solve(my_problem, CubatureJLh(), reltol=1e-3, abstol=1e-3)  # Errors
end
my_solution = my_integration(my_parameters)
my_jakobian = ForwardDiff.jacobian(my_integration, my_parameters)
ForwardDiff.jacobian(x->ForwardDiff.jacobian(my_integration, x), my_parameters)
``````

and the error:

``````ERROR: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6})
Closest candidates are:
(::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat at rounding.jl:200
(::Type{T})(::T) where T<:Number at boot.jl:772
(::Type{T})(::AbstractChar) where T<:Union{AbstractChar, Number} at char.jl:50
...
Stacktrace:
[1] (::Cubature.var"#17#18"{Bool, Bool, Int64, Float64, Float64, Int64, Int32, Ptr{Nothing}, Cubature.IntegrandData{IntegralsCubature.var"#3#15"{IntegralProblem{false, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}, 6}}, typeof(my_function), Float64, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}}}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Int64})()
@ Cubature ~/.julia/packages/Cubature/5zwuu/src/Cubature.jl:215
[2] disable_sigint
@ ./c.jl:473 [inlined]
[3] cubature(xscalar::Bool, fscalar::Bool, vectorized::Bool, padaptive::Bool, fdim::Int64, f::IntegralsCubature.var"#3#15"{IntegralProblem{false, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}, 6}}, typeof(my_function), Float64, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}}, xmin_::Float64, xmax_::Float64, reqRelError::Float64, reqAbsError::Float64, maxEval::Int64, error_norm::Int32)
@ Cubature ~/.julia/packages/Cubature/5zwuu/src/Cubature.jl:169
[4] hquadrature(f::Function, xmin::Float64, xmax::Float64; reltol::Float64, abstol::Float64, maxevals::Int64)
@ Cubature ~/.julia/packages/Cubature/5zwuu/src/Cubature.jl:230
[5] __solvebp_call(::IntegralProblem{false, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}, 6}}, typeof(my_function), Float64, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, ::CubatureJLh, ::Integrals.ReCallVJP{Integrals.ZygoteVJP}, ::Float64, ::Float64, ::Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}; reltol::Float64, abstol::Float64, maxiters::Int64, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
[6] __solvebp(prob::IntegralProblem{false, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}, 6}}, typeof(my_function), Float64, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, alg::CubatureJLh, sensealg::Integrals.ReCallVJP{Integrals.ZygoteVJP}, lb::Float64, ub::Float64, p::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}, 6}}; kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol, Symbol}, NamedTuple{(:reltol, :abstol), Tuple{Float64, Float64}}})
@ Integrals.IntegralsForwardDiffExt ~/.julia/packages/Integrals/9qNWp/ext/IntegralsForwardDiffExt.jl:23
[7] solve(prob::IntegralProblem{false, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}, 6}}, typeof(my_function), Float64, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, alg::CubatureJLh; sensealg::Integrals.ReCallVJP{Integrals.ZygoteVJP}, do_inf_transformation::Nothing, kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol, Symbol}, NamedTuple{(:reltol, :abstol), Tuple{Float64, Float64}}})
@ Integrals ~/.julia/packages/Integrals/9qNWp/src/Integrals.jl:108
[8] my_integration(p::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}, 6}})
@ Main ~/Documents/Research_local/Projects_local/SpatialCompetition/test/main.jl:27
[9] vector_mode_dual_eval!
@ ~/.julia/packages/ForwardDiff/vXysl/src/apiutils.jl:24 [inlined]
[10] vector_mode_jacobian(f::typeof(my_integration), x::Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}, 6, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}, 6}}})
@ ForwardDiff ~/.julia/packages/ForwardDiff/vXysl/src/jacobian.jl:125
[11] jacobian(f::Function, x::Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}, 6, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}, 6}}}, ::Val{true})
@ ForwardDiff ~/.julia/packages/ForwardDiff/vXysl/src/jacobian.jl:21
[12] jacobian(f::Function, x::Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}, 6, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}, 6}}}) (repeats 2 times)
@ ForwardDiff ~/.julia/packages/ForwardDiff/vXysl/src/jacobian.jl:19
[13] #11
@ ~/Documents/Research_local/Projects_local/SpatialCompetition/test/main.jl:31 [inlined]
[14] vector_mode_dual_eval!
@ ~/.julia/packages/ForwardDiff/vXysl/src/apiutils.jl:24 [inlined]
[15] vector_mode_jacobian(f::var"#11#12", x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}})
@ ForwardDiff ~/.julia/packages/ForwardDiff/vXysl/src/jacobian.jl:125
[16] jacobian(f::Function, x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}}, ::Val{true})
@ ForwardDiff ~/.julia/packages/ForwardDiff/vXysl/src/jacobian.jl:21
[17] jacobian(f::Function, x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}}) (repeats 2 times)
@ ForwardDiff ~/.julia/packages/ForwardDiff/vXysl/src/jacobian.jl:19
[18] top-level scope
@ ~/Documents/Research_local/Projects_local/SpatialCompetition/test/main.jl:31
``````

I didn’t run the code, but does it make a difference if the bounds are given as vectors, e.g. `IntegralProblem(my_function, -ones(1), ones(1), p)`? I had similar looking errors once in the past, but not sure…

Thank you @SteffenPL, but unfortunately not. In order to have vector bounds, I vectorized the integrand `my_function = @. x^2 + p[1] * x + p[2]` (since otherwise 1-element vector bounds are not unwrapped automatically), but I got the same error:

``````ERROR: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2})
Closest candidates are:
(::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat at rounding.jl:200
(::Type{T})(::T) where T<:Number at boot.jl:772
(::Type{T})(::AbstractChar) where T<:Union{AbstractChar, Number} at char.jl:50
...
Stacktrace:
[1] (::Cubature.var"#17#18"{Bool, Bool, Int64, Float64, Float64, Int64, Int32, Ptr{Nothing}, Cubature.IntegrandData{IntegralsCubature.var"#3#15"{IntegralProblem{false, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}, 2}}, typeof(my_function), Vector{Float64}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}}}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Int64})()
@ Cubature ~/.julia/packages/Cubature/5zwuu/src/Cubature.jl:215
[2] disable_sigint
@ ./c.jl:473 [inlined]
[3] cubature(xscalar::Bool, fscalar::Bool, vectorized::Bool, padaptive::Bool, fdim::Int64, f::IntegralsCubature.var"#3#15"{IntegralProblem{false, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}, 2}}, typeof(my_function), Vector{Float64}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}}, xmin_::Vector{Float64}, xmax_::Vector{Float64}, reqRelError::Float64, reqAbsError::Float64, maxEval::Int64, error_norm::Int32)
@ Cubature ~/.julia/packages/Cubature/5zwuu/src/Cubature.jl:169
[4] hcubature(f::Function, xmin::Vector{Float64}, xmax::Vector{Float64}; reltol::Float64, abstol::Float64, maxevals::Int64)
@ Cubature ~/.julia/packages/Cubature/5zwuu/src/Cubature.jl:230
[5] __solvebp_call(::IntegralProblem{false, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}, 2}}, typeof(my_function), Vector{Float64}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, ::CubatureJLh, ::Integrals.ReCallVJP{Integrals.ZygoteVJP}, ::Vector{Float64}, ::Vector{Float64}, ::Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}; reltol::Float64, abstol::Float64, maxiters::Int64, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
[6] __solvebp(prob::IntegralProblem{false, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}, 2}}, typeof(my_function), Vector{Float64}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, alg::CubatureJLh, sensealg::Integrals.ReCallVJP{Integrals.ZygoteVJP}, lb::Vector{Float64}, ub::Vector{Float64}, p::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}, 2}}; kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol, Symbol}, NamedTuple{(:reltol, :abstol), Tuple{Float64, Float64}}})
@ Integrals.IntegralsForwardDiffExt ~/.julia/packages/Integrals/9qNWp/ext/IntegralsForwardDiffExt.jl:23
[7] #solve#34
@ ~/.julia/packages/Integrals/9qNWp/src/Integrals.jl:108 [inlined]
[8] my_integration(p::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}, 2}})
@ Main ~/Documents/Research_local/Projects_local/SpatialCompetition/test/main.jl:27
[9] vector_mode_dual_eval!
@ ~/.julia/packages/ForwardDiff/vXysl/src/apiutils.jl:24 [inlined]
[10] vector_mode_jacobian(f::typeof(my_integration), x::Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}, 2, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}, 2}}})
@ ForwardDiff ~/.julia/packages/ForwardDiff/vXysl/src/jacobian.jl:125
[11] jacobian(f::Function, x::Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}, 2, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}, 2}}}, ::Val{true})
@ ForwardDiff ~/.julia/packages/ForwardDiff/vXysl/src/jacobian.jl:21
[12] jacobian(f::Function, x::Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}, 2, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}, 2}}}) (repeats 2 times)
@ ForwardDiff ~/.julia/packages/ForwardDiff/vXysl/src/jacobian.jl:19
[13] #5
@ ~/Documents/Research_local/Projects_local/SpatialCompetition/test/main.jl:31 [inlined]
[14] vector_mode_dual_eval!
@ ~/.julia/packages/ForwardDiff/vXysl/src/apiutils.jl:24 [inlined]
[15] vector_mode_jacobian(f::var"#5#6", x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}})
@ ForwardDiff ~/.julia/packages/ForwardDiff/vXysl/src/jacobian.jl:125
[16] jacobian(f::Function, x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}}, ::Val{true})
@ ForwardDiff ~/.julia/packages/ForwardDiff/vXysl/src/jacobian.jl:21
[17] jacobian(f::Function, x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}}) (repeats 2 times)
@ ForwardDiff ~/.julia/packages/ForwardDiff/vXysl/src/jacobian.jl:19
[18] top-level scope
@ ~/Documents/Research_local/Projects_local/SpatialCompetition/test/main.jl:31
``````
1 Like