`MethodError` when solving `ODEForwardSensitivityProblem` with `autodiff=true` and `Rosenbrock23()`
topolarity opened this issue · 1 comments
topolarity commented
using LinearAlgebra, OrdinaryDiffEq, SciMLSensitivity
function circuit(du, u, p, t)
L, C = p
du[1] = (-u[1] - u[2] + sin(t)) / C
du[2] = (u[1] - u[2])/L
end
u0 = zeros(2)
tspan = (0, 2pi)
p = [1.0, 1.0]
sprob = ODEForwardSensitivityProblem(circuit, u0, tspan, p)
sol = solve(sprob, Rosenbrock23())gives a MethodError when attempting to use autodiff:
julia> sol = solve(sprob, Rosenbrock23())
ERROR: First call to automatic differentiation for the Jacobian
failed. This means that the user `f` function is not compatible
with automatic differentiation. Methods to fix this include:
1. Turn off automatic differentiation (e.g. Rosenbrock23() becomes
Rosenbrock23(autodiff=false)). More details can befound at
https://docs.sciml.ai/DiffEqDocs/stable/features/performance_overloads/
2. Improving the compatibility of `f` with ForwardDiff.jl automatic
differentiation (using tools like PreallocationTools.jl). More details
can be found at https://docs.sciml.ai/DiffEqDocs/stable/basics/faq/#Autodifferentiation-and-Dual-Numbers
3. Defining analytical Jacobians. More details can be
found at https://docs.sciml.ai/DiffEqDocs/stable/types/ode_types/#SciMLBase.ODEFunction
Note: turning off automatic differentiation tends to have a very minimal
performance impact (for this use case, because it's forward mode for a
square Jacobian. This is different from optimization gradient scenarios).
However, one should be careful as some methods are more sensitive to
accurate gradients than others. Specifically, Rodas methods like `Rodas4`
and `Rodas5P` require accurate Jacobians in order to have good convergence,
while many other methods like BDF (`QNDF`, `FBDF`), SDIRK (`KenCarp4`),
and Rosenbrock-W (`Rosenbrock23`) do not. Thus if using an algorithm which
is sensitive to autodiff and solving at a low tolerance, please change the
algorithm as well.
MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 6})
Closest candidates are:
(::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat
@ Base rounding.jl:266
(::Type{T})(::T) where T<:Number
@ Core boot.jl:893
Float64(::IrrationalConstants.Twoπ)
@ IrrationalConstants ~/.julia/packages/IrrationalConstants/vp5v4/src/macro.jl:112
...
Stacktrace:
[1] jacobian!(J::Matrix{…}, f::SciMLBase.UJacobianWrapper{…}, x::Vector{…}, fx::Vector{…}, integrator::OrdinaryDiffEq.ODEIntegrator{…}, jac_config::SparseDiffTools.ForwardColorJacCache{…})
@ OrdinaryDiffEq ~/repos/OrdinaryDiffEq.jl/src/derivative_wrappers.jl:236
[2] calc_J!(J::Any, integrator::Any, cache::Any, next_step::Bool)
@ OrdinaryDiffEq ~/repos/OrdinaryDiffEq.jl/src/derivative_utils.jl:144 [inlined]
[3] calc_W!
@ ~/repos/OrdinaryDiffEq.jl/src/derivative_utils.jl:703 [inlined]
[4] calc_W!
@ ~/repos/OrdinaryDiffEq.jl/src/derivative_utils.jl:633 [inlined]
[5] calc_rosenbrock_differentiation!(integrator::OrdinaryDiffEq.ODEIntegrator{…}, cache::OrdinaryDiffEq.Rosenbrock23Cache{…}, dtd1::Float64, dtgamma::Float64, repeat_step::Bool, W_transform::Bool)
@ OrdinaryDiffEq ~/repos/OrdinaryDiffEq.jl/src/derivative_utils.jl:796
[6] perform_step!(integrator::OrdinaryDiffEq.ODEIntegrator{…}, cache::OrdinaryDiffEq.Rosenbrock23Cache{…}, repeat_step::Bool)
@ OrdinaryDiffEq ~/repos/OrdinaryDiffEq.jl/src/perform_step/rosenbrock_perform_step.jl:149
[7] perform_step!
@ OrdinaryDiffEq ~/repos/OrdinaryDiffEq.jl/src/perform_step/rosenbrock_perform_step.jl:131 [inlined]
[8] solve!(integrator::OrdinaryDiffEq.ODEIntegrator{…})
@ OrdinaryDiffEq ~/repos/OrdinaryDiffEq.jl/src/solve.jl:538
[9] #__solve#740
@ OrdinaryDiffEq ~/repos/OrdinaryDiffEq.jl/src/solve.jl:6 [inlined]
[10] __solve
@ OrdinaryDiffEq ~/repos/OrdinaryDiffEq.jl/src/solve.jl:1 [inlined]
[11] #solve_call#34
@ DiffEqBase ~/.julia/packages/DiffEqBase/SlYdg/src/solve.jl:561 [inlined]
[12] solve_up(prob::ODEProblem{…}, sensealg::Nothing, u0::Vector{…}, p::Vector{…}, args::Rosenbrock23{…}; kwargs::@Kwargs{…})
@ DiffEqBase ~/.julia/packages/DiffEqBase/SlYdg/src/solve.jl:1010
[13] solve_up
@ ~/.julia/packages/DiffEqBase/SlYdg/src/solve.jl:996 [inlined]
[14] solve(prob::ODEProblem{…}, args::Rosenbrock23{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{…})
@ DiffEqBase ~/.julia/packages/DiffEqBase/SlYdg/src/solve.jl:933
[15] top-level scope
@ REPL[1031]:1
caused by: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 6})
Closest candidates are:
(::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat
@ Base rounding.jl:266
(::Type{T})(::T) where T<:Number
@ Core boot.jl:893
Float64(::IrrationalConstants.Twoπ)
@ IrrationalConstants ~/.julia/packages/IrrationalConstants/vp5v4/src/macro.jl:112
...
Stacktrace:
[1] convert(::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 6})
@ Base ./number.jl:7
[2] setindex!(A::Vector{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 6}, i::Int64)
@ Base ./array.jl:969
[3] copyto_unaliased!(deststyle::IndexStyle, dest::AbstractArray, srcstyle::IndexStyle, src::AbstractArray)
@ Base ./abstractarray.jl:1086 [inlined]
[4] copyto!(dest::Vector{Float64}, src::SubArray{ForwardDiff.Dual{…}, 1, Vector{…}, Tuple{…}, true})
@ Base ./abstractarray.jl:1065
[5] (::ODEForwardSensitivityFunction{…})(du::Vector{…}, u::Vector{…}, p::Vector{…}, t::Float64)
@ SciMLSensitivity ~/.julia/packages/SciMLSensitivity/tUYQm/src/forward_sensitivity.jl:117
[6] (::SciMLBase.UJacobianWrapper{…})(du1::Vector{…}, uprev::Vector{…})
@ SciMLBase ~/.julia/packages/SciMLBase/xja2M/src/function_wrappers.jl:30 [inlined]
[7] forwarddiff_color_jacobian!(J::Matrix{…}, f::SciMLBase.UJacobianWrapper{…}, x::Vector{…}, jac_cache::SparseDiffTools.ForwardColorJacCache{…})
@ SparseDiffTools ~/.julia/packages/SparseDiffTools/yhwO4/src/differentiation/compute_jacobian_ad.jl:371
[8] jacobian!(J::Matrix{…}, f::SciMLBase.UJacobianWrapper{…}, x::Vector{…}, fx::Vector{…}, integrator::OrdinaryDiffEq.ODEIntegrator{…}, jac_config::SparseDiffTools.ForwardColorJacCache{…})
@ OrdinaryDiffEq ~/repos/OrdinaryDiffEq.jl/src/derivative_wrappers.jl:234
[9] calc_J!(J::Any, integrator::Any, cache::Any, next_step::Bool)
@ OrdinaryDiffEq ~/repos/OrdinaryDiffEq.jl/src/derivative_utils.jl:144 [inlined]
[10] calc_W!
@ ~/repos/OrdinaryDiffEq.jl/src/derivative_utils.jl:703 [inlined]
[11] calc_W!
@ ~/repos/OrdinaryDiffEq.jl/src/derivative_utils.jl:633 [inlined]
[12] calc_rosenbrock_differentiation!(integrator::OrdinaryDiffEq.ODEIntegrator{…}, cache::OrdinaryDiffEq.Rosenbrock23Cache{…}, dtd1::Float64, dtgamma::Float64, repeat_step::Bool, W_transform::Bool)
@ OrdinaryDiffEq ~/repos/OrdinaryDiffEq.jl/src/derivative_utils.jl:796
[13] perform_step!(integrator::OrdinaryDiffEq.ODEIntegrator{…}, cache::OrdinaryDiffEq.Rosenbrock23Cache{…}, repeat_step::Bool)
@ OrdinaryDiffEq ~/repos/OrdinaryDiffEq.jl/src/perform_step/rosenbrock_perform_step.jl:149
[14] perform_step!
@ OrdinaryDiffEq ~/repos/OrdinaryDiffEq.jl/src/perform_step/rosenbrock_perform_step.jl:131 [inlined]
[15] solve!(integrator::OrdinaryDiffEq.ODEIntegrator{…})
@ OrdinaryDiffEq ~/repos/OrdinaryDiffEq.jl/src/solve.jl:538
[16] #__solve#740
@ OrdinaryDiffEq ~/repos/OrdinaryDiffEq.jl/src/solve.jl:6 [inlined]
[17] __solve
@ OrdinaryDiffEq ~/repos/OrdinaryDiffEq.jl/src/solve.jl:1 [inlined]
[18] #solve_call#34
@ DiffEqBase ~/.julia/packages/DiffEqBase/SlYdg/src/solve.jl:561 [inlined]
[19] solve_up(prob::ODEProblem{…}, sensealg::Nothing, u0::Vector{…}, p::Vector{…}, args::Rosenbrock23{…}; kwargs::@Kwargs{…})
@ DiffEqBase ~/.julia/packages/DiffEqBase/SlYdg/src/solve.jl:1010
[20] solve_up
@ ~/.julia/packages/DiffEqBase/SlYdg/src/solve.jl:996 [inlined]
[21] solve(prob::ODEProblem{…}, args::Rosenbrock23{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{…})
@ DiffEqBase ~/.julia/packages/DiffEqBase/SlYdg/src/solve.jl:933
[22] top-level scope
@ REPL[1031]:1
Some type information was truncated. Use `show(err)` to see complete types.gabo-di commented
Like the output suggests, using autodiff=false does the trick,
using LinearAlgebra, OrdinaryDiffEq, SciMLSensitivity
function circuit(du, u, p, t)
L, C = p
du[1] = (-u[1] - u[2] + sin(t)) / C
du[2] = (u[1] - u[2])/L
end
u0 = zeros(2)
tspan = (0.0, 2*pi)
p = [1.0, 1.0]
sprob = ODEForwardSensitivityProblem(circuit, p, tspan, p)
sol = solve(sprob, Rosenbrock23(autodiff=false))I think it is vaguely related to this issue #35 .
Note that I need to use tspan = (0.0, 2pi) instead of tspan = (0, 2pi)