search

SciML / SciMLSensitivity.jl

A component of the DiffEq ecosystem for enabling sensitivity analysis for scientific machine learning (SciML). Optimize-then-discretize, discretize-then-optimize, adjoint methods, and more for ODEs, SDEs, DDEs, DAEs, etc.
https://docs.sciml.ai/SciMLSensitivity/stable/
Other
332 stars 70 forks source link

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

Open topolarity opened 11 months ago

topolarity commented 11 months ago 
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 10 months ago

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)