Forward differentiation and differential equations

moesphere · February 20, 2019, 7:08am

I am trying to use forward AD in combination with solving an ODE using DifferentialEquations.jl. I want to use an own loss function. However, with the is-in-place version there is a type assertion error. I have come up with following minimal example:

using DifferentialEquations, Optim, LinearAlgebra

function f(dxdt,x,p,t)
    dxdt[1] = -p[1]*x[1]
    dxdt[2] = p[1]*x[1]
end

tspan = (0.,8.)
tdata = range(0.,stop=8.,length=10)
p0 = [1.]
x0 = [10.,0.]

prob = ODEProblem(f,x0,tspan,p0)

sol = DifferentialEquations.solve(prob; saveat=tdata)

xdata = sol[2,:] + randn(length(sol[2,:])) # in-silico data

function loss_function(x)
    prob_new = remake(prob; p=x)
    sol_new = DifferentialEquations.solve(prob_new; saveat=tdata)
    norm(sol_new[2,:] - xdata)^2
end
optimize(loss_function, [2.0], LBFGS(),autodiff=:forward)

Error message is TypeError: in typeassert, expected Float64, got ForwardDiff.Dual{Nothing,Float64,1}.

Is there a solution to this? It works fine when the autodiff keyword is removed, i.e. using finite differences.

pkofod · February 20, 2019, 10:41am

Why would you rob people of the whole stack trace?

julia> optimize(loss_function, [2.0], LBFGS(),autodiff=:forward)
ERROR: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1})
Closest candidates are:
  Float64(::Real, ::RoundingMode) where T<:AbstractFloat at rounding.jl:185
  Float64(::T<:Number) where T<:Number at boot.jl:725
  Float64(::Int8) at float.jl:60
  ...
