Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

MethodError when solving ODEForwardSensitivityProblem with autodiff=true and Rosenbrock23() #944

Open
topolarity opened this issue Nov 30, 2023 · 1 comment
Labels

Comments

@topolarity
Copy link

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.
@topolarity topolarity added the bug label Nov 30, 2023
@topolarity topolarity changed the title MethodError when solving ODEForwardSensitivityProblem with autodiff=true and Rosenbrock23()` MethodError when solving ODEForwardSensitivityProblem with autodiff=true and Rosenbrock23() Nov 30, 2023
@gabo-di
Copy link

gabo-di commented Jan 4, 2024

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)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

No branches or pull requests

2 participants