ODELocalSensitvityFunction is incompatible with algorithm internal autodiff

[aarmey]

As exemplified in the following code, using ODELocalSensitivityProblem results in errors about missing fields that are present in ODEFunction. This code does run if switched to the ODEProblem counterpart.

using DifferentialEquations

function lorenz!(du, u, p, t)
 du[1] = p[1]*(u[2]-u[1])
 du[2] = u[1]*(p[2]-u[3]) - u[2]
 du[3] = u[1]*u[2] - p[3]*u[3]
end

u0 = [1.0;0.0;0.0]
tspan = (0.0, 100.0)
params = [10.0, 28.0, 8.0/3.0]

probS = ODELocalSensitivityProblem(lorenz!, u0, tspan, params)
prob = ODEProblem(lorenz!, u0, tspan, params)

sol = solve(probS, Rodas4())

print(sol)

Gives:

ERROR: LoadError: type ODELocalSensitvityFunction has no field tgrad
Stacktrace:
 [1] getproperty(::Any, ::Symbol) at ./sysimg.jl:18
 [2] has_tgrad(::ODELocalSensitvityFunction{true,ODEFunction{true,typeof(lorenz!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,Nothing,Nothing,DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(lorenz!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},DiffEqDiffTools.ParamJacobianWrapper{ODEFunction{true,typeof(lorenz!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},Tuple{Array{ForwardDiff.Dual{:___jac_tag,Float64,1},1},Array{ForwardDiff.Dual{:___jac_tag,Float64,1},1}},ForwardDiff.JacobianConfig{ForwardDiff.Tag{DiffEqDiffTools.ParamJacobianWrapper{ODEFunction{true,typeof(lorenz!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},Float64},Float64,3,Tuple{Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.ParamJacobianWrapper{ODEFunction{true,typeof(lorenz!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},Float64},Float64,3},1},Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.ParamJacobianWrapper{ODEFunction{true,typeof(lorenz!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},Float64},Float64,3},1}}},SensitivityAlg{0,true,DataType,true},Array{Float64,1},Nothing,Array{Float64,2},LinearAlgebra.UniformScaling{Bool}}) at /Users/asm/.julia/packages/DiffEqBase/8usQ9/src/diffeqfunction.jl:473
 [3] build_grad_config(::Rodas4{0,true,LinSolveFactorize{typeof(LinearAlgebra.lu!)},DataType}, ::Function, ::Function, ::Array{Float64,1}, ::Float64) at /Users/asm/.julia/packages/OrdinaryDiffEq/7mKtS/src/derivative_wrappers.jl:70
 [4] alg_cache(::Rodas4{0,true,LinSolveFactorize{typeof(LinearAlgebra.lu!)},DataType}, ::Array{Float64,1}, ::Array{Float64,1}, ::Type, ::Type, ::Type, ::Array{Float64,1}, ::Array{Float64,1}, ::ODELocalSensitvityFunction{true,ODEFunction{true,typeof(lorenz!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,Nothing,Nothing,DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(lorenz!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},DiffEqDiffTools.ParamJacobianWrapper{ODEFunction{true,typeof(lorenz!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},Tuple{Array{ForwardDiff.Dual{:___jac_tag,Float64,1},1},Array{ForwardDiff.Dual{:___jac_tag,Float64,1},1}},ForwardDiff.JacobianConfig{ForwardDiff.Tag{DiffEqDiffTools.ParamJacobianWrapper{ODEFunction{true,typeof(lorenz!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},Float64},Float64,3,Tuple{Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.ParamJacobianWrapper{ODEFunction{true,typeof(lorenz!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},Float64},Float64,3},1},Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.ParamJacobianWrapper{ODEFunction{true,typeof(lorenz!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},Float64},Float64,3},1}}},SensitivityAlg{0,true,DataType,true},Array{Float64,1},Nothing,Array{Float64,2},LinearAlgebra.UniformScaling{Bool}}, ::Float64, ::Float64, ::Float64, ::Array{Float64,1}, ::Bool, ::Type{Val{true}}) at /Users/asm/.julia/packages/OrdinaryDiffEq/7mKtS/src/caches/rosenbrock_caches.jl:630
...

My actual problem is stiff and compatible with autodiff, but not ParameterizedFunctions, hence not using the macros available there. Thanks for putting together such an amazing set of tools!

[jlumpe]

Same problem (missing invW) in Julia 0.7 with DifferentialEquations v5.2.1.

[ChrisRackauckas]

Yes, I took a look and the overloads are missing from https://github.com/JuliaDiffEq/DiffEqSensitivity.jl/blob/master/src/local_sensitivity.jl#L1-L36 . This only effects the Rosenbrock methods, so other stiff solvers like KenCarp4 should already work. Since the user-specified time gradient function doesn't cover the sensitivity equations, the easiest fix would be to just set tgrad in the ODELocalSensitvityFunction to be nothing, so that the solvers estimate it using autodiff on their own. That shouldn't be a difficult change if someone has the time to just pass the nothing through an add a test.

Same problem (missing invW) in Julia 0.7 with DifferentialEquations v5.2.1.

v0.7 will have other problems here given how much changed internally with the function representations post v1.0, but I'm not worrying about that since our plan with DiffEq was to have it working enough to upgrade old code but not continue to support it beyond that.

[aarmey]

Great—will take a look at making these changes.

[ChrisRackauckas]

I finished up #36 . There is still an issue when the solvers internally autodiff, so now your code works but only like sol = solve(probS, Rodas4(autodiff=false)), so I'm going to relabel it around that.

This issue can be solved by us defining a jacobian and tgrad for ODELocalSensitivityFunction which internally uses the user's jac and tgrad functions. We could then re-use the Jacobian calculation better (though that's not really much of a speedup since the Jacobian inversion is where the real work happens).

[ChrisRackauckas]

Note that we are updating our sensitivity documentation due to the results of https://arxiv.org/abs/1812.01892 . It is instead recommended to use ForwardDiff.jl for this now, and how to do that via Dual seeding will be explicitly documented in the sensitivity part of the docs (PR for that: SciML/DiffEqDocs.jl#180). Since that ends up being faster and more robust to these kinds of issues, fixing this last bit is a lower priority for me unless there's a major application.