Stacktrace:
 [1] convert(::Type{Float64}, ::ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1}) at ./number.jl:7
 [2] setindex!(::Array{Float64,1}, ::ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1}, ::Int64) at ./array.jl:769
 [3] f(::Array{Float64,1}, ::Array{Float64,1}, ::Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1}, ::Float64) at ./REPL[15]:2
 [4] (::ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing})(::Array{Float64,1}, ::Array{Float64,1}, ::Vararg{Any,N} where N) at /home/pkofod/.julia/packages/DiffEqBase/tENmC/src/diffeqfunction.jl:107
 [5] initialize!(::OrdinaryDiffEq.ODEIntegrator{CompositeAlgorithm{Tuple{Tsit5,Rosenbrock23{0,false,LinSolveFactorize{typeof(lu!)},DataType}},AutoSwitch{Tsit5,Rosenbrock23{0,false,LinSolveFactorize{typeof(lu!)},DataType},Rational{Int64},Float64}},true,Array{Float64,1},Float64,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1},Float64,Float64,Float64,Array{Array{Float64,1},1},OrdinaryDiffEq.ODECompositeSolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},Array{Array{Array{Float64,1},1},1},ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem},CompositeAlgorithm{Tuple{Tsit5,Rosenbrock23{0,false,LinSolveFactorize{typeof(lu!)},DataType}},AutoSwitch{Tsit5,Rosenbrock23{0,false,LinSolveFactorize{typeof(lu!)},DataType},Rational{Int64},Float64}},OrdinaryDiffEq.CompositeInterpolationData{ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}},OrdinaryDiffEq.Rosenbrock23Cache{Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Rosenbrock23ConstantCache{Float64,typeof(identity),typeof(identity)},DiffEqDiffTools.TimeGradientWrapper{ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1}},DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1}},LinSolveFactorize{typeof(lu!)},DiffEqDiffTools.JacobianCache{Array{Float64,1},Array{Float64,1},Array{Float64,1},Val{:central},Float64,Val{true}},DiffEqDiffTools.GradientCache{Nothing,Array{Float64,1},Array{Float64,1},Val{:central},Float64,Val{true}}}},AutoSwitch{Tsit5,Rosenbrock23{0,false,LinSolveFactorize{typeof(lu!)},DataType},Rational{Int64},Float64}}}},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}},OrdinaryDiffEq.Rosenbrock23Cache{Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Rosenbrock23ConstantCache{Float64,typeof(identity),typeof(identity)},DiffEqDiffTools.TimeGradientWrapper{ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1}},DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1}},LinSolveFactorize{typeof(lu!)},DiffEqDiffTools.JacobianCache{Array{Float64,1},Array{Float64,1},Array{Float64,1},Val{:central},Float64,Val{true}},DiffEqDiffTools.GradientCache{Nothing,Array{Float64,1},Array{Float64,1},Val{:central},Float64,Val{true}}}},AutoSwitch{Tsit5,Rosenbrock23{0,false,LinSolveFactorize{typeof(lu!)},DataType},Rational{Int64},Float64}},OrdinaryDiffEq.DEOptions{Float64,Float64,Float64,Float64,typeof(DiffEqBase.ODE_DEFAULT_NORM),typeof(opnorm),CallbackSet{Tuple{},Tuple{}},typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN),typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE),typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK),DataStructures.BinaryHeap{Float64,DataStructures.LessThan},DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Nothing,Nothing,Int64,Array{Float64,1},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},Array{Float64,1}},Array{Float64,1},Float64}, ::OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}}) at /home/pkofod/.julia/packages/OrdinaryDiffEq/HJOah/src/perform_step/low_order_rk_perform_step.jl:580
 [6] initialize!(::OrdinaryDiffEq.ODEIntegrator{CompositeAlgorithm{Tuple{Tsit5,Rosenbrock23{0,false,LinSolveFactorize{typeof(lu!)},DataType}},AutoSwitch{Tsit5,Rosenbrock23{0,false,LinSolveFactorize{typeof(lu!)},DataType},Rational{Int64},Float64}},true,Array{Float64,1},Float64,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1},Float64,Float64,Float64,Array{Array{Float64,1},1},OrdinaryDiffEq.ODECompositeSolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},Array{Array{Array{Float64,1},1},1},ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem},CompositeAlgorithm{Tuple{Tsit5,Rosenbrock23{0,false,LinSolveFactorize{typeof(lu!)},DataType}},AutoSwitch{Tsit5,Rosenbrock23{0,false,LinSolveFactorize{typeof(lu!)},DataType},Rational{Int64},Float64}},OrdinaryDiffEq.CompositeInterpolationData{ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}},OrdinaryDiffEq.Rosenbrock23Cache{Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Rosenbrock23ConstantCache{Float64,typeof(identity),typeof(identity)},DiffEqDiffTools.TimeGradientWrapper{ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1}},DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1}},LinSolveFactorize{typeof(lu!)},DiffEqDiffTools.JacobianCache{Array{Float64,1},Array{Float64,1},Array{Float64,1},Val{:central},Float64,Val{true}},DiffEqDiffTools.GradientCache{Nothing,Array{Float64,1},Array{Float64,1},Val{:central},Float64,Val{true}}}},AutoSwitch{Tsit5,Rosenbrock23{0,false,LinSolveFactorize{typeof(lu!)},DataType},Rational{Int64},Float64}}}},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}},OrdinaryDiffEq.Rosenbrock23Cache{Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Rosenbrock23ConstantCache{Float64,typeof(identity),typeof(identity)},DiffEqDiffTools.TimeGradientWrapper{ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1}},DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1}},LinSolveFactorize{typeof(lu!)},DiffEqDiffTools.JacobianCache{Array{Float64,1},Array{Float64,1},Array{Float64,1},Val{:central},Float64,Val{true}},DiffEqDiffTools.GradientCache{Nothing,Array{Float64,1},Array{Float64,1},Val{:central},Float64,Val{true}}}},AutoSwitch{Tsit5,Rosenbrock23{0,false,LinSolveFactorize{typeof(lu!)},DataType},Rational{Int64},Float64}},OrdinaryDiffEq.DEOptions{Float64,Float64,Float64,Float64,typeof(DiffEqBase.ODE_DEFAULT_NORM),typeof(opnorm),CallbackSet{Tuple{},Tuple{}},typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN),typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE),typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK),DataStructures.BinaryHeap{Float64,DataStructures.LessThan},DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Nothing,Nothing,Int64,Array{Float64,1},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},Array{Float64,1}},Array{Float64,1},Float64}, ::OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}},OrdinaryDiffEq.Rosenbrock23Cache{Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Rosenbrock23ConstantCache{Float64,typeof(identity),typeof(identity)},DiffEqDiffTools.TimeGradientWrapper{ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1}},DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1}},LinSolveFactorize{typeof(lu!)},DiffEqDiffTools.JacobianCache{Array{Float64,1},Array{Float64,1},Array{Float64,1},Val{:central},Float64,Val{true}},DiffEqDiffTools.GradientCache{Nothing,Array{Float64,1},Array{Float64,1},Val{:central},Float64,Val{true}}}},AutoSwitch{Tsit5,Rosenbrock23{0,false,LinSolveFactorize{typeof(lu!)},DataType},Rational{Int64},Float64}}) at /home/pkofod/.julia/packages/OrdinaryDiffEq/HJOah/src/perform_step/composite_perform_step.jl:38
 [7] #__init#303(::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}, ::Array{Float64,1}, ::Array{Float64,1}, ::Nothing, ::Bool, ::Nothing, ::Bool, ::Bool, ::Bool, ::Nothing, ::Bool, ::Bool, ::Float64, ::Bool, ::Rational{Int64}, ::Nothing, ::Nothing, ::Int64, ::Rational{Int64}, ::Int64, ::Int64, ::Rational{Int64}, ::Bool, ::Int64, ::Nothing, ::Nothing, ::Int64, ::Float64, ::Float64, ::typeof(DiffEqBase.ODE_DEFAULT_NORM), ::typeof(opnorm), ::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), ::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Int64, ::String, ::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), ::Nothing, ::Bool, ::Bool, ::Bool, ::Base.Iterators.Pairs{Symbol,Bool,Tuple{Symbol},NamedTuple{(:default_set,),Tuple{Bool}}}, ::typeof(DiffEqBase.__init), ::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}, ::CompositeAlgorithm{Tuple{Tsit5,Rosenbrock23{0,false,LinSolveFactorize{typeof(lu!)},DataType}},AutoSwitch{Tsit5,Rosenbrock23{0,false,LinSolveFactorize{typeof(lu!)},DataType},Rational{Int64},Float64}}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /home/pkofod/.julia/packages/OrdinaryDiffEq/HJOah/src/solve.jl:273
 [8] (::getfield(DiffEqBase, Symbol("#kw##__init")))(::NamedTuple{(:default_set, :saveat),Tuple{Bool,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}, ::typeof(DiffEqBase.__init), ::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}, ::CompositeAlgorithm{Tuple{Tsit5,Rosenbrock23{0,false,LinSolveFactorize{typeof(lu!)},DataType}},AutoSwitch{Tsit5,Rosenbrock23{0,false,LinSolveFactorize{typeof(lu!)},DataType},Rational{Int64},Float64}}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at ./none:0
 [9] #__solve#302(::Base.Iterators.Pairs{Symbol,Any,Tuple{Symbol,Symbol},NamedTuple{(:default_set, :saveat),Tuple{Bool,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}}, ::Function, ::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}, ::CompositeAlgorithm{Tuple{Tsit5,Rosenbrock23{0,false,LinSolveFactorize{typeof(lu!)},DataType}},AutoSwitch{Tsit5,Rosenbrock23{0,false,LinSolveFactorize{typeof(lu!)},DataType},Rational{Int64},Float64}}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /home/pkofod/.julia/packages/OrdinaryDiffEq/HJOah/src/solve.jl:6
 [10] (::getfield(DiffEqBase, Symbol("#kw##__solve")))(::NamedTuple{(:default_set, :saveat),Tuple{Bool,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}, ::typeof(DiffEqBase.__solve), ::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}, ::CompositeAlgorithm{Tuple{Tsit5,Rosenbrock23{0,false,LinSolveFactorize{typeof(lu!)},DataType}},AutoSwitch{Tsit5,Rosenbrock23{0,false,LinSolveFactorize{typeof(lu!)},DataType},Rational{Int64},Float64}}) at ./none:0
 [11] #__solve#2(::Bool, ::Base.Iterators.Pairs{Symbol,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},Tuple{Symbol},NamedTuple{(:saveat,),Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}}, ::Function, ::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}, ::Nothing) at /home/pkofod/.julia/packages/DifferentialEquations/UZIcy/src/default_solve.jl:15
 [12] (::getfield(DiffEqBase, Symbol("#kw##solve")))(::NamedTuple{(:saveat,),Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}, ::typeof(solve), ::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}) at ./none:0
 [13] loss_function(::Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1}) at ./REPL[23]:3
 [14] vector_mode_gradient!(::DiffResults.MutableDiffResult{1,Float64,Tuple{Array{Float64,1}}}, ::typeof(loss_function), ::Array{Float64,1}, ::ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1}}) at /home/pkofod/.julia/packages/ForwardDiff/okZnq/src/apiutils.jl:35
 [15] gradient! at /home/pkofod/.julia/packages/ForwardDiff/okZnq/src/gradient.jl:35 [inlined]
 [16] gradient! at /home/pkofod/.julia/packages/ForwardDiff/okZnq/src/gradient.jl:33 [inlined]
 [17] (::getfield(NLSolversBase, Symbol("##14#18")){Float64,typeof(loss_function),ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1}}})(::Array{Float64,1}, ::Array{Float64,1}) at /home/pkofod/.julia/packages/NLSolversBase/LRo6O/src/objective_types/oncedifferentiable.jl:65
 [18] value_gradient!!(::OnceDifferentiable{Float64,Array{Float64,1},Array{Float64,1}}, ::Array{Float64,1}) at /home/pkofod/.julia/packages/NLSolversBase/LRo6O/src/interface.jl:82
 [19] initial_state(::LBFGS{Nothing,LineSearches.InitialStatic{Float64},LineSearches.HagerZhang{Float64,Base.RefValue{Bool}},getfield(Optim, Symbol("##22#24"))}, ::Optim.Options{Float64,Nothing}, ::OnceDifferentiable{Float64,Array{Float64,1},Array{Float64,1}}, ::Array{Float64,1}) at /home/pkofod/.julia/dev/Optim/src/multivariate/solvers/first_order/l_bfgs.jl:158
 [20] optimize(::OnceDifferentiable{Float64,Array{Float64,1},Array{Float64,1}}, ::Array{Float64,1}, ::LBFGS{Nothing,LineSearches.InitialStatic{Float64},LineSearches.HagerZhang{Float64,Base.RefValue{Bool}},getfield(Optim, Symbol("##22#24"))}, ::Optim.Options{Float64,Nothing}) at /home/pkofod/.julia/dev/Optim/src/multivariate/optimize/optimize.jl:33
 [21] #optimize#87 at /home/pkofod/.julia/dev/Optim/src/multivariate/optimize/interface.jl:114 [inlined]
 [22] (::getfield(Optim, Symbol("#kw##optimize")))(::NamedTuple{(:autodiff,),Tuple{Symbol}}, ::typeof(optimize), ::Function, ::Array{Float64,1}, ::LBFGS{Nothing,LineSearches.InitialStatic{Float64},LineSearches.HagerZhang{Float64,Base.RefValue{Bool}},getfield(Optim, Symbol("##22#24"))}, ::Optim.Options{Float64,Nothing}) at ./none:0 (repeats 2 times)
 [23] top-level scope at none:0

