Getting gradient of taylorinteg

General Usage
1

Hi,

I’m quite new to Julia, and currently trying out the gradient provided by ForwardDiff.jl. I tried using it on a function that calls taylorinteg from the TaylorIntegration.jl package ; I was expecting it to give a weird result, or no result at all, but I get an error I’m unable to solve. Is it impossible or useless to try to get a result this way ?

Here is a MWE :

using TaylorIntegration, ForwardDiff

function deriv_x!(dx, x, params, t)
	dx = x
end

function f(x_start)
	t, x = taylorinteg(deriv_x!, x_start, 0.0, 10.,
		15, 1.0e-15, maxsteps=10_000)
	return x[end,:]
end

∇f(x) = ForwardDiff.jacobian(f, x)
print(∇f([1, 2, 3]))

And here is the error message :

ERROR: MethodError: no method matching setindex!(::Taylor1{ForwardDiff.Dual{ForwardDiff.Tag{typeof(f), Int64}, Int64, 3}}, ::ForwardDiff.Dual{ForwardDiff.Tag{typeof(f), Int64}, Float64, 3}, ::Int64)
Closest candidates are:
  setindex!(::Taylor1{T}, ::T, ::Int64) where T<:Number at /home/alseidon/.julia/packages/TaylorSeries/vTIIX/src/auxiliary.jl:86
  setindex!(::Taylor1{T}, ::Vector{T}, ::StepRange{Int64, Int64}) where T<:Number at /home/alseidon/.julia/packages/TaylorSeries/vTIIX/src/auxiliary.jl:99
  setindex!(::Taylor1{T}, ::T, ::StepRange{Int64, Int64}) where T<:Number at /home/alseidon/.julia/packages/TaylorSeries/vTIIX/src/auxiliary.jl:97
  ...
Stacktrace:
  [1] jetcoeffs!(eqsdiff!::typeof(deriv_x!), t::Taylor1{Float64}, x::Vector{Taylor1{ForwardDiff.Dual{ForwardDiff.Tag{typeof(f), Int64}, Int64, 3}}}, dx::Vector{Taylor1{ForwardDiff.Dual{ForwardDiff.Tag{typeof(f), Int64}, Int64, 3}}}, xaux::Vector{Taylor1{ForwardDiff.Dual{ForwardDiff.Tag{typeof(f), Int64}, Int64, 3}}}, params::Nothing)
    @ TaylorIntegration ~/.julia/packages/TaylorIntegration/4IefH/src/explicitode.jl:82
  [2] __jetcoeffs!(#unused#::Val{false}, f::Function, t::Taylor1{Float64}, x::Vector{Taylor1{ForwardDiff.Dual{ForwardDiff.Tag{typeof(f), Int64}, Int64, 3}}}, dx::Vector{Taylor1{ForwardDiff.Dual{ForwardDiff.Tag{typeof(f), Int64}, Int64, 3}}}, xaux::Vector{Taylor1{ForwardDiff.Dual{ForwardDiff.Tag{typeof(f), Int64}, Int64, 3}}}, params::Nothing)
    @ TaylorIntegration ~/.julia/packages/TaylorIntegration/4IefH/src/explicitode.jl:101
  [3] taylorstep!(f!::Function, t::Taylor1{Float64}, x::Vector{Taylor1{ForwardDiff.Dual{ForwardDiff.Tag{typeof(f), Int64}, Int64, 3}}}, dx::Vector{Taylor1{ForwardDiff.Dual{ForwardDiff.Tag{typeof(f), Int64}, Int64, 3}}}, xaux::Vector{Taylor1{ForwardDiff.Dual{ForwardDiff.Tag{typeof(f), Int64}, Int64, 3}}}, abstol::Float64, params::Nothing, parse_eqs::Bool)
    @ TaylorIntegration ~/.julia/packages/TaylorIntegration/4IefH/src/explicitode.jl:237
  [4] taylorinteg(f!::Function, q0::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(f), Int64}, Int64, 3}}, t0::Float64, tmax::Float64, order::Int64, abstol::Float64, params::Nothing; maxsteps::Int64, parse_eqs::Bool)
    @ TaylorIntegration ~/.julia/packages/TaylorIntegration/4IefH/src/explicitode.jl:478
  [5] f
    @ ./REPL[3]:2 [inlined]
  [6] vector_mode_dual_eval!
    @ ~/.julia/packages/ForwardDiff/PBzup/src/apiutils.jl:37 [inlined]
  [7] vector_mode_jacobian(f::typeof(f), x::Vector{Int64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(f), Int64}, Int64, 3, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(f), Int64}, Int64, 3}}})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/PBzup/src/jacobian.jl:148
  [8] jacobian(f::Function, x::Vector{Int64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(f), Int64}, Int64, 3, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(f), Int64}, Int64, 3}}}, ::Val{true})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/PBzup/src/jacobian.jl:21
  [9] jacobian(f::Function, x::Vector{Int64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(f), Int64}, Int64, 3, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(f), Int64}, Int64, 3}}}) (repeats 2 times)
    @ ForwardDiff ~/.julia/packages/ForwardDiff/PBzup/src/jacobian.jl:19
 [10] ∇f(x::Vector{Int64})
    @ Main ./REPL[4]:1
 [11] top-level scope
    @ REPL[5]:1

when I paste this in the REPL.

Thanks !

EDIT : changed ForwardDiff.gradient to ForwardDiff.jacobian, as we are in 2 dimensions.

2

Learned this recently from @ChrisRackauckas :

using TaylorIntegration, ForwardDiff

function deriv_x!(dx, x, params, t)
	dx = x
end

function f(x_start)
	t, x = taylorinteg(deriv_x!, x_start, convert(eltype(x_start), 0.0), convert(eltype(x_start), 10.),
		15, 1.0e-15, maxsteps=10_000)
	return x[end,:]
end

∇f(x) = ForwardDiff.jacobian(f, x)
print(∇f([1, 2, 3]))
4 Likes
3

Thanks, this works indeed ! I will now have to tweak my parameters in order to get rid of the Maximum number of integration steps reached; exiting. warning, but that’s another story. Do you know why setindex! is given the wrong type of variables in the first place ?

4

The problem was that you were writing a Dual number (generated through differentiation with ForwardDiff) to an array, dx, which was allocated as an array of Float64 inside taylorinteg, as the inputs were Float64s. The solution was to change the inputs to be such that an array of Duals is created instead so we can write Dual numbers, and hence derivatives, to it.

2 Likes
5

OK it is much clearer now, thanks !

1 Like
6

I believe taylorinteg could be improved to handle it automatically. But you would have to file an issue.

BTW: this was the thread I forgot to link.

1 Like