Closed hersle closed 4 months ago
Doing it through the whole symbolic system is a bit harder, but it should be possible. You most likely don't want to be doing this though for performance reasons: this would cause extra compilation in the optimization loop. But if you're doing something like reinforcement learning and trying new equations then you might have some interesting applications of this.
Thanks! In case it's helpful, here are some simplified details about my real use case:
I have an ODE model with, say, 3 variables constrained by 1 equation. I want the user to specify initial conditions for any 2 independent variables, and have MTK solve for the third dependent variable. I want all this to work with automatic differentiation.
Here is an example that works with FiniteDiff, but fails with ForwardDiff:
using Test
using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using DifferentialEquations
using ForwardDiff, FiniteDiff
# ODE with 3 independent variables and 1 constraint
@variables x(t) y(t) z(t) s(t)
sys = structural_simplify(ODESystem([
D(x) ~ 0,
D(y) ~ 0,
D(z) ~ 0,
s ~ x + y + z
], t, [x, y, z, s], []; name=:sys))
function solve_instance(ics)
prob = ODEProblem(sys, [s => 1.0; ics], (0.0, 1.0)) # s0 = x0 + y0 + z0 = 1
sol = solve(prob)
return sol
end
# user solves with initial conditions for different subsets of variables
x0(y0, z0) = solve_instance([y => y0, z => z0])(0.0, idxs=x) # analytical: x0 = 1 - y0 - z0
y0(x0, z0) = solve_instance([x => x0, z => z0])(0.0, idxs=y) # analytical: y0 = 1 - x0 - z0
z0(y0, x0) = solve_instance([x => x0, y => y0])(0.0, idxs=z) # analytical: z0 = 1 - x0 - y0
# evaluation passes
@test x0(0.1, 0.2) ≈ y0(0.1, 0.2) ≈ z0(0.1, 0.2) ≈ 0.7
# gradient passes with FiniteDiff, but fails with ForwardDiff
x0(ics::AbstractArray) = x0(ics...)
y0(ics::AbstractArray) = y0(ics...)
z0(ics::AbstractArray) = z0(ics...)
for ∇ in [FiniteDiff.finite_difference_gradient, ForwardDiff.gradient]
@test all(∇(x0, [0.1, 0.2]) .≈ ∇(y0, [0.1, 0.2]) .≈ ∇(z0, [0.1, 0.2]) .≈ [-1, -1])
end
This fails with a very similar error (and the same root cause?):
Test threw exception
Expression: all(∇(x0, [0.1, 0.2]) .≈ ∇(y0, [0.1, 0.2]) .≈ ∇(z0, [0.1, 0.2]) .≈ [-1, -1])
MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{typeof(y0), Float64}, Float64, 2})
Closest candidates are:
(::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat
@ Base rounding.jl:207
(::Type{T})(::T) where T<:Number
@ Core boot.jl:792
Float64(::IrrationalConstants.Loghalf)
@ IrrationalConstants C:\Users\herma\.julia\packages\IrrationalConstants\vp5v4\src\macro.jl:112
...
Stacktrace:
[1] convert(::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{typeof(y0), Float64}, Float64, 2})
@ Base .\number.jl:7
[2] symconvert(::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{typeof(y0), Float64}, Float64, 2})
@ ModelingToolkit C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\parameter_buffer.jl:2
[3] ModelingToolkit.MTKParameters(sys::NonlinearSystem, p::Vector{Pair{SymbolicUtils.BasicSymbolic{Real}, Float64}}, u0::Dict{Any, Any}; tofloat::Bool, use_union::Bool)
@ ModelingToolkit C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\parameter_buffer.jl:114
[4] MTKParameters
@ C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\parameter_buffer.jl:13 [inlined]
[5] process_NonlinearProblem(constructor::Type, sys::NonlinearSystem, u0map::Dict{Any, Any}, parammap::Vector{Pair{SymbolicUtils.BasicSymbolic{Real}, Float64}}; version::Nothing, jac::Bool, checkbounds::Bool, sparse::Bool, simplify::Bool, linenumbers::Bool, parallel::Symbolics.SerialForm, eval_expression::Bool, use_union::Bool, tofloat::Bool, kwargs::@Kwargs{check_length::Bool})
@ ModelingToolkit C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\nonlinear\nonlinearsystem.jl:400
[6] process_NonlinearProblem
@ C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\nonlinear\nonlinearsystem.jl:384 [inlined]
[7] (NonlinearProblem{true})(sys::NonlinearSystem, u0map::Dict{Any, Any}, parammap::Vector{Pair{SymbolicUtils.BasicSymbolic{Real}, Float64}}; check_length::Bool, kwargs::@Kwargs{})
@ ModelingToolkit C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\nonlinear\nonlinearsystem.jl:435
[8] (NonlinearProblem{true})(sys::NonlinearSystem, u0map::Dict{Any, Any}, parammap::Vector{Pair{SymbolicUtils.BasicSymbolic{Real}, Float64}})
@ ModelingToolkit C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\nonlinear\nonlinearsystem.jl:429
[9] NonlinearProblem(::NonlinearSystem, ::Dict{Any, Any}, ::Vararg{Any}; kwargs::@Kwargs{})
@ ModelingToolkit C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\nonlinear\nonlinearsystem.jl:426
[10] NonlinearProblem(::NonlinearSystem, ::Dict{Any, Any}, ::Vararg{Any})
@ ModelingToolkit C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\nonlinear\nonlinearsystem.jl:425
[11] ModelingToolkit.InitializationProblem{true, SciMLBase.AutoSpecialize}(sys::ODESystem, t::Float64, u0map::Vector{Pair{Num, ForwardDiff.Dual{ForwardDiff.Tag{typeof(y0), Float64}, Float64, 2}}}, parammap::SciMLBase.NullParameters; guesses::Dict{Any, Any}, check_length::Bool, warn_initialize_determined::Bool, kwargs::@Kwargs{})
@ ModelingToolkit C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1636
[12] (ModelingToolkit.InitializationProblem{true})(::ODESystem, ::Float64, ::Vararg{Any}; kwargs::@Kwargs{guesses::Dict{Any, Any}, warn_initialize_determined::Bool})
@ ModelingToolkit C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1572
[13] ModelingToolkit.InitializationProblem(::ODESystem, ::Float64, ::Vararg{Any}; kwargs::@Kwargs{guesses::Dict{Any, Any}, warn_initialize_determined::Bool})
@ ModelingToolkit C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1560
[14] process_DEProblem(constructor::Type, sys::ODESystem, u0map::Vector{Pair{Num, ForwardDiff.Dual{ForwardDiff.Tag{typeof(y0), Float64}, Float64, 2}}}, parammap::SciMLBase.NullParameters; implicit_dae::Bool, du0map::Nothing, version::Nothing, tgrad::Bool, jac::Bool, checkbounds::Bool, sparse::Bool, simplify::Bool, linenumbers::Bool, parallel::Symbolics.SerialForm, eval_expression::Bool, use_union::Bool, tofloat::Bool, symbolic_u0::Bool, u0_constructor::typeof(identity), guesses::Dict{Any, Any}, t::Float64, warn_initialize_determined::Bool, build_initializeprob::Bool, kwargs::@Kwargs{check_length::Bool})
@ ModelingToolkit C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:921
[15] process_DEProblem
@ C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:834 [inlined]
[16] (ODEProblem{true, SciMLBase.AutoSpecialize})(sys::ODESystem, u0map::Vector{Pair{Num, ForwardDiff.Dual{ForwardDiff.Tag{typeof(y0), Float64}, Float64, 2}}}, tspan::Tuple{Float64, Float64}, parammap::SciMLBase.NullParameters; callback::Nothing, check_length::Bool, warn_initialize_determined::Bool, kwargs::@Kwargs{})
@ ModelingToolkit C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1085
[17] ODEProblem
@ C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1075 [inlined]
[18] (ODEProblem{true, SciMLBase.AutoSpecialize})(sys::ODESystem, u0map::Vector{Pair{Num, ForwardDiff.Dual{ForwardDiff.Tag{typeof(y0), Float64}, Float64, 2}}}, tspan::Tuple{Float64, Float64})
@ ModelingToolkit C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1075
[19] (ODEProblem{true})(::ODESystem, ::Vector{Pair{Num, ForwardDiff.Dual{ForwardDiff.Tag{typeof(y0), Float64}, Float64, 2}}}, ::Vararg{Any}; kwargs::@Kwargs{})
@ ModelingToolkit C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1062
[20] (ODEProblem{true})(::ODESystem, ::Vector{Pair{Num, ForwardDiff.Dual{ForwardDiff.Tag{typeof(y0), Float64}, Float64, 2}}}, ::Vararg{Any})
@ ModelingToolkit C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1061
[21] ODEProblem(::ODESystem, ::Vector{Pair{Num, ForwardDiff.Dual{ForwardDiff.Tag{typeof(y0), Float64}, Float64, 2}}}, ::Vararg{Any}; kwargs::@Kwargs{})
@ ModelingToolkit C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1051
[22] ODEProblem(::ODESystem, ::Vector{Pair{Num, ForwardDiff.Dual{ForwardDiff.Tag{typeof(y0), Float64}, Float64, 2}}}, ::Vararg{Any})
@ ModelingToolkit C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1050
[23] solve_instance(ics::Vector{Pair{Num, ForwardDiff.Dual{ForwardDiff.Tag{typeof(y0), Float64}, Float64, 2}}})
@ Main .\REPL[9]:2
[24] y0(x0::ForwardDiff.Dual{ForwardDiff.Tag{typeof(y0), Float64}, Float64, 2}, z0::ForwardDiff.Dual{ForwardDiff.Tag{typeof(y0), Float64}, Float64, 2})
@ Main .\REPL[12]:1
[25] y0
@ .\REPL[18]:1 [inlined]
[26] vector_mode_dual_eval!
@ C:\Users\herma\.julia\packages\ForwardDiff\PcZ48\src\apiutils.jl:24 [inlined]
[27] vector_mode_gradient(f::typeof(y0), x::Vector{Float64}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(y0), Float64}, Float64, 2, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(y0), Float64}, Float64, 2}}})
@ ForwardDiff C:\Users\herma\.julia\packages\ForwardDiff\PcZ48\src\gradient.jl:89
...
Would it be more efficient to solve the initialization problem rather than create a new ODE problem? For example:
# continued from above code snippet
x0init(y0, z0) = solve(ModelingToolkit.InitializationProblem(sys, 0.0, [s => 1.0, y => y0, z => z0], []))[x]
@test x0init(0.1, 0.2) ≈ 0.7 # passes
x0init(ics::AbstractArray) = x0init(ics...)
for ∇ in [FiniteDiff.finite_difference_gradient, ForwardDiff.gradient]
@test all(∇(x0init, [0.1, 0.2]) .≈ [-1, -1]) # fails with ForwardDiff
end
This also fails with ForwardDiff with the (seemingly) same error.
Having all this work seamlessly would be absolutely fantastic 😃
Auto-differentiating through
remake(ODEProblem())
works, but not directly throughODEProblem()
:The last test errors with
Could it be unified to work in both ways?
I'm on a fresh updated
master
branch with] status