So without knowing the internals of ODEFunction it’s probably because there is some convert (attempt to take the input and make sure it has a specific type for the rest of the code to be type stable) that is based of off some parametric type. This gets initialized in your constructor, so it’s forced to be Float64 at the time of ForwardDiff’ing - causing the error you see. @ChrisRackauckas do you recognize this problem with ODEProblem?

edit: oooh… the stack trace (though it has a vertical bar) is much easier to read without line wrapping

ChrisRackauckas · February 20, 2019, 11:17am

The automatic differentiation page discusses exactly this issue:

http://docs.juliadiffeq.org/latest/analysis/sensitivity.html#Examples-using-ForwardDiff.jl-1

If you want to propagate derivatives in state, you need to make sure your state variables are dual numbers as well, so you’re missing the type conversion of u0 to be the same as p. Your remake should be:

_prob = remake(prob;u0=convert.(eltype(p),prob.u0),p=p)

moesphere · February 20, 2019, 11:58am

Thanks for the answer!

I should have seen this in the docs.

Matthew_Overlin · September 28, 2020, 1:45pm

I am getting this issue with delay differential equations problem. I define: dde_problem = DDEProblem(…); Then, I remake the problem as: tmp_dde_problem = remake(dde_problem; u0=convert.(eltype(my_parameters), dde_problem.u0),p=my_parameters); Then, I try to solve tmp_dde_problem and I get an error of: ERROR: LoadError: TypeError: in typeassert, expected Float64, got a value of type ForwardDiff.Dual{Nothing,Float64,9}

