forward differentiation using `optimize` for DDE

[Farnazmdi]

Hello,

I am working on an optimization problem for delay differential equation, and I wanted to use forward differentiation using optimize(...; autodif:=forward), but it doesn't work. This is the error we get: MethodError: Cannotconvert` an object of type Discontinuity{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("#residuals#13")){Int64,Array{Float64,2},Array{Float64,2},Array{Float64,1},Array{Float64,1}},Float64},Float64,6}} to an object of type Discontinuity{Float64} Closest candidates are: convert(::Type{S}, !Matched::T<:(Union{CategoricalString{R}, CategoricalValue{T,R} where T} where R)) where {S, T<:(Union{CategoricalString{R}, CategoricalValue{T,R} where T} where R)} at /home/asm/.julia/packages/CategoricalArrays/xjesC/src/value.jl:91 convert(::Type{T}, !Matched::T) where T at essentials.jl:154 Discontinuity{Float64}(::Any, !Matched::Any) where tType at /home/asm/.julia/packages/DelayDiffEq/7Ng3E/src/discontinuity_type.jl:8

Stacktrace: [1] fill!(::SubArray{Discontinuity{Float64},1,Array{Discontinuity{Float64},1},Tuple{UnitRange{Int64}},true}, ::Discontinuity{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("#residuals#13")){Int64,Array{Float64,2},Array{Float64,2},Array{Float64,1},Array{Float64,1}},Float64},Float64,6}}) at ./multidimensional.jl:837 [2] __cat(::Array{Discontinuity{Float64},1}, ::Tuple{Int64}, ::Tuple{Bool}, ::Array{Discontinuity{Float64},1}, ::Vararg{Any,N} where N) at ./abstractarray.jl:1412 [3] _cat_t(::Val{1}, ::Type, ::Array{Discontinuity{Float64},1}, ::Vararg{Any,N} where N) at ./abstractarray.jl:1392 [4] #cat_t#103(::Val{1}, ::Function, ::Type{Discontinuity{Float64}}, ::Array{Discontinuity{Float64},1}, ::Vararg{Any,N} where N) at ./abstractarray.jl:1384 [5] (::getfield(Base, Symbol("#kw##cat_t")))(::NamedTuple{(:dims,),Tuple{Val{1}}}, ::typeof(Base.cat_t), ::Type{Discontinuity{Float64}}, ::Array{Discontinuity{Float64},1}, ::Vararg{Any,N} where N) at ./none:0 [6] typed_vcat(::Type, ::Array{Discontinuity{Float64},1}, ::Discontinuity{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("#residuals#13")){Int64,Array{Float64,2},Array{Float64,2},Array{Float64,1},Array{Float64,1}},Float64},Float64,6}}, ::Discontinuity{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("#residuals#13")){Int64,Array{Float64,2},Array{Float64,2},Array{Float64,1},Array{Float64,1}},Float64},Float64,6}}) at ./abstractarray.jl:1493 continues...`

Here is our notebook. Converting the tspan and u0 also did not work, it throws the same error. I would appreciate your advice on this. @ChrisRackauckas Thank you!

[devmotion]

I'll take a look at it and see if I can figure out the issue.

[Farnazmdi]

Thank you very much @devmotion

[devmotion]

Hej!

I had a look at your notebook, and you're right, there was an issue with ForwardDiff if the delays were parameters. Fortunately, the latest changes in DelayDiffEq had already fixed this error. The next release contains this fix, hopefully it will be available in the next days (as soon as https://github.com/JuliaRegistries/General/pull/2324 is merged).

Generally, I would recommend that you have a look at the documentation for parameter estimation. Some comments to your code:

• The second line in

    prob = DDEProblem(DDEmodel, u0, h, tspan, pp; constant_lags = [pp[3], pp[4]])
    _prob = remake(prob; u0=u0, h=h, tspan=tspan)
  

  is not required since all fields of _prob and prob will be equal.
• The line

       solve(_prob, Vern6())
  

  is incorrect - you have to use a DDE algorithm to solve the DDE.
• The first argument of optimize should be a function with scalar output, but residuals (and resid) returns a matrix.
• DDEsolve should compute the solution only at the discrete time points, otherwise the subtractions in resid do not make sense.
• The function resid will fail if the DDE algorithm was not able to compute a complete solution.
• You use a local optimization algorithm, maybe a global optimization algorithm would be better.
• The parameters, and in particular the delays, can become negative but should be constrained to positive numbers. You could also optimize in log space, especially, if the lower and upper bound of possible parameter values differs by multiple orders of magnitude.

I uploaded three example scripts that show how you could perform global optimization with BlackBoxOptim and NLopt and local optimization with Optim. I also added code for automatic differentiation, but it seemed quite unstable in my experiments and I did not manage to compute a solution. I assume that it might be problematic that the delays are also parameters and hence the dual numbers are not only propagated to the states and derivatives, but also to all discontinuities and time points. The global optimization with BlackBoxOptim seemed to work fine, and hence maybe that's something you could try to use.

[ChrisRackauckas]

Have you tried a direct check on the gradient defined by autodiff? We should add a test that it's computing the right thing in that case since it's quite nested.

(But global optimization is probably better anyways here, I think it's still important for us to make sure the autodiff checks out.)

[Farnazmdi]

@devmotion thank you so much, this is very helpful. Greatly appreciated.