I would expect similar behavior when ForwardDiff is used with an ODEProblem and a DDEProblem. I have had success using ForwardDiff with ODEProblem. Could anyone provide insights as to why it might not be working for DDEProblem?
Thanks.

ChrisRackauckas · September 28, 2020, 1:47pm

I’d need more details since this is covered by tests

https://github.com/SciML/DelayDiffEq.jl/blob/master/test/interface/ad.jl

pierre-haessig · August 31, 2022, 2:48pm

Just searching for infos on the topic of ForwardDiff and type issues, I came across this thread and noticed that the documentation link to juliadiffeq.org is broken. Here are updates:

Archived version of late 2019 Local Sensitivity Analysis (Automatic Differentiation) · DifferentialEquations.jl
The closest document on the new sciml.ai website: Local Sensitivity Analysis (Automatic Differentiation) · DifferentialEquations.jl

However, in the new doc, I see no more any usage of convert inside the call to remake.

ChrisRackauckas · August 31, 2022, 2:50pm

It’s handled automatically these days.

edljk · November 5, 2023, 8:27pm

A related question: as mentioned above and in the ode faq, automatic differentiation of the final state with respect to parameters using ForwardDiff works perfectly.

I am interested in computing all the jacobians associated to many final states with respect to parameters but with different initial states (with the same parameters).

I wonder if there is a way to keep the benefit of the multi-threading provided by an EnsembleProblem in that context?

Thanks.

ChrisRackauckas · November 5, 2023, 8:34pm

Kind of. It’s a bit easier to just use PolyesterForwardDiff if you just want multithreading though, but yeah if you build the duals yourself you can make use of ensembles.

edljk · November 5, 2023, 9:09pm

Thanks for the fast answer. Since I am not enough familiar with ForwardDiff and duals I am gonna try PolyesterForwardDiff.jl which looks great!
Btw, in my setting, PolyesterForwardDiff.jl seems to downgrade many packages like:

⌅ [aae01518] ↓ BandedMatrices v1.1.0 ⇒ v0.17.18
⌃ [052768ef] ↓ CUDA v5.0.0 ⇒ v4.4.1
⌃ [13f3f980] ↓ CairoMakie v0.10.12 ⇒ v0.9.4
⌅ [34da2185] ↓ Compat v4.10.0 ⇒ v3.46.2
⌃ [a93c6f00] ↓ DataFrames v1.6.1 ⇒ v1.3.6
⌃ [0c46a032] ↓ DifferentialEquations v7.11.0 ⇒ v7.2.0
⌃ [e9467ef8] ↓ GLMakie v0.8.12 ⇒ v0.7.4
⌅ [7ed4a6bd] ↓ LinearSolve v2.15.0 ⇒ v1.34.1
…

Is it expected?

ChrisRackauckas · November 5, 2023, 9:14pm

edljk · November 5, 2023, 11:18pm

Probably some difficulties from previous update:

num_threads() is not defined anymore
Defining num_threads() = Threads.nthreads() and adapting the example of the ode faq:

using DifferentialEquations
function func(du, u, p, t)
    du[1] = p[1] * u[1] - p[2] * u[1] * u[2]
    du[2] = -3 * u[2] + u[1] * u[2]
end
function f(p)
    prob = ODEProblem(func, eltype(p).([1.0, 1.0]), (0.0, 10.0), p)
    # Lower tolerances to show the methods converge to the same value
    solve(prob, Tsit5(), save_everystep = false, abstol = 1e-12, reltol = 1e-12)[end]
end
J = zeros(2, 2)
using PolyesterForwardDiff
PolyesterForwardDiff.threaded_jacobian!(f, J, [1.5, 1.0], ForwardDiff.Chunk(8))

gives the error

ERROR: Cannot determine ordering of Dual tags ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, ForwardDiff.Dual{Nothing, Float64, 8}} and Nothing
Stacktrace:
  [1] ≺(a::Type, b::Type)
...

edljk · November 7, 2023, 9:01am

I finally had a try with Dual evaluations using EnsembleProblem and following ForwardDiff.jacobian implementation.
The naive implementation below seems to work on this simple example.
Any comments are welcome.

using DifferentialEquations, ForwardDiff
using TimerOutputs

to_ode = TimerOutput()
#-------------------------------------------------------------------------------
function func(du, u, p, t)
    du[1] = p[1] * u[1] - p[2] * u[1] * u[2]
    du[2] = -3 * u[2] + u[1] * u[2]
end
function f(u0, λ, εode)
    prob = ODEProblem(func, u0, (0.0, 10.0), λ)
    solve(prob, Tsit5(), save_everystep = false, abstol = εode, 
          reltol = εode)[end]
end
#-------------------------------------------------------------------------------
function jacobiansODE(f, prob, Vu0, λ, solver,
    cfg::ForwardDiff.JacobianConfig{T} = ForwardDiff.JacobianConfig(f, λ);
    εode::Float64 = 1e-4, maxit::Int64 = 1_000) where T

    N = ForwardDiff.chunksize(cfg)
    Vjacobians, Vydual = Array[], Vector{ForwardDiff.Dual}[]

    λdual = cfg.duals
    ForwardDiff.seed!(λdual, λ, cfg.seeds)

    # use multi-threading to evaluate dual variables
    function prob_funcdual(prob, i, repeat)
        remake(prob, u0 = convert.(eltype(λdual), Vu0[i]), p = λdual)
    end
    ∇ensemble_prob = EnsembleProblem(prob, prob_func = prob_funcdual)
    sols = solve(∇ensemble_prob, solver, trajectories = length(Vu0),
                 maxiters = maxit, abstol = εode, reltol = εode,  
                 save_everystep = false)
    for sol ∈ sols 
        push!(Vydual, sol[end])
    end

    for ydual ∈ Vydual
        ydual isa AbstractArray || throw(ForwardDiff.JACOBIAN_ERROR)
        result = similar(ydual, ForwardDiff.valtype(eltype(ydual)), 
                        length(ydual), N)
        ForwardDiff.extract_jacobian!(T, result, ydual, N)
        ForwardDiff.extract_value!(T, result, ydual)
        push!(Vjacobians, [result[k, l].value for k = 1:size(result, 1), 
                                                  l = 1:size(result, 2)])
    end
    return Vjacobians
end
#-------------------------------------------------------------------------------
function testjacobians(nbexp::Int64 = 10)
    εode = 1e-8
    λ = [1.5, 1.0]
    allu0 = [(2 * rand(2) .- 1) / 10 .+ 1 for _ = 1:nbexp]
    prob = ODEProblem(func, allu0[1], (0.0, 10.0), λ)
    jacobians1 = Matrix[]
    for u0 ∈ allu0
        @timeit to_ode "ForwardDiff" push!(jacobians1, 
                                   ForwardDiff.jacobian(p -> f(u0, p, εode), λ))
    end
    @timeit to_ode "jacobiansODE" jacobians2 = jacobiansODE(f, prob, allu0, λ, 
                                                        Tsit5(), εode = εode)
    @show norm(jacobians1 - jacobians2, Inf)

    print(to_ode)
end

# call test function 
testjacobians(10)